O Galland

J Scheibert

Analytical model of surface uplift above axisymmetric flat-lying magma intrusions: Implications for sill emplacement and geodesy

Keywords: Sills, Laccoliths, Analytical model, Axi-symmetric, Ground deformation, Mogi

In this paper, we develop a new axisymmetric analytic model of surface uplift upon sills and laccoliths, based on the formulation of a thin bending plate lying on an elastic foundation. In contrast to most former models also based on thin bending plate formulation, our model accounts for (i) axi-symmetrical uplift, (ii) both upon and outside the intrusion. The model accounts for shallow intrusions, i.e. the ratio a/h > 5 where a and h are the radius and depth of the intrusion, respectively. The main parameter of the model is the elastic length l, which is a function of the elastic properties of the bending plate and of the elastic foundation. The model exhibits two regimes depending on the ratio a/l. When a/l < 5, the uplift spreads over a substantial domain compared to that of the intrusion. In contrast, when a/l > 5, the uplift is mostly restricted upon the intrusion.

When the elastic foundation is very stiff, our model converges towards that of a clamped plate. We provide, as supplementary material, a Matlab function that calculates the surface uplift from the set of system and control parameters. We discuss three possible applications of our model: (i) The model can be used to describe sill propagation by introducing a propagation criterion. For realistic values, our model reproduces well the behavior of horizontal intrusions simulated in experiments; (ii) The model can also be used to compute the critical size of saucer-shaped sills. It shows, for instance, that a soft elastic foundation favors the horizontal spreading of sills before they form inclined sheets; (iii) We show that the classical Mogi point source model cannot be used to con-1 strain sill properties from the surface uplift. We thus propose that our model can be used as a valuable alternative to both simple analytical models like Mogi's and more complex numerical models used to analyze ground deformation resulting from sill intrusions in active volcanoes.

Introduction

Surface deformation in active volcanic systems is generally assumed to reflect the dynamics of magma intrusion and transport at depth. Modern monitoring techniques allow good measurements of the deformation of volcanic edifices before, during and after an eruption [START_REF] Cayol | Three-dimensional modeling of the 1983-1984 eruption at Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Froger | Active spreading and regional extension at Mount Etna imaged by SAR interferometry[END_REF][START_REF] Froger | Two scales of inflation at Lastarria-Cordon del Azufre volcanic complex, central Andes, revealed from ASAR-ENVISAT interferometric data[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Galland | Experimental modelling of ground deformation associated with shallow magma intrusions[END_REF]. A good example is the Eyjafjallajökull Volcano, Iceland, which has been monitored for more than a decade using InSAR and GPS data (Fig. 1a).

This data allowed to detect the onset of its unrest [START_REF] Pedersen | Temporal development of the 1999 intrusive episode in the Eyjafjallajökull volcano, Iceland, derived from InSAR images[END_REF], and overall the premises of the 2010 eruption that caused massive disruptions in the air traffic in Europe [START_REF] Sigmundsson | Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption[END_REF].

The mechanical analysis of surface deformation patterns is commonly used to constrain the geometry and dynamics of the magma plumbing systems a posteriori [START_REF] Cayol | Three-dimensional modeling of the 1983-1984 eruption at Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric syntheric aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Masterlark | Magma intrusion and deformation predictions: Sensitivities to the Mogi assumptions[END_REF]. One of the first attempts in analyzing surface deformation on active volcanoes has been performed by Mogi [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surface around them[END_REF], who developed an analytical solution of surface deformation induced by a small spherical over-pressured magma reservoir. The solution of Mogi is valid when the size of the magma reservoir is very small compared to its depth, i.e. when a m /h m << 1, where a m and h m are the radius and the depth of the center of the reservoir, respectively.

Although the so-called "Mogi point source" solution provides good fits with data monitored on some active volcanoes, recent studies show that (1) many magma reservoirs do not consist of spherical chambers but exhibit a flat-lying shape [e.g. [START_REF] Amelung | Widespread uplift and "trapdoor" faulting on Galapagos volcanoes observed with radar interferometry[END_REF][START_REF] Fialko | Deformation due to a pressurized horizontal circular crack in an elastic half space, with applications to volcano geodesy[END_REF][START_REF] Fialko | Finite source modelling of magmatic unrest in Socorro, New Mexico, and Long Valley, California[END_REF][START_REF] Chang | Accelerated uplift and magmatic intrusion of the Yellowstone caldera, 2004 to 2006[END_REF][START_REF] Sigmundsson | Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption[END_REF][START_REF] Woo | Intrusion and deformation at Campi Flegrei, southern Italy: Sills, dikes, and regional extension[END_REF][START_REF] Pedersen | Temporal development of the 1999 intrusive episode in the Eyjafjallajökull volcano, Iceland, derived from InSAR images[END_REF][START_REF] Sigmundsson | Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption[END_REF], and (2) the roof of the reservoir can be very shallow (<3 km) [START_REF] Mctigue | Elastic stress and deformation near a finite spherical magma body: Resolution of the point source paradox[END_REF][START_REF] Brandsdóttir | Thin low-velocity zone within the Krafla caldera, NE-Iceland attributed to a small magma chamber[END_REF][START_REF] Gudmundsson | The crustal magma chamber of the Katla volcano in south Iceland revealed by 2-D seismic undershooting[END_REF][START_REF] Brandsdóttir | Färoe-Iceland Ridge Experiment 2. Crustal structure of the Krafla central volcano[END_REF][START_REF] Brandsdóttir | The seismic structure of Iceland[END_REF][START_REF] Elders | Origin of a rhyolite that intruded a geothermal well while drilling at the Krafla volcano, Iceland[END_REF]],

1 olivier.galland@fys.uio.no, tel:+4722856719, Fax: +4722855101 such that the assumption a m /h m << 1 is not satisfied. The Mogi point source solution is thus not applicable in these conditions, and a more relevant analytical solution is needed to interpret surface deformation data. Surface deformation has not only been observed in active volcanoes as a passive consequence of shallow magma emplacement, but it can exert an active mechanical feedback on the emplacement of laccoliths and saucer-shaped sills [e.g. [START_REF] Gilbert | Report on the geology of the Henry Mountains[END_REF][START_REF] Jackson | Flexure and faulting of sedimentary host rocks during growth of igneous domes, Henry Mountains, Utah[END_REF][START_REF] Jackson | Processes of laccolithic emplacement in the Southern Henry Mountains, Southeastern Utah[END_REF][START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF][START_REF] Polteau | Saucer-shaped intrusions: occurrences, emplacement and implications[END_REF][START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF][START_REF] Galerne | 3d relationships between sills and their feeders: evidence from the Golden Valley Sill Complex (Karoo Basin) and experimental modelling[END_REF]. On seismic images and in the field, it can be observed that saucer-shaped sills, for instance, are closely associated with uplift and bending of their overlying strata, the inclined sheets being located under the edges of the uplifted area [Fig. 1b;[START_REF] Hansen | Saucer-shaped sill with lobate morphology revealed by 3D seismic data: implications for resolving a shallow-level sill emplacement mechanism[END_REF][START_REF] Muirhead | Interconnected sills and inclined sheet intrusions control shallow magma transport in the ferrar large igneous province, antarctica[END_REF]. This relationship has been interpreted as a result of the mechanical interaction between the bending of the overlying strata and the spreading of the sills: the differential uplift at the edges of the domes generates stresses that interact with the leading edges of the sills, which in turn are deflected towards the surface [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF][START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF][START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF]. This mechanism producing inclined sheets substantially contributes to magma ascent through sedimentary basins [START_REF] Cartwright | Magma transport through the crust via interconnected sill complexes[END_REF][START_REF] Muirhead | Interconnected sills and inclined sheet intrusions control shallow magma transport in the ferrar large igneous province, antarctica[END_REF] In order to [START_REF] Abdelmalak | Fracture mode analysis and related surface deformation during dyke intrusion: Results from 2D experimental modelling[END_REF] better predict surface deformation due to the emplacement of shallow flat-lying intrusions and (2) better quantify how the bending of strata affects the emplacement of laccoliths and saucer-shaped sills, one needs to better constrain the mechanics of surface deformation. In this paper, we develop a new analytical model of surface deformation above shallow axially symmetric flat-lying intrusions. Our model is based on the theory of a thin bending plate lying on a deformable elastic foundation. After the theoretical development, we discuss the limitations and the effects of the parameters on the model. Subsequently, we discuss some geological applications, notably for sill propagation and for saucer-shaped sills. We also briefly discuss the potential application of our model to the analysis of ground deformation due to sill intrusion in active volcanoes. We provide the code of our model as a Matlab function available as supplementary material.

Existing solutions

Several models have been developed to calculate the deformation field associated with flat-lying intrusions. Because sills are sheet intrusions, they can be considered as horizontal fluid-filled cracks within an elastic medium. The strains in the elastic medium are considered to be small everywhere except close to the tip of the sill, and the formulation of the problem can be achieved using the theory of linear elasticity.

A first approach attempted to develop models tending toward a complete description of cracks in an elastic half-space. [START_REF] Sun | Theoretical size of hydraulically induced horizontal fractures and corresponding surface uplift in an idealized medium[END_REF] described the deformation of a free surface above a fluid-pressurized crack by developing a 3D approximate solution of the vertical and horizontal displacements above a circular crack. To do that, he superimposed (i) the solution for the displacements due to a crack in an infinite elastic medium and (ii) an auxiliary stress function that satisfies the zero traction boundary conditions at the free surface.

However, such superposition generates significant errors when the crack's radius-to-depth ratio a/h becomes greater than 1 [START_REF] Fialko | Finite source modelling of magmatic unrest in Socorro, New Mexico, and Long Valley, California[END_REF]. [START_REF] Pollard | On the mechanical interaction between a fluid-filled fracture and the Earth's surface[END_REF] developed the 2D equations that account for the surface deformation and the stress intensity factor at the tip of an arbitrarily oriented crack contained in an elastic half-space. More recently, [START_REF] Fialko | Evidence for on-going inflation of the Socorro magma body, New Mexico, from interferometric synthetic aperture radar imaging[END_REF] derived the axi-symmetrical solutions for the stress and displacements associated with a horizontal circular crack in an elastic half-space. Although these approaches are powerful tools, obtaining the solutions requires numerical integrations, which can be tricky to implement, like for most analyses of elastic layers under stress [see e.g. 74].

A second, classical approach for describing surface uplift due to a sill or a laccolith is based on thin elastic plate theory [START_REF] Timoshenko | Theory of plates and shells[END_REF]. This theory can easily be adapted to investigate the deflection of sedimentary strata above a magma-filled horizontal sill or laccolith [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Scaillet | The Gangotri granite (Garhwal Himalaya): Laccolithic emplacement in an extending collisional belt[END_REF][START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF][START_REF] Michaut | Dynamics of magmatic intrusions in the upper crust: Theory and applications to laccoliths on Earth and the Moon[END_REF]. This theory accounts only for shallow sills, i.e. the radius of the crack a is large with respect to its depth h (typically a/h > 5), which is the case for many sills. This approach, developed both in 2D and 3D, has been extensively used and generally accepted, such that it is presented as a classical model in textbooks [START_REF] Turcotte | Geodynamics[END_REF] and is used in the mechanics community [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 2. Propagation[END_REF][START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 1. Methods and observations[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF][START_REF] Bunger | Near-surface hydraulic fracture[END_REF][START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF].

The application of thin plate theory alone to sills and laccoliths takes into account the weight of the overburden (q 0 ), heterogeneous magma pressure distributions in the intrusion (P (x)), and the mechanical layering of the overburden [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF]. Very recently, the thin plate formulation has been coupled with the equations for viscous fluid flow to model the emplacement of viscous magma into sills and laccoliths [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF][START_REF] Michaut | Dynamics of magmatic intrusions in the upper crust: Theory and applications to laccoliths on Earth and the Moon[END_REF]. In addition, it has been extended to derive a simple criterion for the upward propagation of saucershaped sills [START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF]. Nevertheless, the formulations developed in these papers assume that the bending plate is clamped to a rigid foundation at the tips of the intrusions [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF]. In other words, the bending of the overburden is restricted upon the intrusion only, which is not realistic in geological systems: for instance Figure 1a shows that the uplift measured at the Eyjafjallajökull volcano, resulting from the emplacement of a sill, is significant in a larger area than that of the sill [Fig. 1a; 77]. Similarly, Figure 1b shows a seismic profile that images a sill and the associated doming of its overburden in the Rockall Basin, offshore Scotland [START_REF] Hansen | The three-dimensional geometry and growth of forced folds above saucer-shaped igneous sills[END_REF]: the doming of the overburden is larger than the associated sill.

A more realistic formulation of surface deformation above laccoliths has been proposed by [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF], based on the theory of a thin bending plate lying on an elastic foundation [START_REF] Timoshenko | Theory of plates and shells[END_REF] (p. 259-269). This 2D model considers the continuous bending of a plate across two domains, one above the intrusion and one around the intrusion. Above the intrusion, the formulation of the bending plate is the same as that with no elastic foundation. In contrast, outside the intrusion the plate is attached to a deformable elastic foundation, and the plate can also deform. Thus, the bending of the plate is not only restricted above the intrusion, and displacement of the plate at the tip of the intrusion is allowed.

In this paper, we improve the 2D solution of [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF] in several ways: we derive an axi-symmetrical analytical solution of surface uplift above a sill or laccolith and we use a generic set of boundary conditions. Note that we consider possible heterogeneous distributions of internal pressure into the sill, as [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF], but we do not restrict our analysis to vanishing internal pressure at the sill periphery. Like [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] and [START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF], the magmatic pressure P 0 and the radius a of the intrusion are known a priori and serve as control parameters. This general formulation allows to characterize, for any values of the input parameters, the instantaneous state of an intrusion in the mechanical equilibrium, i.e. without considering propagation. A propagating intrusion corresponds instead to a quasistatic mechanical system, which requires the definition of a propagation criterion. This will be developed as an application of our model in section 5.1.

Theoretical development

The model

We consider the system sketched in Figure 2: an axisymmetric flat intrusion of radius a lying under a linear elastic strata of thickness h, Young modulus E, Poisson ratio ν and mass density ρ (see Table 1 for definitions of the model parameters). We assume that the intrusion is shallow (a/h > 5), so that the strata can be considered as a thin plate with a bending stiffness D = Eh 3 12(1-ν 2 ) . Above the intrusion (radial distance r < a), the plate is submitted to a radial pressure profile of the form P = P 0 -(P 0 -P a )(r/a) n , in which P 0 and P a are the pressure values at the center (r=0) and periphery (r=a) of the intrusion, respectively (Fig. 2d). We consider a heterogeneous pressure to account for the viscous drag due to the flow of viscous magma into the intrusion. Note that a homogeneous pressure is obtained for n = ∞ and/or P a =P 0 . Outside the intrusion (r > a), the plate is attached to an elastic foundation of elastic modulus k. Such an elastic foundation can accommodate substantial elastic deformation, allowing for displacement outside the intrusion without propagation of the intrusion tip. At all points of the model, the strata is also submitted to the lithostatic stress q 0 = ρgh. Note that below the intrusion and the elastic foundation, the basement is considered to be infinitely rigid, like in the analyses of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF], [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF], [START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF] and [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF]. This assumption is validated by geological observations [START_REF] Gilbert | Report on the geology of the Henry Mountains[END_REF][START_REF] Jackson | Processes of laccolithic emplacement in the Southern Henry Mountains, Southeastern Utah[END_REF][START_REF] Johnson | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, i. Field observations, Gilbert's model, physical properties and flow of the magma[END_REF], geophysical observations [START_REF] Hansen | The three-dimensional geometry and growth of forced folds above saucer-shaped igneous sills[END_REF][START_REF] Polteau | Saucer-shaped intrusions: occurrences, emplacement and implications[END_REF] and modelling results [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF], which show that deformation associated with shallow flat-lying intrusions mostly affect the overburden, but not the substratum.

From thin plate theory, we can write the equilibrium equations of the system as:

D∆ 2 w = q 0 -P 0 + (P 0 -P a )(r/a) n , 0 < r < a, (1)
D∆ 2 w + kw = q 0 , r > a, (2)
where w is the vertical displacement of the plate and ∆ 2 is the bilaplacian operator.

In the following sections, we will refer to w 1 and w 2 for the displacements upon (0 < r < a) and outside (a < r) the sill, respectively. Equation ( 1), when taken in axisymmetric form with abscissa r, has a general solution of the form [see 85, page 54, equation 60]:

w 1 = (q 0 -P 0 )r 4 64D + C 1 r 2 4 + C 2 + C 5 log r a + (P 0 -P a )r n+4 Da n (n + 2) 2 (n + 4) 2 (3)
We set C 5 = 0 because the logarithm would lead to a displacement singularity at r = 0.

We are left with only two unknown constants C 1 and C 2 .

The general solution of Equation (2) when the right member is 0, and when taken in axisymmetric form, writes [see 85, p266, equation h]:

w 2 = C 3 kei 0 (x) + C 4 ker 0 (x) + C 6 ber 0 (x) + C 7 bei 0 (x) (4)
with x = r l , l = 4 D k , and ber ν , bei ν , ker ν , kei ν are Kelvin functions [START_REF] Timoshenko | Theory of plates and shells[END_REF]. We can set C 6 and C 7 to 0 because lim r→∞ ber 0 (r) = ∞ and lim r→∞ bei 0 (r) = ∞, which would yield unphysical infinite displacements far from the sill. Equation (2) also has a constant solution, w 0 = q 0 /k, which must be added to Equation ( 4) to obtain the complete solution. Note that adding this term provides similar boundary condition at r → ∞ as that of [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF], i.e. the displacement w 0 = lim r→∞ w 2 is not zero and corresponds to the effect of the weight of the plate on the elastic foundation.

We are left with the following two equations, with C 1 , C 2 , C 3 and C 4 being four unknown constants:

w 1 = (q 0 -P 0 )r 4 64D + C 1 r 2 4 + C 2 + (P 0 -P a )r n+4 Da n (n + 2) 2 (n + 4) 2 0 < r < a (5)
w 2 = C 3 kei 0 r l + C 4 ker 0 r l + q 0 k r > a (6)
We therefore need four boundary conditions to solve the mathematical problem. Continuity of the displacement w and its three first derivatives with respect to r at r = a yield these four boundary conditions, which write:

w 1 (a) = w 2 (a), (7)
w 1 (a) = w 2 (a), (8)
w 1 (a) = w 2 (a), (9)
w 1 (a) = w 2 (a), (10)
with prime standing for derivation with respect to r. This system of four equations allows for the determination of the four constants.

Equation [START_REF] Bradley | Intrusion of major dolerite sills[END_REF] writes:

(q 0 -P 0 )a 4 64D + C 1 a 2 4 + C 2 + (P 0 -P a )a 4 D(n + 2) 2 (n + 4) 2 = C 3 kei 0 a l + C 4 ker 0 a l + q 0 k (11)
Equation ( 8) writes:

(q 0 -P 0 )a 3 16D + C 1 a 2 + (P 0 -P a )a 3 D(n + 2) 2 (n + 4) = C 3 √ 2l kei 1 a l -ker 1 a l + C 4 √ 2l kei 1 a l + ker 1 a l (12)
Equation ( 9) writes:

3(q 0 -P 0 )a 2 16D + C 1 2 + (P 0 -P a )a 2 (n + 3) D(n + 2) 2 (n + 4) = C 3 2l 2 ker 0 a l -ker 2 a l + C 4 2l 2 kei 2 a l -kei 0 a l (13)
Equation [START_REF] Brandsdóttir | The seismic structure of Iceland[END_REF] writes:

3(q 0 -P 0 )a 8D + (P 0 -P a )a(n + 3) D(n + 2)(n + 4) = C 3 4 √ 2l 3 3ker 1 a l -ker 3 a l + 3kei 1 a l -kei 3 a l + C 4 4 √ 2l 3 3ker 1 a l -ker 3 a l -3kei 1 a l + kei 3 a l (14)
These equations constitute a system of four coupled linear equations, which can be written matricially as :

A.C = B, (15)
with

A =          a 2 4 1 kei 0 ( a l ) ker 0 ( a l ) a 2 0 kei1( a l )-ker1( a l ) √ 2l kei1( a l )+ker1( a l ) √ 2l 1 2 0 ker0( a l )-ker2( a l ) 2l 2 kei2( a l )-kei0( a l ) 2l 2 0 0 3ker1( a l )-ker3( a l )+3kei1( a l )-kei3( a l ) 4 √ 2l 3 3ker1( a l )-ker3( a l )-3kei1( a l )+kei3( a l ) 4 √ 2l 3          , B =          q0 k -(q0-P0)a 4 64D -(P0-Pa)a 4 D(n+2) 2 (n+4) 2 -(q0-P0)a 3 16D -(P0-Pa)a 3 D(n+2) 2 (n+4) -3(q0-P0)a 2 16D -(P0-Pa)a 2 (n+3) D(n+2) 2 (n+4) -3(q0-P0)a 8D -(P0-Pa)a(n+3) D(n+2)(n+4)          and C =          C 1 C 2 C 3 C 4         
The solution matrix C has an analytic solution which is given in the Appendix A. It allows us to obtain, for any set of system parameters (h, E, ν, ρ, k, P a and n) and for any couple of control parameters (a and P 0 ), the analytical expression of the radial profile of the vertical displacement w(r).

It is important to note that there are two length scales in the model: h and l. The thickness h of the elastic strata is a parameter related to the geometry of the intrusion.

Our model is based on Equations 1 and 2, which are only valid when a/h >> 1. In the following, we will therefore only consider values of a such that a/h >5, with 5 being an arbitrarily chosen limit for the validity of the thin plate formulation, already used by e.g. [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] and [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF]. The elastic length l = 4 D k is an intrinsic length scale of the model, which represents the lateral distance, beyond the sill periphery, over which significant displacements are found. In the following, we will see that the model exhibits two regimes, the transition between which is controlled only by a/l, and not by a/h. Note however that h is involved in the value of l, via D.

One can also notice that positive displacements w are defined downward, meaning that upward displacement of the plate would be negative. Notice also that before the intrusion forms, the weight of the plate already pushes down on the elastic foundation, so that there is already a homogeneous displacement w 0 = q0 k . We will consider this equilibrium state as the initial condition when the intrusion starts forming. Consequently, in order to calculate the displacement due to the intrusion, one needs to calculate the differential displacement w i = w -w 0 = w -w(r → ∞). For practical reason, in the figures of the next sections, we plot the uplift induced by the emplacement of the intrusion, i.e. -w i .

Validation of the model

In order to test our model, we first compare it to the existing model of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF], for the same parameters (Fig. 3). The only difference is that we consider an elastic foundation of stiffness k, whereas [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] considered a perfectly rigid foundation, i.e. k=∞. This latter condition imposes no displacement outside the intrusion. In contrast, the additional elasticity introduced in our model with the elastic foundation implies that w i and its derivatives are significantly different from 0 at the intrusion periphery (Fig. 3a).

Our model is expected to converge towards that of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] when k becomes very large.

Figure 3, which shows the evolution of our model when the stiffness k of the elastic foundation is varied (the overpressure ∆P 0 = P 0 -q 0 is kept constant), demonstrates that this is indeed the case. In order to quantify how our solution converges towards that of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF], we define the relative difference χ GSP J = ( w i -w iP J )/ w iP J between our model (GS) and that of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] (PJ). In the range shown in Fig. 3b, this relative difference is observed to decrease roughly as a power law of exponent -1/4. Numerically, our analytic solution is not practically computable for any set of system and control parameters. The reason is that the Kelvin functions of decay exponentially with the value of the argument a/l and can become smaller than the numerical accuracy that Matlab can achieve. The effect of this limitation is shown in Figure 4a, in which we plot the relative uplift difference χ GSP J between our model and the one of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] as a function of a. For each given h, the analytic solution becomes unstable for high values of a, as evidenced from the peaks on the right of the curves. The main reason is that the denominator D 1 derived from the matrix inversion of Equation [START_REF] Chang | Accelerated uplift and magmatic intrusion of the Yellowstone caldera, 2004 to 2006[END_REF] (see the complete formulation in Appendix A) becomes so small that the ratio becomes unstable. The collapse of the curves displayed in Figure 4b shows that the boundary of the stability domain of our model corresponds to a critical value of the ratio a/l. This critical value of a/l is 450. In practice, therefore, we can compute the analytical solution only for a/l <450.

It is however possible to overcome this numerical limitation by using the asymptotic forms of the Kelvin functions. For large x, all Kelvin functions scale as

π 2x e -x/ √ 2 f (x),
where f is a sine or cosine function and x is the argument of the Kelvin function (see Appendix A). By cancelling out all exponentials in the inversion of the matrix of Equation ( 15), a solution is found in which no exponential is involved, and thus in which the numerical accuracy is practically never reached (see Appendix A). The relative uplift difference between this asymptotic model and that of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] is shown as the black dashed curve in Figure 4b. The asymptotic solution never becomes unstable and prolongates very nicely the analytic solution for large a/l (a/l >450). In the following, we will therefore use the analytic solution for a/l <450 and the asymptotic solution for a/l ≥450.

Figure 5 shows the relative uplift difference χ AnAsym = ( w iAn -w iAsym )/ w iAsym between the analytic and the asymptotic solutions. As was already visible in Figure 4b, the asymptotic solution is a very good approximation of the analytic solution even for relatively small values of a/l. The relative uplift difference is smaller than 1% for a/l >3.7; it is smaller than 0.1% for a/l >9.1.

Both the analytic (used when a/l <450) and the asymptotic (used when a/l ≥450) solutions have been implemented in a Matlab code, which is provided as supplementary material online.

One limitation of our model is that it is valid only when a/h > 5. In active volcanoes or volcanic systems, this is not always the case [START_REF] Amelung | Widespread uplift and "trapdoor" faulting on Galapagos volcanoes observed with radar interferometry[END_REF][START_REF] Fialko | Finite source modelling of magmatic unrest in Socorro, New Mexico, and Long Valley, California[END_REF][START_REF] Pedersen | Temporal development of the 1999 intrusive episode in the Eyjafjallajökull volcano, Iceland, derived from InSAR images[END_REF][START_REF] Chang | Accelerated uplift and magmatic intrusion of the Yellowstone caldera, 2004 to 2006[END_REF][START_REF] Sigmundsson | Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption[END_REF][START_REF] Woo | Intrusion and deformation at Campi Flegrei, southern Italy: Sills, dikes, and regional extension[END_REF]. Volcanic edifices, nevertheless, mostly consist of piles of strong lava flows intercalated between weak scoria layers. Each flow exhibits a typical thickness of a few meters. [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] and [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF] noticed that slip occurred between the layers of the overburden, indicating that the overburden behaved like a stack of thin plates of thicknesses h i for which a/h i >> 1. Therefore, one can calculate an equivalent elastic stiffness of the overburden

D e = n i=1 D i h 3 i /12,
which can replace D in Equation [START_REF] Chang | Accelerated uplift and magmatic intrusion of the Yellowstone caldera, 2004 to 2006[END_REF]. In these conditions, strata can slide over one another, and our model might become valid to analyze ground deformation due to sill emplacement in given active volcanoes even when a/h < 5 [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Scaillet | The Gangotri granite (Garhwal Himalaya): Laccolithic emplacement in an extending collisional belt[END_REF].

Behaviour of the model

In this section we illustrate the behavior of the model against changes in control and system parameters.

Effect of the control parameters a and P 0

In Figure 6a we show the evolution of the maximum uplift -w imax = w(r → ∞)w(r = 0) as a function of the overpressure ∆P = P 0 -q 0 , for different values of the intrusion radius a. For all a, -w imax increases linearly with ∆P . The slopes -w imax /∆P of the curves thus represent the effective elastic compliance of the strata. In Figure 6b, the values of this elastic compliance are plotted as a function of a for various values of E, i.e. of the bending stiffness D.

Figure 6 can be made completely non-dimensional by rescaling both axis (Fig. 7):

-w imax is divided by (q0-P0)a 4

64D

as a function of a/l. Doing so, the rescaled displacement is constant for large a/l, meaning that, in this regime, the first term in Equation 5dominates, i.e. -w imax ∼ a 4 . For small a/l, the dominating term is the second one in Equation 5, which corresponds to -w imax ∼ a 2 and yields a slope of -2 (Fig. 7). A clear cross-over is observed, for a/l 5, between the two regimes with power laws of exponent -2 and 0 for small and large a/l respectively. Physically, the regime where -w imax ∼ a 2 corresponds to a very soft elastic foundation compared to the bending plate; it means that uplift can occur over a large area outside the intrusion. In contrast, the regime where -w imax ∼ a 4 corresponds to a relatively stiff elastic foundation compared to the bending plate; this means that uplift mainly occurs upon the intrusion, with a small uplifted zone outside the intrusion. Interestingly, plotting geological values (Table 2) in Figure 7 shows that natural systems lie around the transition between these two regimes (black curve).

Effect of the heterogeneity of the pressure within the intrusion

Up to now, we have only considered situations in which the pressure is homogeneous within the intrusion (P a = P 0 ). The radial pressure distribution P (r) = P 0 -(P 0 -P a )(r/a) n can be tuned by varying the values of P 0 , P a and n (see Fig. 2d). The pressure at the center (respectively periphery) of the intrusion is P 0 (respectively P a ).

n=0 corresponds to a homogeneous pressure P a . A higher value of n yields a sharper decrease in the pressure when r approaches a (Fig. 2d). n=∞ would correspond to a homogeneous pressure P = P 0 . Outside the intrusion, the pressure is set to zero.

To test the effect of those parameters, we compare uplift profiles for distinct values of P 0 , P a and n (Figure 8). For a consistent comparison, we impose that the force F applied to the bending plate is the same for all curves. We have:

F = 2π 0 a 0 [∆P 0 -(∆P 0 -∆P a )(r/a) n ] r dr dθ = π∆P 0 a 2 1 - 2 n + 2 1 - ∆P a ∆P 0 . ( 16
)
where ∆P a = P a -q 0 . This force F will be taken equal to that yielded by a constant overpressure ∆P h =3 MPa, i.e. πa 2 ∆P h .

Figure 8a shows uplift profiles for varying ratios ∆P a /∆P 0 between 0 and 1, and for constant n = 1 (see pressure distributions in Fig. 8b). Note that cases in which ∆P a < ∆P 0 will be discussed in section 4.3. Although the total force applied to the bending plate is the same, the maximum uplift -w imax increases when ∆P a /∆P 0 decreases. This suggests that large pressures applied to the centre of the sill is more proficient to lift up the overburden of the sill.

Figure 8c shows uplift profiles for varying n between 0.5 and 20, and for constant ∆P a = 0 (see pressure distributions in Fig. 8d). The maximum uplift -w imax increases when n decreases. This also shows that larger pressures applied at the centre of the sill enhance the uplift of the sill overburden, although the total force is the same.

Discussion

Comparison with former models

Various models of sill and laccolith emplacement have been formerly designed. Most of them are based on the formulation of a clamped thin plate [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF][START_REF] Michaut | Dynamics of magmatic intrusions in the upper crust: Theory and applications to laccoliths on Earth and the Moon[END_REF], although our model shows that uplift outside the sill can be substantial. Thus, uplift outside the sill needs to be taken into consideration. [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF] and [START_REF] Michaut | Dynamics of magmatic intrusions in the upper crust: Theory and applications to laccoliths on Earth and the Moon[END_REF] solved the hydrodynamics of magma flow within the intrusion, taking into account the viscosity of, and the body forces within the magma. The pressure distribution within the intrusion is thus predicted by the calculation, and not prescribed. Nevertheless, such approaches require numerical implementation, whereas our purpose it to propose a purely analytical model. Section 4.3, however, shows how the numerical results of [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF] and [START_REF] Michaut | Dynamics of magmatic intrusions in the upper crust: Theory and applications to laccoliths on Earth and the Moon[END_REF] can be considered to constrain the pressure distribution within the intrusion.

[51] developed a 2D analytical model of laccolith emplacement based on thin bending plate laying upon an elastic foundation. Our model expands their formulation to axisymmetrical form, allowing the calculation of the volume of the intrusions. In contrast to our model, that of [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF] did not prescribe a priori the length a of the sill. This difference arises from a different boundary condition for the displacement w(r = a) at the tip of the intrusion. In our model, the only constraint on the displacement is the continuity between the two domains upon and outside the intrusion (Eq. 7). In contrast, [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF] also imposed the value of the displacement w(r = a) = 0, i.e. -w i (r = a) = q 0 /k. This boundary condition is equivalent to an implicit propagation criterion of the intrusion based on a critical displacement at the tip of the sill. Such a criterion does not take into account the mechanical properties of the system, and other propagation criteria cannot be used. In this respect, our model is more general as it can be combined to any mechanical propagation criterion separately. Such a propagation criterion establishes a relationship between the radius a of the intrusion and the other parameters like the volume of the intrusion V . Doing so, one can reduce the number of control parameters in our model from two (a, ∆P ) to only one (V ), as in e.g.[]Kerr1998. An example is developed in section 5.1.

Note that in our model, flat-lying intrusions are emplaced along a rheological boundary. This strong assumption is in agreement with geological and geophysical observations [e.g., [START_REF] Hawkes | The Sandfell Laccolith and 'Dome of Elevation[END_REF][START_REF] Johnson | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, i. Field observations, Gilbert's model, physical properties and flow of the magma[END_REF][START_REF] Thomson | Determining magma flow in sills, dykes and laccoliths and their implications for sill emplacement mechanisms[END_REF][START_REF] Hutton | Insights into magmatism in volcanic margins: bridge structures and a new mechanism of basic sill emplacement -theron mountains, antarctica[END_REF], which show that sills were mostly emplaced along weak layers such as coal, shale or unconsolidated tuff. This assumption is also in good agreement with results from laboratory experiments [e.g., [START_REF] Kavanagh | An experimental investigation of sill formation and propagation in layered elastic media[END_REF][START_REF] Menand | The mechanics and dynamics of sills in layered elastic rocks and their implications for the growth of laccoliths and other igneous complexes[END_REF][START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF][START_REF] Galland | Experimental modelling of ground deformation associated with shallow magma intrusions[END_REF], which show that layering is required to turn a vertical dyke to horizontal sill. Compressional stresses can also control the formation of horizontal magma conduits [START_REF] Hubbert | Mechanics of hydraulic fracturing[END_REF][START_REF] Sibson | Brittle-failure controls on maximum sustainable overpressure in different tectonic regimes[END_REF], and that buoyancy can play a role on the emplacement of flat-lying intrusions [e.g., [START_REF] Bradley | Intrusion of major dolerite sills[END_REF][START_REF] Francis | Magma and sediment -i. emplacement mechanism of late carboniferous tholeiite sills in northern britain[END_REF][START_REF] Taisne | Dike propagation through layered rocks[END_REF]. Recently, laboratory experiments suggested that the high pore fluid pressures contained in sedimentary rocks can also control the formation of sills [START_REF] Gressier | Control of pore fluid pressure on depth of emplacement of magmatic sills: An experimental approach[END_REF]. These latter mechanisms, however, can be important for the initiation of flat-lying conduits, whereas our model accounts for the evolution of an already initiated intrusion. Therefore, these parameters are not critical in our model.

Physical meaning of k

Sills and laccoliths are often emplaced in layered rocks, and notably along weak layers [e.g. [START_REF] Hawkes | The Sandfell Laccolith and 'Dome of Elevation[END_REF][START_REF] Duffield | Deformation of poorly consolidated sediment during shallow emplacement of a basalt sill, Coso Range, California[END_REF][START_REF] Dixon | Centrifuge modelling of laccolith intrusion[END_REF][START_REF] Román-Berdiel | Analogue models of laccolith formation[END_REF][START_REF] Kavanagh | An experimental investigation of sill formation and propagation in layered elastic media[END_REF][START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF]]. In the model presented in this paper, we propose that the elastic foundation represents a weak sedimentary layer, such as clay, tuff, or unconsolidated sediments. We thus expect this layer to deform elastically more than the overlying strata at the vicinity of an intruding sill.

An expression of k can be derived by assuming that the deformation of the elastic foundation corresponds to uniaxial strain in the case of an infinite plate, such that x = y = 0 and σ x = σ y = ν wl 1-ν wl σ z . Hooke's law thus writes [START_REF] Turcotte | Geodynamics[END_REF]:

-z = w h wl = (1 + ν wl )(1 -2ν wl ) E wl (1 -ν wl ) σ z , (17)
where h wl and E wl and ν wl are the thickness, Young's modulus and Poisson's ratio of the weak layer, respectively.

By definition, the displacement w corresponds to the response of the elastic foundation under any vertical stress σ z , such that w = σ z /k. Replacing this relation into Equation [START_REF] Daniels | The shapes of dikes: Evidence for the influence of cooling and inelastic deformation[END_REF] leads to the expression of k:

k = E wl (1 -ν wl ) (1 + ν wl )(1 -2ν wl )h wl . ( 18
)
k is a linear function of the Young's modulus E wl and is an inverse function of h wl , meaning that the effective stiffness of the elastic foundation becomes smaller if the weak layer is thicker.

From Equation 18, values of k can be estimated. The Young's modulus E wl can range given by [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF].

One can notice that for low values of k, the displacement at the tip of the sills -w i (r = a) can be substantial (Fig. 3). This has two implications: [START_REF] Abdelmalak | Fracture mode analysis and related surface deformation during dyke intrusion: Results from 2D experimental modelling[END_REF] for these inelastic processes associated with the emplacement of dykes and sills.

4.3. Tabular shape of sills and laccoliths [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF] highlight that the classical bell shapes predicted by the models of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] and [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF] are not fitting with the shapes of natural sills and laccoliths, which are tabular. In the models of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] and [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF], the pressure distribution in the sills were prescribed a priori.

Thus, the only pressure distributions considered were either homogeneous or decreasing from the centre to the tips (Fig. 2d). In contrast, [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF] implemented numerically the viscous flow and the body forces of the magma through the Navier-Stokes equations and calculated a pressure profile at the base of the deforming plate. Counterintuitively, the over-pressure at the centre of the sill ∆P 0 is negative and smaller than the over-pressure at the tip of the sill ∆P a , which is positive [Fig. 7 of 12]. In addition, the shape of their pressure profile resembles that of a function of the form of P = P 0 -(P 0 -P a )(r/a) n , where P 0 < P a and n < 1.

In Figure 9, we compare two uplift profiles calculated with our model for (1) a homogeneous pressure distribution, and (2) a pressure distribution similar to that of [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF]. For comparison purpose, the forces applied at the bottom of the deforming plate are the same (see Section 3.3.2). The homogeneous pressure gives the typical bell-shape as calculated by [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] and [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF]. In contrast, the other pressure distribution indeed produces tabular morphology, which is compatible with those observed for sills and laccoliths. Therefore, the general formulation of our model allows to take into account a physically realistic pressure distribution to calculate geologically relevant intrusion shapes. Nevertheless, in order to predict the correct pressure distribution within the sill, one would need to combine our model to the hydrodynamic equations, in a similar way as that presented by [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF], and solve the problem numerically. Note that although the force applied to the deforming plate are the same in both cases of Figure 9, the maximum uplift are substantially different.

Applications

Propagation of sills

The dynamics of sill emplacement have been studied extensively during the last decade using experimental methods [START_REF] Bunger | Near-surface hydraulic fracture[END_REF][START_REF] Kavanagh | An experimental investigation of sill formation and propagation in layered elastic media[END_REF][START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF]. These experimental studies have shown that [START_REF] Abdelmalak | Fracture mode analysis and related surface deformation during dyke intrusion: Results from 2D experimental modelling[END_REF] the fluid pressure decreases hyperbolically during emplacement at constant volumetric flow rate [START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF][START_REF] Galland | Rise and emplacement of magma during horizontal shortening of the brittle crust: Insights from experimental modeling[END_REF][START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF], and (2) the propagation velocity of the sill front increases during emplacement at constant fluid overpressure [START_REF] Kavanagh | An experimental investigation of sill formation and propagation in layered elastic media[END_REF]. These results have been confirmed

numerically [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF] and theoretically [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Pollard | On the mechanical interaction between a fluid-filled fracture and the Earth's surface[END_REF][START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF]. They indicate that the propagation of sills, or of horizontal fluid-filled fractures in general, results from a complex interplay between the pressure in the fluid, the size of the sill and the mechanical properties of the surrounding material.

Experiments of fluid injection at constant flow rate Q have shown that the fluid pressure P 0 ∝ V β , where V = Qt is the injected volume of fluid and is known at any time t. Simple scaling analysis shows that β is a function of the fracture propagation criterion [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF][START_REF] Galland | Rise and emplacement of magma during horizontal shortening of the brittle crust: Insights from experimental modeling[END_REF]. A fracture propagation criterion is in general a relationship between ∆P , the geometry of the system and the mechanical properties of the fracturing medium;

it expresses how the solid responds to the loading due to an over-pressurized fluid within a crack.

Detailed analysis of pressure data during experiments of fluid injection in solids have thus been performed to derive the dynamics of fracture formation [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF][START_REF] Galland | Rise and emplacement of magma during horizontal shortening of the brittle crust: Insights from experimental modeling[END_REF]. The goal is to fit the pressure and uplift evolution curves with the results of a mechanical model to calculate the intrinsic physical parameters of the fracturing process. [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF], for example, used the model of a clamped plate of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] to study the development of shallow flat-lying cracks in soils. In his analyses, [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF] defined a fracture propagation criterion based on the mode I stress intensity factor K I derived for a thin clamped plate attached on a rigid foundation, which writes:

K I = ∆P a 2 3 32h 3 1/2 . ( 19
)
When K I exceeds a critical value K IC (also called the fracture toughness, which is a material property), the fracture tip becomes unstable and propagates.

Here we propose to develop a similar approach to those of [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF] but using our analytical model. In our model, the propagation of an intrusion occurs when the elastic foundation fails, i.e. when it reaches a given mechanical threshold, e.g. a critical stretching. We made the assumption that the failure of the elastic foundation can be expressed by a critical value of the stress intensity factor. Equation 19, however, is valid only for the clamped model of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] and [START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF], but not for our model. Nevertheless, the derivation of an analytic formula of the stress intensity factor for our model proved to be very challenging, and we propose to calculate it numerically (see Appendix B). Such a procedure produces the relationship ∆P ∝ 1/a 2 (Fig. 10), which is similar to that of Equation 19.

In order to compare our model with the analyses of [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF] and [START_REF] Galland | Rise and emplacement of magma during horizontal shortening of the brittle crust: Insights from experimental modeling[END_REF], we need to derive the volume V of the intrusion, which writes in cylindrical coordinate system:

V = - 2π 0 a 0 rw i (r) dr dθ = -2π a 0 rw i (r)dr. ( 20
)
Combining Equations ( 5) and ( 20) and integrating, the expression of the volume of the sill at any time t writes:

V (t) = -2π (q 0 -P 0 )a 6 384D + C 1 a 4 16 + C 2 a 2 2 - q 0 a 2 2k . ( 21
)
For simplicity, we consider here the case of homogeneous pressure within the sill (P a = P 0 ). According to Equations ( 32) and [START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF] given in the Appendix A, C 1 and C 2 are complicated functions of a and ∆P . Combining the relationship between ∆P and a (Fig. 10) with Equation ( 21) implies that the volume V of the sill can be calculated numerically as a function of a or ∆P only, such that a and ∆P can be plotted as functions of V (Fig. 11). Considering a constant injection flow rate Q, V = Qt can directly be converted into time t.

Figures 11a, b, andc show the evolution of ∆P , a and -w imax with V , respectively.

Qualitatively, the curves of Figures 11a,b, and c match well the experimental results of shallow horizontal intrusions of [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 1. Methods and observations[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF] and [START_REF] Galland | Experimental modelling of ground deformation associated with shallow magma intrusions[END_REF]. In log-log representations (Fig. 11), the data produce linear trends, thus showing that ∆P , a and -w imax are all proportional to V α . The slopes in Figure 11 are -1/2, 1/4 and 1/2, respectively, showing that

∆P ∝ 1/ √ V , a ∝ 4 √ V and -w imax ∝ √ V .
The very same relationships have been obtained experimentally by [START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF], suggesting that our model is relevant for describing the emplacement of shallow flat-lying intrusions.

Notice that the domain of validity of the thin plate approximation is for a/h > 5.

Some authors, however, used this formulation for small values of a/h [e.g., [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Scaillet | The Gangotri granite (Garhwal Himalaya): Laccolithic emplacement in an extending collisional belt[END_REF]. We recall that, for the sake of rigor, we chose to use our model only for a/h >5, and to represent only such data points in all figures. If one would be interested in the case where a/h < 5, e.g. for small sills, further tests and comparisons with elastic half-space models [e.g., [START_REF] Sun | Theoretical size of hydraulically induced horizontal fractures and corresponding surface uplift in an idealized medium[END_REF][START_REF] Fialko | Deformation due to a pressurized horizontal circular crack in an elastic half space, with applications to volcano geodesy[END_REF] would be required.

[51] also derived an analytical solution of the surface uplift as a function of V . In their approach, nevertheless, the mathematical problem was over-determined, such that a crack propagation criterion was implicitly introduced in their boundary conditions: the fixed displacement imposed at the tip of the sill yielded a relationship between a and P [Eq. ( 23) of 51], playing a role analogous to our relationship between a and ∆P in Figure 10. We emphasize that our approach is different and more general: we solved the mathematical problem in the most generic manner for any values of the system and control parameters, and we subsequently used an additional physical definition of a crack propagation criterion (Fig. 10) that links the pressure and the radius of the intrusion.

With this criterion, the radius a of the intrusion is not an independent control parameter anymore (like in the analysis of [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF]), and the surface uplift can be calculated as a function of e.g. the volume V only. Because the application of a crack propagation criterion is separated from the main analysis, any crack propagation criterion can be used with our model, in contrast to the model of [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF].

Critical size of saucer-shaped sills

Saucer-shaped sills result from the mechanical interactions between a growing horizontal sill and the deforming overburden [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF][START_REF] Polteau | Saucer-shaped intrusions: occurrences, emplacement and implications[END_REF][START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF]. The magma overpressure within the sill lifts up the overburden to form a gentle dome, at the rim of which stresses interact with the leading edge of the sill [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF]. When these stresses reach a critical value, the sill tip is deviated from the horizontal, and it propagates upward to form inclined sheets.

Such a phenomenon is fundamental in sedimentary basins as it controls the formation of The analysis developed in this paper can be used to predict under which condition the sill-to-inclined sheet transition occurs [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF]. The bending of the overriding plate generates a strain distribution along the plate due to outer arc and inner arc deformation (Fig. 12). In the formulation of a thin bending plate, the vertical normal strain z = ∂w/∂z within the plate and upon the sill is zero, because w is a function of r only.

In contrast, r = ∂u/∂r is not zero because the radial displacement u is a function of r. Notably, the tip of the sill (r = a) is submitted to tensile strains ra (Fig. 12a). We assume that the overriding plate fails at the tip of the intrusion when ra reaches a critical value c (radial tensile strain at failure), leading to the formation of an inclined sheet due to mode I fracture (Fig. 12b). This reasoning is in good agreement with experiments [Fig. 12b; 69] and geological observations [START_REF] Muirhead | Interconnected sills and inclined sheet intrusions control shallow magma transport in the ferrar large igneous province, antarctica[END_REF], which show that upward open fractures, i.e. inclined sheets, form at the tips of sills. Nevertheless, we cannot rule out that shear deformation also affects the overriding plate.

The above analysis has been applied by [START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF] using the formulation of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF], i.e. that of a clamped plate. They found the expression of the critical size a c at which sills turn to inclined sheets a c = 2h

E c 3(1-ν 2 )∆P .
Here we apply a similar analysis using our more general solution. To this end, one needs to calculate the strain at the tip of the sill, i.e.

ra . The expression of the radial strain r at any point of the bending plate is given by:

r = ∂u ∂r = -z ∂ 2 w ∂r 2 , ( 22
)
where u is the radial displacement and z is the distance from the neutral surface. Combining Equations ( 5) and ( 22) yields an expression of r :

r = -z 3(q 0 -P 0 )r 2 16D + C 1 (a, P 0 ) 2 + (n + 3)(P 0 -P a )r n+2 Da n (n + 2) 2 (n + 4) , r < a. ( 23
)
The points at the contact between the sill and the plate are at z = h/2. The transition from a horizontal sill to an inclined sheet is assumed to occur when the bending plate above the tip of the sill fails, i.e. when ra = c at z = h/2. The strain at the tip of the sill can be calculated from Equation [START_REF] Fialko | Evidence for on-going inflation of the Socorro magma body, New Mexico, from interferometric synthetic aperture radar imaging[END_REF],

ra = - h 2 3(q 0 -P 0 )a 2 16D + C 1 (a, P 0 ) 2 + (n + 3)(P 0 -P a )a 2 D(n + 2) 2 (n + 4) . (24)
Here, we will consider the case of homogeneous pressure within the sill, i.e. P 0 = P a .

Equation ( 24) thus becomes:

ra = - h 2 3(q 0 -P 0 )a 2 16D + C 1 (a, P 0 ) 2 , (25)
with the constant C 1 given by Equation [START_REF] Galland | Rise and emplacement of magma during horizontal shortening of the brittle crust: Insights from experimental modeling[END_REF]. We can calculate the strain r at the tip of the sill for any values of P 0 and a. Here we consider that E = 10 10 Pa and ν = 0.35.

Figure 13 illustrates how ra evolves versus a, for given overpressures ∆P . The transition from horizontal sill to inclined sheet occurs when ra reaches the critical value c . [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF] have suggested that c may lie in the range 2 × 10 -3 to 2 × 10 -2 . Thus, the size a c of the sill at which the transition occurs is the solution of:

c = - h 2 3(q 0 -P 0 )a 2 c 16D + C 1 (a c , P 0 ) 2 . (26)
In the case of a clamped plate, C 1 is a simple function of a and P 0 (C 1 = P 0 a 2 /16D), such that a c can easily be found analytically. Here, Equation ( 26) is not directly solvable because C 1 is a complicated function of both a and P (see Eq. 32). Nevertheless, Equation ( 26) can be solved numerically, as shown on Figure 13. Figure 14 shows a map of the calculated values of a c versus ∆P and h. For clarity, we plot profiles of this map (white lines), which are shown on Figures 15 and16.

The values of a c calculated with our model, by using typical geological parameters, are typically tens of kilometers for sill depths of a few kilometers. This result is in good agreement with sills observed in the field [e.g., 54, 52, 70, 33, 30] and on seismic profiles [e.g., [START_REF] Thomson | Volcanic features of the North Rockall Trough: application of visualisation techniques on 3D seismic reflection data[END_REF][START_REF] Thomson | Geometry and growth of sill complexes: insights using 3D seismic from the North Rockall Trough[END_REF][START_REF] Hansen | Saucer-shaped sill with lobate morphology revealed by 3D seismic data: implications for resolving a shallow-level sill emplacement mechanism[END_REF][START_REF] Thomson | Determining magma flow in sills, dykes and laccoliths and their implications for sill emplacement mechanisms[END_REF][START_REF] Hansen | Miocene igneous activity in the Northern Subbasin, offshore Senegal, NW Africa[END_REF]. This suggests that our model is a reasonable approximation of reality.

Figure 15 shows the difference between the clamped solution of [START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF] and the results of our model for three different values of k, with a constant h = 2000 m. The clamped solution of [START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF] follows a hyperbolic function of ∆P of the form 1/ √ ∆P . When k = 10 10

Pa m -1 , the elastic foundation is very rigid, almost all the deformation is accommodated by the bending plate and consequently, our solution provides results very similar to the clamped solution. In contrast, when k decreases the elastic foundation is softer and it accommodates part of the deformation. As a consequence, the bending of the plate at the tip of the sill is less pronounced. Therefore, although a c follows a hyperbolic trend similar to that of the clamped solution, the calculated values are higher (Fig. 15).

This has the important physical meaning that an elastic foundation favors the horizontal spreading of sills. In addition, for soft elastic foundations, the linear behavior in the log-log plot of Figure 15b is lost, showing that a c cannot be described by a hyperbolic function of ∆P anymore.

Figure 16 shows the relationship between a c and h for different values of k. According to [START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF], a c is a linear function of h in the case of a clamped plate. This is consistent with geological and geophysical observations [START_REF] Polteau | Saucer-shaped intrusions: occurrences, emplacement and implications[END_REF]. Similarly to Figure 15 and for the same reasons, the values of a c calculated from our solution deviates from the clamped solution. In addition, we can notice for the case k = 10 5 Pa m -1 that the relationship between a c and h in a log-log plot is not perfectly linear anymore, but the curve slightly bends downward. This latter behavior has been observed in numerical simulations [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF] and in experiments [START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF]. This result once again suggests that the solution presented in this paper is more realistic than the simple clamped solution.

Comparison with the model of [58]

In active volcanoes, the shape and dynamics of intrusions can be inferred from geophysical [START_REF] Brandsdóttir | Thin low-velocity zone within the Krafla caldera, NE-Iceland attributed to a small magma chamber[END_REF][START_REF] Brandsdóttir | The seismic structure of Iceland[END_REF] and geodetic data [START_REF] Amelung | Widespread uplift and "trapdoor" faulting on Galapagos volcanoes observed with radar interferometry[END_REF][START_REF] Fialko | Deformation due to a pressurized horizontal circular crack in an elastic half space, with applications to volcano geodesy[END_REF][START_REF] Battaglia | Evidence for fluid migration as the source of deformation at Campi Flegrei caldera (Italy)[END_REF][START_REF] Pedersen | Temporal development of the 1999 intrusive episode in the Eyjafjallajökull volcano, Iceland, derived from InSAR images[END_REF][START_REF] Chang | Accelerated uplift and magmatic intrusion of the Yellowstone caldera, 2004 to 2006[END_REF][START_REF] Vasco | Crustal deformation and source models of the Yellowstone field from geodetic data[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Sigmundsson | Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption[END_REF][START_REF] Woo | Intrusion and deformation at Campi Flegrei, southern Italy: Sills, dikes, and regional extension[END_REF]. In order to estimate their geometrical and dynamic characteristics, a standard procedure consists of (1) making an assumption on the shape of the intrusion responsible for the monitored ground deformation, and (2) fitting ground deformation data monitored from geodetic survey (GPS, InSAR) with an analytical solution or a numerical model relevant to the assumed intrusion shape. The parameters of the model that provide the best fit are considered to characterize the shape of the intrusion.

Analytical solutions are convenient to use as they often do not require heavy computational skills. A classical analytical model used for analyzing ground deformation in volcanic systems is that of [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surface around them[END_REF]. The "Mogi point source" model has been used extensively to invert ground deformation data measured on active volcanoes when the shape of the magma source was not known a priori [e.g. [START_REF] Masterlark | Magma intrusion and deformation predictions: Sensitivities to the Mogi assumptions[END_REF][START_REF] Bathke | Inflation and deflation at the steep-sided Llaima stratovolcano (Chile) detected by using InSAR[END_REF][START_REF] Ji | Transient signal detection using GPS measurements: Transient inflation at Akutan volcano, Alaska, during early 2008[END_REF][START_REF] Ofeigsson | Deep magma storage at Hekla volcano, Iceland, revealed by InSAR time series analysis[END_REF]. This model accounts for the vertical (w M ) and horizontal (u M ) displacements due to an over-pressured spherical magma body of radius a m at a depth h m >> a m :

u M = 3a 3 m ∆P m 4µ r (h 2 m + r 2 ) 3/2 , (27)
w M = 3a 3 m ∆P m 4µ h m (h 2 m + r 2 ) 3/2 , ( 28
)
where r is the radial distance on the surface, and µ = Eν/(1 + ν)(1 -2ν).

In general, the horizontal displacements produced by models based on thin bending plates are not considered [e.g., [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF][START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF]. One reason is that u = 0 at the neutral line at the centre of the plate, i.e. at z = 0, because the radial strain r = 0 (Fig. 12a).

Nevertheless, above and below the neutral line, r = 0 as outer-arc stretching and innerarc shortening occurs (Fig. 12a). Combining Equation 22with Equations 5 and 6 lead to the expressions of u at the surface of the model, i.e. at z = -h/2:

u(r < a) = h 2 (q 0 -P 0 )r 3 16D + C 1 2 r + (P 0 -P a )r n+3 Da n (n + 2) 2 (n + 4) , (29)
u(r > a) = h 2 C 3 √ 2l kei 1 r l -ker 1 r l + C 4 √ 2l kei 1 r l + ker 1 r l . ( 30
)
We compare the vertical (-w i ) and horizontal (u) displacements computed from our analytical model with standard values (h = 2000 m, a = 10000 m) with the Mogi point source model [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surface around them[END_REF] by fitting this latter to our model (Fig. 17). When the fitting is applied to the domain 0 < r < 2a, the shapes of the curves exhibit some differences (Fig. 17a), but the overall shapes are similar. Interestingly, the result of the fit gives an estimate of the depth of the Mogi point pressure source h m = 7550 m, i.e. the Mogi point source model greatly over-estimates the depth of the magma source (by more than a factor of 3.5). One can notice that the ratio -w imax /u max is larger in our model than in the model of [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surface around them[END_REF], i.e. the horizontal displacements u calculated from our model are proportionally smaller than u M .

Notice that the failure of the Mogi model is mainly due to the geometry of the intrusion, rather than to the details of the sill model. As an illustration, in Figure 17b, we fitted the Mogi point source model to the clamped model of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF]. The result also

shows that the Mogi model over-estimates the depth of the magma source (h m = 6681 m versus h = 2000 m), but to a lower extent than with our model.

Because ground deformation data are often produced from InSAR measurements [e.g.

3, 66, 53, 77], only the vertical displacement w can be calculated. We thus fitted the vertical displacements -w i calculated from our model with the vertical component w M of the Mogi model (Fig. 17c) on the domain 0 < r < 2a. The fit appears much better than in Figure 17a, but the calculated depth of the pressure source is h m = 7439 m, i.e.

the same depth as that calculated from fitting both -w i and u (Fig. 17a). Applying the Mogi model to analyze ground deformation induced by the emplacement of a sill, and monitored with InSAR also leads to large over-estimates of the depth of the pressure source (by a factor of 3.5).

All these tests show that the systematic use of the Mogi point source model [e.g., 78] would lead to substantial misunderstanding of the structure and dynamics of volcanic systems if the magma reservoir is shallow and flat-lying, such as in many volcanoes [e.g., [START_REF] Amelung | Widespread uplift and "trapdoor" faulting on Galapagos volcanoes observed with radar interferometry[END_REF][START_REF] Fialko | Finite source modelling of magmatic unrest in Socorro, New Mexico, and Long Valley, California[END_REF][START_REF] Sigmundsson | Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption[END_REF][START_REF] Woo | Intrusion and deformation at Campi Flegrei, southern Italy: Sills, dikes, and regional extension[END_REF]. This large error can have important implications for, e.g. interpreting petrological data for inferring the depth of magma reservoirs. Constraining properly the depth of magma reservoir would also avoid surprises such as at Krafla Volcano, Iceland, where drilling unexpectedly ended up directly in the magma chamber at about 2 km depth [e.g., 21].

Our analysis suggests that inversion of ground deformation data based on analytical models may not consider the Mogi point source model only, but should systematically complement this latter with another simple analytical model that takes into account a flat-lying shape, like our model. In this respect, a systematic test of the suitability of our analytical model for ground deformation inversion, and especially a comparison with the solutions of [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF] and [START_REF] Fialko | Deformation due to a pressurized horizontal circular crack in an elastic half space, with applications to volcano geodesy[END_REF] would be very interesting to perform.

Conclusions

In this paper, we developed a new axi-symmetrical model of surface uplift due to flat-lying magma intrusions. The model is based on the formulation of a thin bending plate lying on a deformable elastic foundation. In contrast to former models of sills and laccoliths considering the thin bending plate formulation [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF][START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF][START_REF] Scaillet | The Gangotri granite (Garhwal Himalaya): Laccolithic emplacement in an extending collisional belt[END_REF], our model allows uplift both upon and outside the intrusion. The model proposed here is an improved version of that proposed by [START_REF] Kerr | Toward more realistic formulations for the analysis of laccoliths[END_REF], as we extend it to axi-symmetrical and we set generic boundary conditions. The pressure distribution within the intrusion is P = P 0 -(P 0 -P a )(r/a) n , where P 0 and P a are the pressure at the center (r=0) and periphery (r=a) of the intrusion, respectively.

We propose a fully analytical solution of the model. Numerically, however, the solution become unstable for a/l > 450, with a the radius of the sill and l the elastic length emerging from the model. In these cases, we propose an asymptotic solution which prolongates the analytic solution for large values of a/l. Both the analytic and asymptotic solutions are provided in a Matlab file as supplementary material.

The uplift calculated from our model depends on both the elastic properties of the bending plate and of the elastic foundation. We show that when the elastic foundation is very stiff, i.e. when k → ∞, our model converges towards that of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] for a clamped plate, as expected.

The model exhibits two regimes depending on the ratio a/l. When a/l < 5, the maximum uplift at the center of the intrusion -w imax evolves as a 2 (Fig. 7). In this regime, the elastic foundation is soft compared to the bending plate, the calculated uplift spreads over a large domain compared to the intrusion area. In contrast, when a/l > 5,

the maximum uplift at the center of the intrusion evolves as a 4 . In this latter regime, the bending plate is soft compared to the elastic foundation, the uplift extends mostly upon the sill. Interestingly, geological systems are expected to be found on both sides of the transition between the two regimes.

By introducing a relevant fracture propagation criterion, we show that the model can be used to describe sill propagation. For realistic values of the model parameters, our solution reproduces well the temporal evolution of shallow horizontal intrusions simulated in experiments [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF][START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF][START_REF] Galland | Experimental modelling of ground deformation associated with shallow magma intrusions[END_REF].

The model can also be used to predict the critical size a c at the transition from inner sill to inclined sheet in saucer-shaped sills as a function of h and ∆P . a c approaches linear functions of h and 1/ √ ∆P , as observed in the field [START_REF] Polteau | Saucer-shaped intrusions: occurrences, emplacement and implications[END_REF], in experiments [START_REF] Galland | Experimental modelling of shallow magma emplacement: Application to saucer-shaped intrusions[END_REF] and in numerical simulations [START_REF] Malthe-Sørenssen | Formation of saucer-shaped sills[END_REF]. We also show that soft elastic foundations favours the horizontal spreading of sills before they form inclined sheets.

We compared the vertical and horizontal displacements calculated from our model with the Mogi point source model [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surface around them[END_REF], and we show that this latter can strongly overestimate the depth of a flat-lying intrusion by more than a factor of 5. We thus propose that our model can be useful to analyze ground deformation resulting from sill intrusions in active volcanoes and offers a simple alternative to numerical models.

Appendix A

The elements C i of the solution matrix of Eq. 15 are such that:

C 1 D i 2a = 4b 6(h i m i -g i n i ) + 3(g i q i -h i p i )a + (n i p i -m i q i )a 2 + (31) s (h i m i -g i n i )(24 + 26n + 9n 2 + n 3 ) + (g i q i -h i p i )(12 + 7n + n 2 )a + (n i p i -m i q i )(4 + n)a 2 C 2 - q 0 k -5b -s 1 - (n + 4)(n + 3) 2 a 4 D i a = (32)
c i h i -f i g i + a (f i m i -c i n i ) - a 2 2 (h i m i -g i n i ) (24b + s(2 + n)(3 + n)(4 + n)) + a 2 c i q i -f i p i - a 2 2 (m i q i -n i p i ) (8b + s(4 + n)(2 + n)) C 3 D i a = 8b 3(n i a -h i ) -q i a 2 + s (24 + 26n + 9n 2 + n 3 )(n i a -h i ) -(8 + 6n + n 2 )q i a 2 (33) C 4 D i a = (g i -m i a) [24b + s(2 + n)(3 + n)(4 + n)] + p i [8b + s(4 + n)(2 + n)] a 2 (34) with b = q0-P0 64D , s = P0-Pa D(n+2) 2 (n+4) 2 , D i = h i p i -g i q i + (m i q i -n i p i )a
and the following definition of the coefficients for the analytic (i = 1) or the asymptotic (i = 2) solutions:

c 1 = -kei 0 ( a l ) (35)
f 1 = -ker 0 ( a l ) (36)
g 1 = - kei 1 ( a l ) -ker 1 ( a l ) √ 2l (37)
h 1 = - kei 1 ( a l ) + ker 1 ( a l ) √ 2l (38)
m 1 = - ker 0 ( a l ) -ker 2 ( a l ) 2l 2 (39)
n 1 = - kei 2 ( a l ) -kei 0 ( a l ) 2l 2 (40)
p 1 = - 3ker 1 ( a l ) -ker 3 ( a l ) + 3kei 1 ( a l ) -kei 3 ( a l ) 4 √ 2l 3 (41)
q 1 = - 3ker 1 ( a l ) -ker 3 ( a l ) -3kei 1 ( a l ) + kei 3 ( a l ) 4 √ 2l 3 (42)
and

c 2 = πl 2a sinσ (43)
f 2 = -πl 2a cosσ (44)
g 2 = √ π 4l 2 ( a l ) 3 2 2acosσ -( √ 2l + 2a)sinσ (45)
h 2 = √ π 4l 2 ( a l ) 3 2 2asinσ + ( √ 2l + 2a)cosσ (46)
m 2 = √ π( a l ) 3 2 8a 4 l(3 √ 2l + 4a)sinσ -4a(l + √ 2a)cosσ (47)
n 2 = - √ π( a l ) 3 2 8a 4 l(3 √ 2l + 4a)cosσ + 4a(l + √ 2a)sinσ (48)
p 2 = √ πa l 3 16a 4 2a(4a 2 + 6 √ 2la + 9l 2 )cosσ + (8a 3 -18l 2 a -15 √ 2l 3 )sinσ (49)
q 2 = √ πa l 3 16a 4 2a(4a 2 + 6 √ 2la + 9l 2 )sinσ -(8a 3 -18l 2 a -15 √ 2l 3 )cosσ (50)
with σ = a l √ 2 + π 8 .
with α a constant independent of ∆P . Combining this simple relationship with Eq. 51 leads to:

U = - V 0 V dV α = V 2 2α = 1 2 ∆P V. (53)
Combining Eqs. 53 and 21 leads to:

U = π∆P q 0 k -C 2 a 2 2 - C 1 16 a 4 - q 0 -P 0 384D a 6 . (54)
The stress intensity factor is given by the following expression:

K I = √ GE , (55)
where G = -dU/dA, E = E/(1 -ν 2 ), and A is the area of the intrusion. G is a function of a and P , and can be rewritten as:

G(a, P ) = - dU da 1 2πa . (56)
U being a very complex function of a through C 1 and C 2 , G cannot be derived analytically. However, for a given pressure, G can be calculated numerically as:

G(a, P ) = - U (a + da) -U (a) da 1 2πa . (57)
Combining Equations 55, 57 and 54, and giving K I a critical value K IC , i.e. the fracture toughness, leads to a relationship between P and a, which can be calculated numerically (Fig. 10). Combining this numerical relationship with Equation 21 allows one to calculate the evolution of a, P and -w imax as a function of V ( Fig 11).

In order to test the relevance of this approach, we compare it with the clamped plate problem with homogeneous pressure, in which:

C 1 4 = 2∆P a 2 64D , (58)
q0 k -C 2 = ∆P a 4 64D . (59)
Combining Eqs. 58 and 59 with Eq. 54 leads to an expression of the input energy:

U = π∆P 2 a 6 384D . ( 60
)
Using Eq. 60 in the definition of G (Eq. 56) leads to an expression of G:

G = 3∆P 2 a 4 32E h 3 . (61)
Finally, the stress intensity factor can be derived from Eq uation 55:

K I = ∆P a 2 3 32h 3 . ( 62
)
This expression is identical to that derived by [START_REF] Murdoch | Hydraulic fracturing of soil during laboratory experiments. Part 3. Theoretical analysis[END_REF] and [START_REF] Murdoch | Mechanical analysis of idealized shallow hydraulic fracture[END_REF] (see also Eq. [START_REF] Duffield | Deformation of poorly consolidated sediment during shallow emplacement of a basalt sill, Coso Range, California[END_REF], thus validating our calculation. a. Unloaded system. A plate of thickness h is located above an elastic foundation of stiffness k. Positive vertical axis is downward. The lithostatic stress of the plate, q 0 , is not yet applied on the elastic foundation. b. Gravity loading only. Application of the weight of the plate yields a homogeneous vertical displacement w 0 = q 0 /k. c. System after a sill intrusion. A (possibly heterogeneous) pressure distribution P (r) is applied within an axisymmetric sill of radius a. The displacement w is calculated with respect to state a (dashed line). In contrast, the uplift w i due to the sill is calculated with respect to state b, i.e. w i = w -w 0 . d. Plot of the radial distribution of the pressure P (r) in the intrusion, normalized by P 0 , as a function of r/a, for various values of n (n increases from 0 to 10 by steps of 1). The pressure at the intrusion periphery was chosen, for this illustration, such that P a = P 0 /4. Figure 3: a. Plots of the radial profile of the surface topography -w i for a=5000 m, ∆P 0 = P 0 -q 0 =1 MPa, h=1000 m, E=10 GPa, ν=0.35, ρ=2500 kg m -3 and P a =P 0 .

Dashed black line: model of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF]. Solid grey lines: this model for the same parameters, and with k=10 5 , 10 6 , 10 8 , 10 10 Pa m -1 . Horizontal arrow with r = a locates the sill. b. Log-log plot of the relative uplift difference between this model and the model of [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] χ GSP J = ( w i -w iP J )/ w iP J , where GS refers to this model (Galland and Scheibert) and PJ refers to [START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF]. The dashed line has a slope of -1/4. and P a =P 0 , but any other set of parameters would have produced the same curve (cf. collapse in Fig. 4b). The slope of the line being -3 shows that the asymptotic solution converges quickly towards the analytical solution. P 0 and P a such that ∆Pa ∆P0 varies from 0 to 1 with a step of 0.1. Bottom line: varying n = [0.5, 0.7, 1, 2, 50, 10, 20]. P 0 is calculated so that the force on the plate is the same constant as in the top line. Note that we consider P a = q 0 , such that ∆P a = 0 Pa. Solid grey line: for ∆P = P 0 -(P 0 -P a )(r/a) n -q 0 , where ∆P 0 = -0.5 MPa, ∆P a = 5

MPa, n=5, a=5000 m, h=1000 m, ρ=2500 kg m -3 , E=10 GPa, ν=0.35, and k=5 MPa m -1 . Note that ∆P 0 <0, as calculated by [START_REF] Bunger | Modeling the growth of laccoliths and large mafic sills: Role of magma body forces[END_REF]. Black dashed line: for homogeneous pressure ∆P h . Note that the value of ∆P h = P 0 -2(P 0 -P a )/(n + 2) -q 0 , where P 0 and P a are the values from the heterogeneous case of a, ensuring that the total force applied to the bending plate is the same in both cases. b. Plots of over-pressure distributions for both cases presented in a). This study (see Section 4.2) ∆P 10 6 -10 7 Pa e.g. [START_REF] Rubin | Propagation of magma-filled cracks[END_REF]

from 10 5

 5 Pa for unconsolidated tuff or sediments, or a paleosoil [e.g. 63, 4, 2] to 10 9 Pa for consolidated tuff of more competent rocks. Considering realistic thickness h wl between 1 to 10 meters for the weak layer, and considering ν = 0.35, this leads to a geological range for k between 10 4 and 10 9 Pa m -1 . These values have been used in Figure 7 when defining the geologically relevant model parameters. This range is consistent with the value of k = 2-5×10 7 Pa m -1

( 1 )

 1 numerous magmatic sills in the Karoo Basin [e.g.,[START_REF] Chevallier | Morpho-tectonics and mechanism of emplacement of the dolerite rings and sills of the western Karoo, South Africa[END_REF][START_REF] Polteau | Saucer-shaped intrusions: occurrences, emplacement and implications[END_REF][START_REF] Goulty | Implications of simple flexure theory for the formation of saucershaped sills[END_REF][START_REF] Galerne | 3d relationships between sills and their feeders: evidence from the Golden Valley Sill Complex (Karoo Basin) and experimental modelling[END_REF], offshore UK and Norway [e.g.,[START_REF] Hansen | The three-dimensional geometry and growth of forced folds above saucer-shaped igneous sills[END_REF][START_REF] Thomson | Determining magma flow in sills, dykes and laccoliths and their implications for sill emplacement mechanisms[END_REF], and (2) many sand intrusions [e.g.,[START_REF] Duranti | Fluidization and injection in the deep-water sandstones of the Eocene Alba Formation (UK North Sea)[END_REF][START_REF] Huuse | Eocene sandstone intrusions in the Tampen Spur area (Norwegian North Sea Quad 34) imaged by 3D seismic data[END_REF][START_REF] Szarawarska | Three-dimensional seismic characterisation of large-scale sandstone intrusions in the lower Palaeogene of the North Sea: completely injected vs. in situ remobilised sandbodies[END_REF].

Figure 1 :

 1 Figure 1: Geological examples of uplift due to sill emplacement. a. Satellite image of Eyjafjallajökull Volcano, Iceland, showing InSAR and seismic data during the Spring 2010 eruption [after 77]. The data were monitored between 25 September 2009 and 20 March 2010. They were related to pre-eruptive period, which corresponded to the emplacement of a ≈10-km large sill at a few kilometers depth [transparent white surface; Figure 3e of 77]. Black arrows show the satellite flight path (downward arrow) and look directions (leftward arrow). Black dots show earthquake epicenters during this period. The red stars locate the eruption localities. The yellow triangles locate GPS stations that monitored continuously flank deformation. The colored fringes represent ground displacement calculated from TerraSAR-X interferograms from descending satellite orbits. Each fringe corresponds to line-of-sight, i.e. distance from the satellite, change of 15.5 mm. The total displacement can thus be several tens of centimeters. Background is shaded topography. The rounded patterned line on the right of the image locates the caldera of Katla Volcano. Note that the uplifted area is at least twice wider than the underlying sill. b. Seismic profile illustrating the relationships between a saucer-shaped sill and the structure in its overburden, Rockall Basin, offshore Scotland [modified after 39]. Vertical scale is the time for seismic wave travel (in seconds). The profile shows that sill overburden is bent, forming a dome structure (uplift). This profile shows that the dome is about 1.3 times wider than the sill.

Figure 2 :

 2 Figure 2: Schematic diagram of the model developed in this paper [modified after 51].

Figure 4

 4 Figure4: a. Plot of the relative uplift difference χ GSP J between this model and that of[START_REF] Pollard | Mechanics of growth of some laccolithic intrusions in the Henry Mountains, Utah, II. Bending and failure of overburden layers and sill formation[END_REF] as a function of intrusion radius a, for various strata thickness h (h is varied from 25m to 5000m). The lower bound of h is geologically unrealistic in purpose. See definition of χ GSP J in caption of Fig.3. ∆P 0 =2 MPa, E=10 11 Pa, ν=0.35, ρ=2500 kg m -3 , k=1 MPa m -1 and P a =P 0 . Note that for each curve, the data becomes numerically unstable above a critical value of a. b. The same data is plotted in solid grey line as a function of a/l. All curves collapse on a single curve. The dashed line shows the asymptotic form of our model. Notice that collapse of the dashed lines make them appear as a solid black line. Note that only points corresponding to a/h >5 have been considered. A lower thin dashed line illustrates the slope of -1. Note that the analytical values of χ GSP J become numerically unstable for the same critical value of a/l 450

Figure 5 :

 5 Figure 5: Plot of the relative uplift difference χ AnAsym between the analytic solution and the asymptotic solution as a function of a/l. Curve obtained for ∆P 0 = 2 MPa, h ranging between 25 m and 5000 m, E=10 11 Pa, ν=0.35, ρ=2500 kg m -3 , k=1 MPa m -1

Figure 6 :

 6 Figure 6: a. Plot of the maximum uplift -w imax as a function of ∆P for various a (a varies from 5000m to 10000m by steps of 500 m). h=1000 m, k=1 MPa m -1 , E=100GPa, ν=0.35, ρ=2500 kg m -3 and P a =P 0 . b. Plot of the effective elastic compliance of the strata -w imax /∆P as a function of a, for various E (E=10 i Pa with i increasing from 9 to 12 by steps of 0.5). h=1000 m, k=1 MPa m -1 , ν=0.35, ρ=2500 kg m -3 and P a =P 0 . Note that the pressure is constant, so that ∆P =∆P 0 .

2 Figure 7 : 7 Figure 8 :

 2778 Figure 7: Plot of -64Dw imax /∆P a 4 as a function of a/l, for various E (E=10 i Pa with i an integer from 7 to 16), h varying between 10 and 10000 m, a varying between 100 and 10 8 m, k varying between 10 3 and 10 11 Pa m -1 , ν=0.35, ρ=2500 kg m -3 , ∆P varying between 10 6 and 10 7 Pa, and P a =P 0 . We chose in purpose range values that are not relevant for geological systems in order to test the generic behavior of the model (grey curve). All the data scale on one master curve, which exhibits two regimes. Geologically relevant values of the model parameters are indicated with the black curve, and fall right at the transition between the two regimes. Geological values are ∆P between 10 6 and 10 7 Pa, h between 1000 and 5000 m, a between 1000 and 10 5.5 m, E between 10 9 and 10 11 Pa, and k between 10 4 and 10 9 Pa m -1 (see Section 4.2). These values gives values of a/l ranging between 0.0539 m and 3203 m. Only points corresponding to a/h >5 have been considered. Note that the pressure is constant, so that ∆P =∆P 0 .

Figure 9 :

 9 Figure 9: a. Plot of the uplift profile -w i for two pressure distributions within the sill.

Figure 10 : 2 Figure 11 :Figure 12 :Figure 13 :

 102111213 Figure 10: a. Plots of the overpressure ∆P as a function of a for h = 500 m (solid line), h = 2000 m (dashed line), and h = 5000 m (dashed-dotted line). Here we chose K IC = 10 1 0 Pa m 1/2 , E = 10 10 Pa and ν = 0.35, ρ = 2000 kg m -3 , k=0.5 MPa m -1 . b.Same data in a log-log plot: the slope is -2, showing that ∆P ∝ 1/a 2 . This is similar to the relation between ∆P and a in Equation 19 for the clamped model.

Figure 14 :Figure 15 : 5

 14155 Figure 14: Diagram representing the calculated critical size a c of a sill as a function of ∆P and h, for c = 2 × 10 -3 . The color scale is in meters. The critical size a c of sills is calculated numerically, as illustrated in Fig. 13. The calculations are done with a constant k = 5 MPa m -1 . Horizontal and vertical white lines locate the cuts shown in Figs. 15 and 16, respectively. E = 100 GPa, ν = 0.35 and ρ = 2000 kg m -3 .

Figure 16 :

 16 Figure 16: Plot of the calculated critical size a c of a sill as a function of h for k = 10 5 Pa m -1 (solid line), k = 5 × 10 6 Pa m -1 (dashed line), and k = 10 10 Pa m -1 (dashed-dotted line). The curve with k = 5 × 10 6 Pa m -1 corresponds to the vertical cut of Fig. 14. The bold black curve represents the analytical solution of [35], which would correspond to k = ∞. ∆P = 5 × 10 6 Pa, E = 10 11 Pa, ν = 0.35, ρ = 2000 kg m -3 , c = 2 × 10 -3 .

- 2 k 1 2 c 57 Table 2 : 9 - 5 -10 9 4 - 10 9

 212572959410 Stiffness of elastic foundation, Pa m -1 l Characteristic length of the system, m n Parameter of magma pressure distribution P Magma pressure distribution in the sill, Pa P 0 Magma pressure at centre of the sill, Pa P a Magma pressure at tip of the sill, Pa P m Overpressure in magma reservoir in the model of [58], Pa Q Injection volumetric flow rate, m 3 s -1 q 0 Weight of the plate per unit surface, Pa r Lateral coordinate, m t Time, s u Horizontal displacement, m u M Horizontal displacement calculated from the model of [58], m w Vertical displacement, m w 0 Vertical displacement before intrusion forms, m w Vertical displacement above the intrusion, m w 2 vertical displacement outside the intrusion, m w i Uplift due to the intrusion, m w imax Uplift at r = 0, m w M Vertical displacement calculated from the model of [58], m w P J Uplift calculated from the solution of [69], m V Volume of sill, m 3 z Distance from the neutral plane of the plate, m α Coefficient of proportionality between V and ∆P in Appendix B, m 3 Pa -1 β Exponent of the power law relationship between P 0 and V ∆P Magma overpressure, Pa ∆P 0 Magma overpressure at centre of the sill, Pa ∆P a Magma overpressure at tip of the sill, Pa r Radial strain ra Radial strain at the tip of the intrusion (r = a, z = h/Radial tensile strain at failure z Vertical strain ν Poisson ratio ρ Density of overburden, kg m -3 Range of geological values. Parameter Range of geological values References a 10 3 -10 5.5 m e.g. [70] E 10 10 11 Pa e.g. [86] E wl 10 Pa e.g. [63, 2] h 1000 -5000 m e.g. [70] h wl 1 -10 m e.g. [43] k 10 Pa m -1

Table 1 :

 1 Units and symbols. Radius of the sill at sill-to-inclined sheet transition, m a m Radius of spherical magma reservoir in the model of [58], m C i Integration constants, variable D Bending stiffness of the plate, Pa m 3 D e Equivalent bending stiffness of a stack of thin plates, Pa m 3 D i Bending stiffness of a thin plate i, Pa m 3 E Young Modulus, Pa E wl Young Modulus of the weak layer, Pa F Force applied at the base of bending plate, N F h Force applied at the base of bending plate for homogeneous pressure, N h Depth of sill, i.e. thickness of the elastic plate, m h i Thickness of a plate i, m h m Depth of spherical magma reservoir in the model of [58], m h wl Thickness of weak layer, m g Acceleration due to gravity, m s

	Parameter Definition and dimension
	a	Radius of the sill, m
	a c

Acknowledgments

This study was supported by a Center of Excellence grant from the Norwegian Research Council to PGP. This research was supported by a Marie Curie FP7-Intra-European Fellowship (PIEF-GA-2009-237089) and a Marie Curie FP7-Reintegration-Grants (PCIG-GA-2011-303871) within the 7th European Community Framework Programme. The authors acknowledge the constructive reviews of two anonymous reviewers and the suggestions of the editor Joan Marti.

Appendix B

The derivation of the expression of the stress intensity factor is based on the mechanical energy U of the system due to the intrusion. It is equal to the work of the overpressure on the roof of the sill as it is inflated:

The expression of the volume of the intrusion for any pressure distribution is given by:

When P a = P 0 , this Eq. 52 leads to Eq. 21. In addition, Eqs. 32 and 33 show that when P a = P 0 , C 1 ∝ ∆P and C 2 -q 0 /k ∝ ∆P . Therefore, Eq. 21 shows that V = α∆P ,