Mauricio M Mora
email: mmmora@geologia.ucr.ac.cr

Philippe Lesage

Fabien Albino

Gerardo J Soto

Guillermo E Alvarado

Continuous subsidence associated to the long lasting eruption of Arenal volcano (Costa Rica) observed by dry tilt stations

Arenal volcano is a small (~1750 m a.s.l., ~10 km 3 ) stratovolcano that continuously erupted between July 1968 and October 2010. During this long-lasting eruption (over 42 years), a large volume of material -about 5.6 x 10 8 m 3 of dense rock equivalent -has been extruded and has produced a thick and extended lava field, mainly on the western flank of the edifice. Measurements of ground deformation obtained using a network of dry-tilt stations are presented for the period 1986-2000. They show a continuous subsidence of the volcano with maximal amplitude on the western side. The load effect of the lava field is calculated and explains the largest part of the observed tilts. Once the data are corrected by this load effect, pressure source models are not supported by the observations and by quality criteria on the models. Although the dry tilt data from Arenal volcano give limited constraint on the deformation models, they are representative of a long period of activity that cannot be recovered by other means. Moreover, the corresponding interpretative model is consistent with results obtained by geotechnical studies and modern ground deformation methods like InSAR.

INTRODUCTION

Ground deformation study is, along with seismology and remote sensing of volcanic gases, one of the most common and widely used disciplines applied to research and monitoring of volcanoes. It can give important constrains on many volcanological problems such as magma storage and transport to the surface and instability of the volcano flanks. Geodetic techniques used on volcanoes have greatly evolved from man dependent methods (levelling, electronic distance measurement, dry tilt) to automatic ones (GPS, electronic tiltmeters), offering possible real-time measurements, and to satellite remote sensing techniques, like Synthetic Aperture Radar interferometry (InSAR), that allow measurement of the overall deformation field. Application of these techniques depends, of course, on economic and logistical constraints. For example, many active volcanoes are located in tropical environments which pose many difficulties to geodetic techniques, either on land or from space, because of the weather conditions (i.e. frequent cloud cover), vegetation and high erosion rates.

One of the first ground deformation studies applied to active volcanoes in Costa Rica, is linked to the onset of Arenal eruption in July 29 th 1968 and the interest of Costa Rican Institute of Electricity (ICE) to monitor the volcano due to its proximity to a hydroelectrical power plant site. Arenal is a small (~1750 m a.s.l., ~10 km 3 ) and young (~7000 years old) stratovolcano that has been continuously erupting since the series of explosions, beginning on July 29 th 1968, that formed 3 new craters (A, B and C, from lower to upper). Those craters together with the former crater D, lie on an E-W fissure system (Fig. 1). From September 1968 to 1974, effusive activity from crater A produced a basaltic-andesitic lava field on the western flank. In 1974, the activity migrated to crater C, forming a new cone on the western flank of the edifice. Since then, there was a lava pool at crater C [START_REF] Alvarado | Pyroclastic flow generated by crater-wall collapse and outpouring of the lava pool of Arenal Volcano[END_REF], where rare and transient lava domes and hornitos (or a combination of both) were formed. Aa lava flows have also been erupted almost continuously from the lava pool, changing to blocky lava flows downslope [START_REF] Borgia | On the Mechanisms of Lava Flow Emplacement and volcano Growth: Arenal, Costa Rica[END_REF][START_REF] Wadge | The magma budget of Volcán Arenal, Costa Rica, from 1968 to 1980[END_REF]. The lava field extruded since 1974 from crater C has partly overlapped the previous lava field produced by crater A. Up to now, Arenal has erupted for over 42 years, totaling ~0.70 km 3 of lava and tephra (~0.56 km 3 Dense Rock Equivalent -DRE), weighing ~1.4 x 10 12 kg, over a surface of about 7.5 km 2 with a maximum thickness of ~150-200 m. Up to the year 2000 (the end of our period of analysis), the volume and weight were about 0.54 km 3 DRE and 1.35 x 10 12 kg, respectively [START_REF] Wadge | The magma budget of Volcán Arenal, Costa Rica from 1980 to[END_REF].

Strombolian activity began in 1984 [START_REF] Borgia | Structural, stratigraphic and petrologic aspects of the Arenal-Chato volcanic system, Costa Rica: Evolution of a young stratovolcanic complex[END_REF][START_REF] Barquero | Arenal Volcano (Costa Rica) Premonitory Seismicity[END_REF]. However, since 1987 it changed to Vulcanian eruptions with infrequent pyroclastic flows originating from column collapse, lava front collapse and lava pool collapse. The major pyroclastic flow events occurred in June 1975, August 1993, May 1998, August 2000, and March 2001[START_REF] Alvarado | Pyroclastic flow generated by crater-wall collapse and outpouring of the lava pool of Arenal Volcano[END_REF] and smaller ones in September 2003, September 2007, June 2008, June-July 2009and March 2010. Since 1998, the activity at crater C changed from very explosive to more effusive. Consequently, this crater grew several tens of meters and the lava pool was replaced by a viscous lava crust with dome-like structures leading to lava flows, up to October 2010. Since then, volcanic activity decreased substantially with practically no explosive events and minimal seismic activity.

The first measurements carried out to detect and observe volcanic deformation at Arenal were performed in 1969 (from September 11 th to December 14 th ) using a diamagnetic tiltmeter that was located 3.5 km NNW of the active crater A [START_REF] Sawdo | Tiltmeter installation at Arenal volcano in Costa Rica: Report Arthur D. Little Inc[END_REF][START_REF] Simon | Preliminary results from tiltmeter recording at Arenal volcano[END_REF]. Tiltgrams showed deflation and inflation on the order of 10 µrad. After that, ICE built four optical leveling (dry tilt) stations along a radial line on the west flank of Arenal. The data obtained from those dry tilt stations for the period 1976 to 1978, were analyzed by [START_REF] Melson | Arenal volcano: Results of dry tilt measurements[END_REF]. This dry tilt network was extended in 1985 by six additional stations installed on the northern and eastern flanks. Although some stations have been destroyed by lava flows, pyroclastic flows or erosion, very long-duration data series have been obtained. The long-term deformation revealed by this network correspond to a deflation of the volcano. This trend has been confirmed by other geodetic studies, such as that of [START_REF] Van Der Laat | Deformaciones asociadas al vulcanismo activo de Costa Rica entre 1981 y 1988[END_REF] from July 1982to 1988. Later, Hagerty et al. (1997) carried out the first continuous GPS measurements at Arenal, observing a shortening of 7.5 ± 0.4 mm/yr (mm yr -1 ) from May 1995 to March 1997 along a N-S baseline.

Two main hypotheses have been proposed to explain the long term subsidence of Arenal. First, following the most common model to explain deflation on volcanoes, [START_REF] Simon | Preliminary results from tiltmeter recording at Arenal volcano[END_REF], [START_REF] Melson | Arenal volcano: Results of dry tilt measurements[END_REF], [START_REF] Van Der Laat | Deformaciones asociadas al vulcanismo activo de Costa Rica entre 1981 y 1988[END_REF] and [START_REF] Hagerty | Observations at Costa Rican volcano offer clues to causes of eruptions: Eos[END_REF] interpreted the deformation as the result of the depressurization of a very shallow magma chamber about 1 km in diameter, with depths between 0.8-1.5 and 4 km. Second, [START_REF] Wadge | The magma budget of Volcán Arenal, Costa Rica, from 1968 to 1980[END_REF], Alvarado et al. (1988[START_REF] Soto | Análisis de la inclinometría seca en el volcán Arenal[END_REF], [START_REF] Mora | Étude de la structure superficielle et de l'activité sismique du volcan Arenal[END_REF], and [START_REF] Alvarado | Estabilidad del cono y comportamiento de la fundación del edificio volcánico del Arenal (Costa Rica)[END_REF] proposed as an alternate explanation, the response to the load produced by the lavas emplaced since 1968. Indeed, post-eruptive deformation related to emplacement of lava fields have been detected by geodetic and remote sensing measurements at Sakurajima volcano [START_REF] Ishihara | Lava flows at Sakurajima Volcano (I) -volume of the historical lava flows: Annals of Disaster Prevention Research Institute[END_REF], Piton de la Fournaise [START_REF] Delorme | Apport des déformations à la compréhension des mécanismes éruptifs: le Piton de la Fournaise[END_REF], and Etna [START_REF] Murray | The influence of loading by lavas on the siting of volcanic eruption vents on Mt Etna[END_REF][START_REF] Briole | Post-eruptive deformation associated with the 1986-87 and 1989 lava flows of Etna detected by radar interferometry[END_REF]. Moreover, [START_REF] Grapenthin | Pressure sources versus surface loads: Analyzing volcano deformation signal composition with an application to Hekla volcano, Iceland[END_REF] demonstrated that it is important to correct the observed deformation from the load effect prior to any estimation of pressure source parameters. Recently, two studies brought new insights on Arenal deformation. [START_REF] Alvarado | Stability of the cone and foundation of Arenal volcano, Costa Rica[END_REF] speculated on an incipient deformation stage generated by spreading of the basement composed of weathered volcanic (mostly epiclastics) rocks strongly lateritized. [START_REF] Ebmeier | Steady downslope movement on the western flank of Arenal volcano[END_REF] interpreted InSAR measurements as a steady downslope movement on the western flank.

In this paper, we present the dry tilt network and the observations obtained from 1986 to 2000. Following the approach of [START_REF] Grapenthin | Pressure sources versus surface loads: Analyzing volcano deformation signal composition with an application to Hekla volcano, Iceland[END_REF], we estimate the deformation produced by the load of lava emplaced from 1988 to 2000 and compare them to the tilt measured in the same interval. Then, by using numerical methods, we investigate the ability of the pressure source model to explain the observed deformations corrected by the load effect. Finally, we open a discussion on the possible origins of the long-term subsidence at Arenal volcano.

DESCRIPTION OF THE TILT STATIONS AND OBSERVATIONS

The dry tilt network ICE adopted the dry tilt technique developed at the Hawaiian Volcano Observatory in the late 1960s [START_REF] Kinoshita | The measurement of crustal deformation related to volcanic activity at Kilauea Volcano, Hawaii[END_REF] and following studies [START_REF] Sylvester | The dry tilt method of measuring crustal tilt[END_REF]Yamashita 1981). The first 4 dry-tilt stations (A, B, C and D) were set up on the western flank of Arenal in October 1974 (Fig. 1). Station A was covered by a lava flow on January 29th 1977, then in 1985, a year after the beginning of the explosive phase, six other stations (E, F, G, I, J and K) were installed on the eastern and northern flanks of the volcano, completing a network of 9 dry-tilt stations.

The stations are composed of benchmarks forming 40-m-side triangles (stations E, F, G, I, J and K) and squares (stations C and D) with a reference point at the center. These last stations consist of two triangular sub-arrays called C1 and C2, and D1 and D2 (see Table 1). The benchmarks are stainless steel rods (2.5 cm in diameter and 3-6 m long) placed in holes cemented with concrete, leaving only 5 cm of the rod above the surface. The top of each rod is rounded and protected to prevent corrosion.

On each flank, the different sets of stations were built following radial lines from the active crater. Volcanic and human activity, though, destroyed part of the network. Station B collapsed by erosion along a contiguous gully in 1991, but serious damage was already visible by late 1990. Station C was damaged by the August 28 th 1993 pyroclastic flow and was repaired some months later. However, one of the benchmarks was eroded and measurements at sub-array C2 stopped in 1995. Station G was unintentionally destroyed by local people in 1999 and station E was covered by the August 23 rd 2000 pyroclastic flow.

Dry tilt data set

Measurements were usually carried out every three to four months using an optical level and a 3 m long Kern invar rod. However, logistical and economic problems sometimes increased this measurement interval and, for some years, only two or one measurement is available. The local tilt is calculated following the method proposed by Yamashita (1981). The standard error of the tilt is about 2 µrad, but frequently 3-5 µrad according to experience in Iceland although larger errors ocurr if the rods are not well coupled to the surrounding rock or soil [START_REF] Tryggvason | Furnas volcano ground deformation measurements[END_REF]. At Arenal, all the stations are constructed on tephra and epiclastic deposits because of the absence of solid rock, which may induce local instabilities.

The total accumulated tilt vectors obtained from the ICE's dry tilt data set from 1988 to 2000, show an overall pattern of tilt down towards the west flank of Arenal volcano, except for station E, which is anomalous in direction and magnitude (Fig. 1). The amplitude of the down-tilt observed is larger at stations close to the summit and those at western flank (Fig 1, Table 1). In fact, the tilt vectors at stations J, G, F, B and D are mainly directed to the NW flank of the volcano. Figures 2 and3 display the temporal variations of the tilt vectors and amplitudes, respectively. It can be observed that the direction of the vectors at stations E and F are quite constant. For stations G, I, J, and K, the tilt directions are also relatively constant, although the signal to noise ratios are much lower. At station C, the directions obtained for the two sub-arrays are similar and constant until 1997. After that, the direction of the vector at C1 is strongly modified. Sub-array D2 also appears to be perturbed, while D1 is stable with a constant tilt direction. At almost all the stations, the tilt amplitude varies approximately linearly with the cumulate volume of magma emitted (Fig. 3) indicating a progressive source process. Hence, the tilt rates tend to progressively decrease together with the extrusion rate. The highest tilt rate is observed at station B until it was damaged in 1990.

Superposed on the long-term down-tilt, there are short-term variations that are not coherent among the stations. Those variations can result from perturbations produced by non-volcanic processes (i.e., charge and discharge of water, temperature changes). However, the lack of complementary continuous deformation and meteorological data as well as detailed knowledge of the response of each station site makes it difficult to discriminate external effects from the volcanic signals. Hence those variations are not considered in this work.

DEFORMATION MODELLING

Quality criteria of models

To evaluate the goodness of models, several complementary criteria are used. The first one is the significance of the models from a physical and volcanological point of view. Unrealistic models must be discarded even if they explain well the data. The second criterion is a misfit function between the observed and calculated tilts. Here we use the error function defined as:

∑ =         - = N i i i o i c t t 1 2 2 σ χ (1)
where N is the number of data (twice the number of dry tilt stations as tilts are twocomponents vectors), i o t and i c t are the observed and calculated tilts, and i σ is the standard error assigned to the datum. This error includes measurement uncertainties, different kinds of perturbation such as monument instabilities, and modelling error due for example to oversimplification of the models. Here we take µrad i 20 0 = = σ σ for all stations and components. In order to compare models that have different numbers of parameters, we use also the Akaike Information Criteria (AIC - [START_REF] Akaike | A new look at the statistical model identification[END_REF] to test their statistical significance and measure the improvement when more parameters are added. The AIC is defined as

L p AIC ln 2 2 - =
, where p is the number of parameters and L is the likelihood of the model:

( )         - ⋅ = - 2 exp 2 2 2 2 0 χ πσ N L
(2).

The AIC can be rewritten as:

( ) 2 2 0 2 ln 2 χ πσ + + = N p AIC
(3).

In our case, where the number of data, N, is small, it is convenient to use a corrected Akaike Information Criteria [START_REF] Burnham | Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach[END_REF]:

1 ) 1 ( 2 - - + + = p N p p AIC AICc
(4).

Mechanical parameters

The numerical models used below require as inputs several mechanical characteristics of the medium such as density ρ, Poisson ratio ν, and Young's modulus E. On the basis of field observations and measurements on samples, the global porosity of the lava field material was estimated in the range 30 to 35%. Thus, taking a density of 2500 kg m -3 for an andesitic dense rock, an average density of 1700 kg m -3 was obtained for the lava field [START_REF] Mora | Shallow velocity structure and seismic site effects at Arenal volcano, Costa Rica[END_REF]. The real medium is a heterogeneous stratified structure. Its elastic coefficients and seismic velocities strongly vary with depth. [START_REF] Mora | Shallow velocity structure and seismic site effects at Arenal volcano, Costa Rica[END_REF] obtained velocity models at Arenal by using the SPAC method [START_REF] Aki | Space and time spectra of stationary stochastic waves, with special reference to microtremors[END_REF]. The estimated S-wave velocities vary from about 300 m 1 s -1 , in the first 10 meters of the structure, to ~1500 m 1 s -1 for depth greater than 300 m. The Young's modulus can be estimated from Vs with the relation

) 1 ( 2 2 ν ρ + = S V E .
Taking ν = 0.25, this yields a range of values for E from 1 to 14 GPa. Since this parameter is poorly known, we will explore the model space using 4 values of E (1.2, 4, 9, and 14 GPa) in the models. These values are consistent with the estimations of E for pyroclastic rocks of Arenal [START_REF] Alvarado | Stability of the cone and foundation of Arenal volcano, Costa Rica[END_REF].

Effect of lava load

As noted above, Arenal has produced large quantities of pyroclastic material and lava that lie on the flanks. [START_REF] Wadge | The magma budget of Volcán Arenal, Costa Rica from 1980 to[END_REF] obtained relatively precise isopach maps for several periods since 1968. The period from 1988 to 2000 was taken for this work because the best data set is complete for this period, both from lava production (Table 2) and deformation.

From 1968 to 1988, 42.6 x 10 7 m 3 DRE were extruded, while 11.4 x 10 7 m 3 DRE were emitted in 1988-2000 (Table 2; [START_REF] Wadge | The magma budget of Volcán Arenal, Costa Rica from 1980 to[END_REF][START_REF] Hofton | Quantifying recent pyroclastic and lava flows at Arenal Volcano, Costa Rica, using medium-footprint lidar[END_REF]. Following the procedure of [START_REF] Grapenthin | Pressure sources versus surface loads: Analyzing volcano deformation signal composition with an application to Hekla volcano, Iceland[END_REF], it is thus possible to estimate and take into account the ground deformations produced by the surface load of lava. We adopted the approach presented by [START_REF] Pinel | Discriminating volcano deformation due to magma movements and variable surface loads: application to Katla subglacial volcano, Iceland[END_REF] to study surface displacements generated by ice load variations at Katla volcano, Iceland. The vertical and horizontal radial displacements of the surface, induced by a unit point mass applied on an elastic half-space are given by:

r E g r U r E g r U r z 1 ) 2 1 )( 1 ( 2 ) ( 1 1 ) ( 2 ν ν π ν π - + - = - = ( 5
)
where g is gravity and r is the distance to the source [START_REF] Sneddon | Fourier transforms[END_REF]. The components of the corresponding tilt vector in Cartesian coordinates are:

       - - = ∂ ∂ - - = ∂ ∂ 3 2 3 2 1 1 r y E g y U r x E g x U z z ν π ν π (6) with 2 2 y x r + =
. These derivatives can be considered as Green functions. In order to obtain the tilt produced by a lava flow load, the Green functions are convolved by the mass distribution:

       ′ ⋅ ′ ⋅ ′ ⋅ ′ - ∂ ∂ = ′ ⋅ ′ ⋅ ′ ⋅ ′ - ∂ ∂ = ∫ ∫ r d r h r r r y U r t r d r h r r r x U r t S z y S z x             ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ρ ρ (7)
where ρ and h are the density and thickness of the lava deposit, respectively, and the integration is carried out over the whole area S of the lava field.

Table 3 shows the misfit χ 2 and AICc obtained for the tilts associated to surface load. The best fit is obtained for E = 1.2 GPa. For higher values of E, the load effect decreases and does not account for the observations. We verified also that lower values of E yield worse fitting. Figure 4 displays the observed tilts vectors and those calculated with E = 1.2 GPa.

Pressure source model

The deformation produced by pressure variations in a spherical source embedded in an elastic half-space can be calculated analytically [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformation of the ground surfaces around them[END_REF]. However, caution has to be taken in presence of topography which may have a strong influence on tilt when the flanks of the volcano are steep. [START_REF] Cayol | Effects of topography on the interpretation of the deformation field of prominent volcanoes -Application to Etna[END_REF] showed that, for slopes greater than 20°, the direction of the tilt can even be reversed and hence, it is important to take topography into account when modeling the deformations. At Arenal, only tiltmeters D, G and K, are located at sites with slope less than 20°.

We calculated the strain and stress perturbations by using the Finite Element Method. An axisymetrical geometry includes an azimuthally average topography of Arenal. In order to avoid border effects, the computing domain is much larger than the volcano. It is a cylinder with radius of 100 km and depth of 100 km including about 10 4 triangular elements. The source is a cavity at depth in which a pressure variation is applied. It is spherical or ellipsoidal. Figure 5 presents the geometry of the model. We compared the results of the FEM calculations with the analytical solution of [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformation of the ground surfaces around them[END_REF] model for plane topography. For distances larger than 2000 m, the two models give similar tilts, while they significantly differ at shorter distance from the crater.

In the case of a spherical cavity, the pressure source model is characterized by 2 parameters, its depth H c and radius R c . The parameter space is explored in the ranges 0 to 5000 m for H c and 100 to 600 m for R c . For each value of the point (H c , R c ), the pressure variation that produces the best fit between the observations and the sum of the tilts associated to the surface load and to the chamber model is estimated. Figure 6 and Table 3 shows the results obtained for the 4 values of E. It can be observed that the misfit does not depend on the radius and only varies slightly with depth. Although minimum misfits are obtained at depths of 500 m below sea level (i.e., 2000 m below the crater) for E ≥ 4 GPa, this parameter is poorly constrained. For E = 1.2 GPa, the depth of the cavity is greater than 3000 m. The optimal pressure variation P ∆ is always negative (depressurization) and strongly depends on the radius. This is similar to the [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformation of the ground surfaces around them[END_REF] model in which ground deformations are proportional to P R c ∆ 3 . In the case E = 1.2 GPa, very low pressure variations are obtained which corresponds to very small surface deformations. Thus, the pressure source barely improves the fit. When E increases, the medium is more and more rigid, the load effect alone no longer account for the observed tilts and a chamber with larger and larger pressure variations must be added to fit the data. The best fitting (χ 2 = 37, Table 3) is obtained for E = 1.2 GPa. The misfit χ 2 increases when E becomes larger. For high Young's modulus, the superposition of the lava load and pressure source models yields better fits than the load effect alone. However, for E = 1.2 GPa, although χ 2 is slightly better, the AICc increases from 122 to 126 due to higher number of parameters. Furthermore, the pressure decrease cannot exceed the lithostatic pressure at the level of the magma chamber. This rough constraint imposes a lower limit on the source radius (Fig. 6). For example, for E = 9 GPa, the radius must be greater than 500 m for an optimal depth of 2000 m below the crater. Thus, from the physical and volcanological points of view, the pressure source model appears to be unrealistic. Some calculations with an ellipsoidal chamber have also been carried out. For an oblate cavity, for example, with horizontal axis larger than the vertical axis, the pressure variation is slightly reduced, allowing a smaller radius, and the data fitting is somewhat improved. However the corresponding AICc is higher due to the increasing number of parameters of the model. Thus, in any of the cases, following most of the defined criteria, a pressure source is not required by the data. In summary, once the surface lava load is taken into account, the existence and effects of a shallow magma reservoir are really not supported by the observed tilts.

DISCUSSION

Our results, based on 12 years of dry tilt measurements, demonstrate that the long lasting eruption of Arenal volcano is associated with a continuous subsidence with larger amplitude on the western flank. This general trend is also confirmed by recent GPS [START_REF] Del Potro | Enhanced three-dimensional stochastic adjustment form combined volcano geodetic networks[END_REF] and InSAR [START_REF] Ebmeier | Steady downslope movement on the western flank of Arenal volcano[END_REF] measurements carried out on shorter time scales. In order to understand the long term evolution of the ongoing eruption, it is thus important to identify the physical processes that can be related to the subsidence.

Considering the large volume and weight of material erupted over 42 years, we first analyze the effect of the surface load of those deposits. Using the precise isopach maps obtained for the period 1988-2000 by [START_REF] Wadge | The magma budget of Volcán Arenal, Costa Rica from 1980 to[END_REF], a density ρ = 1700 kg m -3 and small Young modulus E = 1.2 GPa, we obtain a relatively good fit between the calculated and observed tilts.

The other common source of deformation on volcanoes is the pressure variation in the magma chamber (e.g., [START_REF] Rymer | Volcanic eruption prediction: Magma chamber physics from gravity and deformation measurements[END_REF][START_REF] Lisowski | Analytical volcano deformation source models[END_REF][START_REF] Kohno | Pressure sources beneath Unzen Volcano inferred from leveling and GPS data[END_REF]. Therefore, we also modelled a pressure source embedded in the structure, taking topography into account. We see that for low values of Young's modulus (E = 1.2 GPa), the pressure source does not improve the fit already obtained by the load effect and is discarded by using the Akaike Information Criteria. For higher values of E, data fit is improved by considering the depressurization of a shallow reservoir (2000 m below the crater) with a radius larger than 350 to 550 m. Furthermore, this model requires a large pressure decrease, compared with the lithostatic pressure, which would have induced a strong reduction of the extrusion rate, in contradiction with observations. Consequently, our analysis demonstrates that a deflation of a shallow magma chamber is not the origin of the measured tilts. This result is consistent with geological [START_REF] Wadge | The magma budget of Volcán Arenal, Costa Rica, from 1968 to 1980[END_REF][START_REF] Alvarado | Interpretación preliminar de las deformaciones asociadas al Volcán Arenal (Costa Rica)[END_REF][START_REF] Soto | Análisis de la inclinometría seca en el volcán Arenal[END_REF], petrological and geochemical considerations on a relatively deep chamber [START_REF] Reagan | Changes in magma composition at Arenal volcano, Costa Rica, 1968-1985: Real-time monitoring of open system differentiation[END_REF][START_REF] Sachs | Mafic metaigneuos lower crust beneath Arenal volcano (Costa Rica): evidence from xenoliths[END_REF][START_REF] Streck | Mineral inventory of continuously erupting basaltic andesites at Arenal volcano, Costa Rica: Implications for interpreting monotonous, crystal-rich, mafic arc stratigraphies[END_REF]Williams-Jones et al., 2001).

Using equation 5 and the lava thickness for the period 1988-2000, we also calculated the displacement resulting from the surface load. At about 1800 m west of the crater, close to dry tilt station C, the vertical displacement ranges from 8 cm (using E = 4 GPa) to 26 cm (E = 1.2 GPa). This corresponds to an average subsidence rate of 0.65 to 2.2 cm yr -1 for this period. This rate can be compared to the downward displacement of about 2 cm yr -1 estimated by InSAR data for the period from September 2007 to January 2008 at approximately the same area by [START_REF] Ebmeier | Steady downslope movement on the western flank of Arenal volcano[END_REF]. As the extrusion rate was much higher in 1988-2000 than during the later period, the real subsidence rate was probably larger than that calculated with the surface load. This rough comparison suggests that the load effect alone does not account completely for the observed subsidence. Non-elastic delayed deformation processes can also contribute to it.

For example, a geological and geotechnical study carried out by [START_REF] Alvarado | Estabilidad del cono y comportamiento de la fundación del edificio volcánico del Arenal (Costa Rica)[END_REF][START_REF] Alvarado | Stability of the cone and foundation of Arenal volcano, Costa Rica[END_REF] to characterize and understand the deformation processes of the basement and slopes of Arenal edifice showed that, at the sites of maximum lava thickness, the rate of subsidence would be 60 mm yr -1 , for a cumulated subsidence of ~2 m that represents only 18% of consolidation of the underlying units. They also suggested, according to their calculations, that the total subsidence will be about 19 m once the loading effect ceases (this may happen in one to several centuries according to the consolidation coefficient). At a larger scale, [START_REF] Alvarado | Estabilidad del cono y comportamiento de la fundación del edificio volcánico del Arenal (Costa Rica)[END_REF][START_REF] Alvarado | Stability of the cone and foundation of Arenal volcano, Costa Rica[END_REF] estimated that the Arenal edifice (1.1 km high, 10 km 3 of volume, ~20 x 10 3 kPa) has produced a total subsidence of its basement up to 45 m for the last 7000 years. The large amplitude and long duration of the subsidence associated to the consolidation process suggests that this phenomenon could explain part of the present observations. In particular, the discrepancy between observed and calculated tilt vectors could be due to the consolidation of heterogeneous underlying layers at some stations, located at different active tectonic blocks. [START_REF] Ebmeier | Steady downslope movement on the western flank of Arenal volcano[END_REF] detected a downslope movement of the upper western flank of the volcano, with large horizontal components at the highest parts of the studied zone. They interpreted these deformations as creep along a shallow sliding plane. Along an E-W profile, their interferograms extend from about 800 m to 1900 m west from the crater, including the site of dry tilt station C. At this location, the displacements estimated by InSAR are almost vertical and downward, in agreement with our observations and models. This is also consistent with the studies of [START_REF] Alvarado | Estabilidad del cono y comportamiento de la fundación del edificio volcánico del Arenal (Costa Rica)[END_REF][START_REF] Alvarado | Stability of the cone and foundation of Arenal volcano, Costa Rica[END_REF] who concluded that spreading processes due to dispersion can affect Arenal volcano. It is thus probable that most of the vertical displacements detected by InSAR results from loading effects.

Finally, some discrepancies between the tilt vector calculated with the surface load model and the observed one, especially in their directions (Fig. 4), could result from local perturbations due to: monument instabilities, heterogeneities of the geological structure, proximity of a fault or proximity to a lava flow. Station E, for example, was very close to the N-S fault crossing the edifice (Fig. 1). As the corresponding tilt direction is constant and the amplitude is approximately proportional to the emitted volume, the observed deformation at station E can be interpreted as the result of a perturbation of the general tilt field by the fault. At some stations, nearby lava flows could also produce large misfit if its thickness and geometry are not precisely known.

CONCLUSIONS

The long time series of dry-tilt data presented in this study clearly demonstrate a progressive subsidence of Arenal volcano with maximal amplitude on the western flank. We conclude it is related to the long-lasting eruption, that has produced about 5.6 x 10 8 m 3 of lava and tephra, from which about 1.1 x 10 8 m 3 DRE of material were produced between 1988 and 2000, with a total weight of 0.26 x 10 12 kg. Thanks to the good knowledge of the distribution of this material, the effect of its load could be calculated by using an elastic structure. It appears to be predominant in the observed tilts. Models of pressure source have been superimposed on the load model; however they not improve significantly the fit between observed and calculated tilts and are not supported by the Akaike Information Criteria. Furthermore, although poorly constrained by the data, the corresponding features of a magma chamber -depth, radius, pressure variation -are not consistent with many volcanological data. Other phenomena, such as large-scale sliding, visco-elastic response to loading of the structure and basement of the volcano, settlement of material and local perturbations by structural heterogeneities and faults, have probably some influence on the deformations. Long term deformation measurements after the end of the extrusion would be useful to follow the evolution of subsidence and to confirm, or not, our results and interpretation. Although the dry tilt data from Arenal volcano are of moderate quality and give limited constraint on the deformation models, they are representative of a long period of activity that cannot be recovered by other means. Moreover, the corresponding interpretative model is consistent with results obtained by geotechnical studies and modern ground deformation methods like InSAR. The use of other monitoring technologies, such as GPS, would help identifying the different sources of deformation at Arenal and give a new and wider insight on volcanic processes and hazards.

Williams-Jones, G., Stix, J., Heiligmann, M., Barquero, J., Fernandez, E., Corresponding misfit values between observed tilts and tilts calculated with both load and pressure source models are displayed as black solid lines. Models with pressure variation higher than lithostatic pressure, indicated by black dotted lines, are unrealistic.

Figure 2 .

 2 Figure 2. Path of tilt vectors for three consecutive periods. For C and D tiltmeters, black and gray traces correspond to measurements carried out at C1 and D1, and C2 and D2 subarrays, respectively.

Figure 3 .

 3 Figure 3. Amplitude of cumulative tilt vectors. Bottom panels are: first, effusion rate based on data obtained from Soto and Arias (1998) and Wadge et al. (2006), second, cumulated volume (DRE) of magma calculated from effusion rate estimation.

Figure 4 .

 4 Figure 4. Observed tilt vectors for period 1988-2000 (thin black arrows) and tilts calculated with: A) surface load model alone (thick black arrows), using E = 1.2 GPa, B) combining load and pressure source models (thick black arrows). Vectors obtained from both models are very close to each other due to great similarity of two sets of calculated tilts. Light gray contours: isopach map of the lava flow effused from 1988 to 2000 with contour levels every 20 m.

Figure 5 :

 5 Figure 5: Sketch of pressure source model. Spherical cavity with radius R C embedded at elevation H C over sea level. Axi-symetrical representation of Arenal topography used in mechanical model.

Figure 6 .

 6 Figure6. Pressure source models for different values of Young's modulus. Pressure variation that produces best fitting is plotted as a function of elevation (a.s.l.) and radius of magma chamber. Gray color scale is logarithmic and it is different for case E = 1.2 GPa. Corresponding misfit values between observed tilts and tilts calculated with both load and pressure source models are displayed as black solid lines. Models with pressure variation higher than lithostatic pressure, indicated by black dotted lines, are unrealistic.

  and Duarte Gonzalez, E. (2001) A model of degassing and seismicity at Arenal volcano, Costa Rica; Journal of Volcanology and Geothermal Research, v. 108, p. 121-141. Yamashita, K.M., 1981, Dry tilt: a ground deformation monitor as applied to the active volcanoes of Hawaii: U.S. Geological Survey Open-File Report 81-523. Arenal volcano in its tectonic framework (inset). Main faults labeled: ARF, Aguacaliente River Fault; DF, Danta Fault. Dotted lines: faults covered by the lava field. Craters, 1968-2009 lava field, location of the dry tilt meters and corresponding tilt vectors used in this work. Local coordinates Lambert Costa Rica North.

	Figure captions
	Figure 1:

TABLE 1 .

 1 CUMULATED TILTS FROM 1986[START_REF] Alvarado | Los volcanes de Costa Rica: geología, historia y riqueza natural: San José[END_REF]. Notes: a: B was destroyed in 1990. b: C2 data are not included after its reconstruction in 1994. c: E was destroyed in 2000. d: G was destroyed in 1999.

	Station B	Radial distance to active crater (km) 1.5	Down-tilt amplitude (µrad) 235	Azimuth 52	Duration of observation (months) 45.4 a
	C1 C2	1.9 1.9	253 196	113 81	178.2 106.2 b
	D1	3.1	49	73	178.2
	D2 E	3.1 2.0	103 390	61 286	178.2 170.6 c
	F G	2.7 3.3	129 79	220 237	178.2 158.0 d
	I	2.3	97	308	178.2
	J	2.8	49	245	178.2
	K	3.5	58	326	178.2

TABLE 2 .

 2 VOLUME (DENSE ROCK EQUIVALENTS), WEIGHT AND PROPORTION OF LAVA EMITTED AT DIFFERENT PERIODS (DATA ADAPTED FROM WADGE ET AL., 2006).

	Period	Volume (10 6 m 3 , DRE)	Weight (10 12 kg)	%
	1968 -1988	426	1.07	76
	1988 -2000	114	0.29	21
	2000 -2009	20	0.05	3
	Total	560	1.41	100

TABLE 3 .

 3 MISFIT, χ 2 , AND CORRECTED AKAIKE INFORMATION CRITERIA, AICC, OBTAINED FOR 1) SURFACE LAVA LOAD AND 2) LAVA LOAD + PRESSURE SOURCE MODELS, USING DIFFERENT YOUNG´S MODULUS VALUES. OPTIMAL ELEVATION, HC, AND MINIMAL RADIUS, RC, OF MAGMA CHAMBER ARE ALSO GIVEN.

	Young's Modulus	Lava load			Lava load + pressure source
	E (GPa)	χ 2	AICc	χ 2	AICc	Hc (m)	Rc min (m)
	1.2	42	122	37	126	-3000	200
	4	125	205	51	141	-700	350
	9	161	241	56	146	-500	500
	14	172	252	58	147	-500	550

Mauricio M. Mora -Figure 2 -Figure 2.ai

Mauricio M. Mora -Figure3-Figure3.ai

Mauricio M. Mora -Figure5-Figure5.ai

ACKNOWLEDGMENTS

We are grateful to ICE for supporting the logistical work and the data acquisition. The financial support at different stages of the study was obtained from projects: 113-A4-501, 113-A6-503, 830-A7-511 and 113-B1-230 from University of Costa Rica. Support was also obtained by the European Commission 6th Framework Project VOLUME (Contract N° 018471). We thank Delioma Oramas Dorta for providing the isopach data of Arenal volcano lavafield. We thank Virginie Pinel for useful discussions and for providing some computer programs. We thank the contributions and guidance of Tim Dixon, Geoff Wadge, Rodrigo del Potro and Maurizio Battaglia during earlier stages of the manuscript. Comments and contributions from Glyn Williams-Jones and Hugo Delgado Granados greatly improve the present manuscript. We are also grateful to the Arenal Conservation Area staff for the facilities provided during the field work.

Mauricio M. Mora -Figure 1 -Figure 1.cdr


Julien Lavalle
email: lavalle@in2p3.fr

David Maurin
email: dmaurin@lpsc.in2p3.fr

Antje Putze
email: putze@lapth.cnrs.fr

Direct constraints on diffusion

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

The theoretical understanding of cosmic-ray (CR) transport relies on diffusion of charged particles off magnetic turbulences and has been established for decades [START_REF] Berezinskii | Astrophysics of cosmic rays (Amsterdam[END_REF][START_REF] Schlickeiser | Cosmic ray astrophysics[END_REF][START_REF] Strong | [END_REF]. In this picture, CRs are confined in an extended region that encompasses the Galactic disk, which can be assumed as a homogeneous magnetic cylinder at first order. Therein the diffusion tensor is reduced to a rigiditydependent scalar (homogeneous and isotropic diffusion). Yet, it is only very recently that we have been able to start probing the fine structure of CR phenomenology. With the advent of space experiments like PAMELA [4][5][6][7][8][9], Fermi [10,11], and more recently AMS-02 [12], the physics of Galactic cosmic rays (CRs) has just entered the precision era. CR measurements also provide very interesting probes of exotic physics. In particular, the survey of antimatter CR species may unveil traces of dark matter (DM) annihilation or decay in the Galaxy (e.g. [13][14][15]).

The background to DM searches mostly comes from secondary CRs, i.e. those CRs produced from nuclear interactions between the CR nuclei and the interstellar gas (so-called spallation). This secondary component is used to constrain the CR transport model parameters, as the ratio of secondary-to-primary CR nuclei depends very little on the properties of the primaries at their sources, while very strongly on the transport history [16,17]. The most widely used ratio is B/C [18][19][20][21][22][23][24][25], although other ratios like 2 He/1 H and 3 He/ 4 He are equally powerful [26]. Once the transport parameters are set (from B/C analysis), one can fully predict the fluxes of the other secondary species provided the relevant production cross sections are known. Such calculations have been done for secondary positrons [27][28][29][30], antiprotons [31][32][33][34], and antideuterons [35]. For all but the positron case, for which energy losses play a major role in contrast to nuclei, these computations are poorly sensitive to the theoretical uncertainties affecting the transport parameters in spite of the large degeneracies induced by the B/C analysis [19,20]. Indeed, secondary nuclei experience the same propagation history as boron nuclei.

In 2-zone diffusion models, the most critical uncertainty for DM searches stems from the degeneracy between the normalization of the diffusion coefficient 1 K 0 and half the vertical extent of the diffusion halo L. In-deed, the B/C data mostly constrain the CR escape time ∝ L/K 0 -as all other secondary-to-primary ratios of stable nuclei originating from the Galactic disk. In contrast, DM-induced CRs are produced all over the diffusion halo (and outside), and its size L has a strong impact on the signal predictions: their flux roughly scales like ∼ L 2 /K 0 (assuming a constant DM density, a fairly good approximation for qualitative understanding). This picture is valid whenever the transport is dominated by spatial diffusion.

In Ref. [36], the Authors proposed two extreme configurations to bracket the theoretical uncertainties on the DM signal predictions, dubbed min and max, relying on the B/C analysis performed in [20]. The former (latter) is featured by a very small (large) diffusion zone with L = 1 (15) kpc, and is associated with low (large) signal predictions. In practice, the min model is usually invoked to minimize the antiproton constraints on DM candidates, while max is used to promote detectability -the relative difference between the two almost reaches two orders of magnitude in terms of flux predictions. Such a large range for L strongly affects the antimatter CRs as reliable probes of the DM parameter space. This is particularly important in the light-intermediate WIMP mass range (10-100 GeV), where antiprotons could be used to place severe constraints on WIMPs annihilating or decaying into quarks [37].

There are serious hints that L should be larger than ∼1 kpc, but no stringent bounds so far. Radioactive species are insensitive to L at low energy because their lifetime is shorter than the diffusion time to reach the halo boundary. Using for instance 10 Be/ 9 Be breaks the K 0 /L degeneracy and sets constraints on L [19,21]. However, it was shown in several studies that this method is very sensitive to the modeling of the local interstellar medium (ISM), and strongly affected by the presence of a local underdensity, known as the local bubble [38,39]: this relaxes the lower bound on L, depending on the size of the underdense region [24,40]. There are also other, while more indirect, hints for larger values for L, coming e.g. from calculations of the diffuse Galactic gamma-ray [41] or radio emissions [42][43][44]. Nevertheless, predictions of these observables rely on more ingredients (line-of-sight integrals depending on the astrophysical source, ISM, and/or magnetic field distributions).

In this paper, we propose to use the low-energy secondary CR positrons as an additional direct constraint on L. We exploit the fact that the propagation history of positrons is different from that of nuclei, due to energy losses. This typically shortens the mean free path of positrons, and the dependence on L is milder. In particular, the secondary positron flux roughly scales like ∼ 1/ √ K 0 , allowing us to place a lower bound on L from the current positron data, assuming the B/C-induced relation between K 0 and L. We will only rely on secondary positrons, though it is clear that a primary component is also expected from recent measurements of the positron fraction [4,10,12]. Our approach, suggested in [45], is complementary to the study carried later in Ref. [44], though with a different propagation treatment, for which the main constraints came from the diffuse radio emission data. We first briefly discuss the propagation modeling and relevant parameters, then sketch our statistical analysis method, before going to the results and conclusion.

II. TRANSPORT OF COSMIC-RAY POSITRONS

We wish to constrain small 2D diffusion halo models, with typical vertical sizes of L ∼ 1 kpc. In this context, as observers located at 8 kpc from the Galactic center but still far enough from the radial border located at a distance ∼ 10-15 kpc from us, we can safely neglect the radial escape. We then restrict ourselves to a much simpler 1D problem along the vertical axis. The source of secondary positrons originates in the Galactic disk from spallation processes induced by the primary cosmic rays scattering off the interstellar gas. This source term is well constrained within a radius of ∼ 1 kpc around the Earth, since the cosmic-ray flux barely varies over such a spatial scale, and the average gas density confined in the disk is well estimated. It can safely be approximated to

Q(E, x) = 2 h n ism δ(z) Q 0 (E), where n ism = 1 cm -3
is the ISM gas density, h = 100 pc is half the disk width. The energy dependence is carried by Q 0 (E) = 4 π n ism i,j f j dT (dφ cr,i (T )/dT ) (dσ ij→e + (E)/dE), which convolves the CR flux (species i) with the ISM gas (species j, featuring a fraction f j ).

In the following, we will stick to the formalism presented in [28][29][30] for the calculation of secondary positrons. In [28] (see their Fig. 10), it was shown that transport configurations implying both reacceleration and convection resulted in a prominent low-energy bump around 1 GeV in the secondary positron flux. To be conservative in our comparison with the data, we consider that positrons are only driven by spatial diffusion and energy losses. This actually significantly reduces the computational time. We can define an energy-dependent propagation scale λ from the diffusion coefficient and the energy loss term b(E) = -dE/dt as follows:

λ 2 (E, E s ) = 4 Es E dE K(E ) b(E ) (1) = (3.56 kpc) 2 K 0 10 -2 kpc 2 /Myr τ l 10 16 s × (E/E 0 ) (δ-α+1) α -δ -1 1 -(E/E s ) α-δ-1 , (2)
where E (E s ≥ E) is the observed (injected) positron energy, and where we have used the expression b(E) = (E 0 /τ l )(E/E 0 ) α for the energy losses. For GeV positrons, losses are dominated by inverse Compton and synchrotron processes and the Thomson approximation holds, such that α = 2 and τ l 10 16 s (for E 0 = 1 GeV) to a very good approximation. The propagation scale λ allows to write the Green function for the positron transport:

G 1D (E, z ← E s , z s ) = e -| z-zs λ | 2 b(E) √ π λ 2 z,zs→0 -→ 1 b(E) √ π λ 2 . ( 3
)
For small-L models, one has to account for the boundary effect. When both the observer and the source term are confined into the disk (z, z s = 0), the following series expansions can be used in the regime λ L, most relevant in the low energy range [46]:

G 1D,Helm (E ← E s ) = 1 b(E) L n=+∞ n=1 e -| (2 n-1) λ 4 L | 2 (4) ≈ e -| λ 4 L | 2 + e -| 3 λ 4 L | 2 b(E) L ,
If we express the B/C correlation between K 0 and L as K 0 = κ L, with κ ∼ 10 -3 kpc/Myr [20], we can read off the dependence of the CR positron density on L from the leading term of the series:

G 1D,Helm (E ← E s ) ≈ e -|κ τ l f (E,Es)| 2 b(E) L , where f (E, E s ) ≈ (E/E 0 ) (δ-1) /(1 -δ) ≈ 2 E/(GeV).
Interestingly, the exponential term no longer depends on K 0 and L (because of the cancellation in the ratio λ/L ∝ K 0 /L ∝ κ), and the leading pre-factor is therefore proportional to 1/(b(E) L), where it appears clearly that small values of L lead to large secondary positron fluxes.

So far, we have mostly discussed the roles of K 0 and L. The impact of δ can be understood from the predicted spectral shape for the secondary positron flux [29], which roughly scales as E -γ , with the spectral index γ γ+(1+δ)/2, where γ ≈ 2.7 is the source index associated with secondary positron production. Therefore, the smaller δ the harder the secondary positron spectrum.

Finally, we also emphasize that since we focus on the GeV energy range, our results are sensitive to solar modulation effects. We use the one-parameter force-field model [47] to take it into account, and discuss its effects into more details afterward.

III. STATISTICAL ANALYSIS

Given the measured positron flux φ data (E) and associated statistical error σ φ data (E), we wish to constrain only those transport models which lead to secondary positron fluxes in excess with respect to the data. For each data point at kinetic energy E i , the number of standard deviations in a one-sided hypothesis test of a Gaussian variate, i.e. the Z-score, is calculated:

Z i = φ model (E i ) -φ data (E i ) σ φ data (E i ) , (5)
where φ model (E i ) is the modulated flux estimated for a given parameter set {K 0 , L, δ}. In a subsequent step, the individual p-value p i is estimated only for data points with a positive Z i :

p i = 1 -Φ(Z i ) = 1 -erf(Z i / √ 2) 2 , (6)
where Φ(Z) is the cumulative distribution function of the Gaussian distribution. The independent p i values of a given model are eventually combined into a single test statistic X using Fisher's method [48]:

X = -2 n i log p i . (7)
X follows closely a χ 2 2n with 2n degrees of freedom from which a p-value for the global hypothesis can be easily obtained:

p = 1 - γ(n, χ 2 2n /2) Γ(n) , (8)
where γ and Γ are the lower incomplete and complete gamma functions, respectively.

In the following, we will use a 3-σ exclusion criterion of p < 0.001345, which corresponds to excluding only models that exceed the data by 3-σ or more.

IV. RESULTS AND DISCUSSION

For the predictions of the secondary positron flux, we use the framework defined in [28,29] and associated energy loss parameters. As we disregard reacceleration and convection, a transport model is defined by the normalization of the diffusion coefficient K 0 , its slope δ, and the halo height L, that we vary in the ranges K 0 /L ∈ [10 -3 ,10 -2 ] kpc/Myr, δ ∈ [0.2-0.9], consistent with the B/C constraints [20].

To cover the great variety of possible diffusion model parameters found in the literature, 500 000 parameter sets have been uniformly drawn in the ranges defined above. We have used the PAMELA data [9] associated with a reference solar modulation potential of 520 MV taken from the recent analysis in Ref. [49], and corresponding to the data taking period. We also checked our method with the unpublished AMS-02 data presented at ICRC 2013 (and available from online plots), but will only comment on the trend waiting for the AMS-02 collaboration to confirm their results. A 2D projection in the log (K 0 /L)-L plane of (i) the models parameters drawn, (ii) those leading to positron fluxes not exceeding the PAMELA data (∀ i : Z i ≤ 0), and (iii) those leading to positron fluxes in excess with respect to one or more data points (∃ i : Z i > 0), are shown in left, middle, and right panel of Fig. 1, respectively. The AMS-02 data from ICRC 2013 would lead to the same trends. It can already be seen that the secondary positrons constrain extensively low values of L and log(K 0 /L) independent of the spectral index δ (white area in the middle plot).

The influence of the slope δ can be understood from Fig. 2. In the left (middle) panel, we show the predicted fluxes associated with Z > 0 while not excluded (excluded) by our analysis, for different slices in δ (the Z < 0 models, not shown here, give fluxes that spread over a large area below the data points). We see that the larger δ, the more important the constraint from low-energy data points, as it can be expected from the flux dependence in δ (see Sec. II). In contrast, fluxes predicted from a more gradual slope of δ ∼ 0.3 follow the experimental data more closely (in the GeV range) and are hence equally constrained by the PAMELA data points 1 to 14 ranging from 1.64 to 33.1 GeV. In the right panel of Fig. 2, we display the data points which have individual positive Z-score (Z i > 0) for the models excluded from our nominal analysis, for different slices in δ. The trend explained above is explicit in this plot.

In Fig. 3, we show the 3-σ exclusion curves we obtain in the log(K 0 /L) -L plane. This is the main result of this paper. Plain lines correspond to the exclusion curves associated with different values of δ. Dashed lines show how the results are changed if the secondary positron flux predictions are increased by a global factor of 30%, a very rough way to account for a possible contamination of primary positrons at low energy. In the top panels, we adopt nominal values for the solar modulation, i.e. 520 MV for the PAMELA data taking period, and use all data points to compute the Z-score. In the bottom panel, we take a very conservative viewpoint to secure the results against systematic effects potentially coming from the solar modulation modeling or any other low energy effects, (i) by imposing an analysis threshold of ∼2 GeV (lower energy data points are not considered), and (ii) by significantly increasing the solar modulation parameter φ up to 700 MV.

δ K0 L Vc VA -(kpc 2 Myr -1 ) (kpc) (km s -1 ) (km s -
We report on the same plots the constraints set from B/C as obtained in [20], in the form of bands corresponding to different diffusion slopes δ. There are several reasons as for why most recent studies are not used. First, although powerful statistical tools have been employed since [24][25][26] to sample the most probable regions of the parameter space, we are interested in a wider and more conservative range such as given in [20]. This allows us to include possible systematic uncertainties on the transport parameters that plague their precise determination [24]. Second, the benchmark models (min, med, max) of Ref. [36], widely used in the literature, are based on the parameters found in Ref. [20]. We have placed the points corresponding to the min and med models in both panels of Fig. 3, while the max model is featured by a too large value of L = 15 kpc, outside the scope of this paper. Each point belongs to a B/C band associated with a given slice in δ. The detailed values of the parameters are given in Table I.

An important feature of Fig. 3 is that the positron and B/C constraints are almost orthogonal in the log(K 0 /L)-L plane, underlying their complementarity. This is definitely a strength of this approach as it can complement the one relying on radioactive species, which has completely different systematics. We first concentrate on the left panel, which shows our nominal result. As we can see, the min model is completely excluded. This is an important result per se because this model is very often used in the literature to deal with the theoretical uncertainties affecting the dark matter signal predictions. More generally, models with large values of δ (∼ 0.8), generically associated with small diffusion halos (L ∼ 1 kpc), are excluded by our analysis. Smaller values of δ are fea-tured by larger K 0 /L ratios. For example, we see that for the B/C band that contains the med model (δ ∼ 0.7), the positron data exclude diffusion halos with L 3 kpc, while a less severe constraint of L 2.3 is obtained for smaller diffusion slopes. The dashed lines indicate how these limits would move if a primary positron component contributed an additional 30% to the positron flux at low energy, leading to much larger values of L. Nevertheless, though we know that a primary component is present in the positron data, its origin is still to be confirmed, and its intensity at low energy can hardly be predicted. We emphasize that the local positron flux provides a direct constraint, which relies on less assumptions (e.g. on the distribution of the Galactic magnetic field or the interstellar radiation field) than, for example, limits set from the diffuse radio and gamma-ray emissions, where electrons and positrons are also strongly involved.

We remind the reader that we have already followed a conservative line by neglecting diffusive reacceleration, which would have increased our predictions around 1.5-2 GeV as δ decreases [28]; we have also neglected any primary component. We can still investigate further the systematic effects. As we deal with low-energy positrons, the exact value of the solar modulation level has significant impact on the limit derived, as well as its modeling. The force-field approximation is known to be a simplistic modeling of the solar modulation, but is expected to be accurate within ∼10-20% above 100 MeV [50]. In our nominal analysis, we took a modulation potential of 520 MV, and exploited all PAMELA data points down to 1.64 GeV -the lowest energy data points have strong impact on large δ propagation models. In order to check our results against systematic effects, we adopt the radical option to (i) increase the modulation potential up to 700 MeV and (ii) to remove the first two data points from the analysis, setting a threshold at ∼2 GeV (more precisely 2.38 GeV, 3 rd PAMELA data point). Such an option is extreme, so the corresponding results are to be taken as super-conservative. They are displayed in the right panel of Fig. 3. The limits on L are generically weaker, but the min model remains excluded, as propagation models with large values of δ (∼ 0.8).

V. CONCLUSION

In this paper, we have investigated the low-energy positron constraints ( 10 GeV) on 2-zone cosmic-ray propagation models. We have shown that the resulting bounds are almost orthogonal to the B/C constraints in the L-log(K 0 /L) plane, which makes them particularly attractive. The main result is that we exclude the min benchmark propagation model and more generally large diffusion indices δ ( 0.8). We also strongly disfavor small diffusion halo models with L 3 kpc, the constraint weakening as the diffusion slope δ decreases. This has important consequences for dark matter studies which are often addressed in the frame of 2-zone models. Indeed, this pushes the DM signal predictions toward larger values, which has significant impact on the discovery/exclusion potential of current and future experiments. This will be of particular interest for the searches in the antiproton and antideuteron channels, with AMS-02 and GAPS [51][52][53].

The strength of the proposed analysis, complementary to the B/C or radioactive studies, will significantly im-prove when the AMS-02 data are released (more data points with smaller error bars). Preliminary calculations based on the preliminary AMS-02 data presented at the ICRC-2013 conference already give slightly stronger constraints on L. The next step of this work, that we currently prepare, will be to implement a full study including all low-energy effects and combining the coming AMS-02 data on positrons and PAMELA and AMS-02 data on B/C. We note that the PAMELA collaboration has just released its B/C data [54]. With an analysis based on a very limited number of parameters, the Authors constrain the diffusion slope δ to be in the range ∼ [0.38 -0.42]. However, the choice of the free parameters (wind or no wind, low energy dependence of the diffusion coefficient) has a strong impact on the result, and in particular on δ [24,55]. It is important to perform a full scan of the parameter space with these new B/C data (awaiting for the AMS-02 positron data and final B/C analysis). Better constraints on δ should increase the discrimination power of positrons.

FIG. 1 :FIG. 2 :

 12 FIG. 1: 2D histograms for L versus log(K 0 /L) for δ ∈ [0.2-0.9].

FIG. 3 :

 3 FIG.3: Constraints on propagation parameters in the log(K0/L in [kpc/Myr])-L [kpc] plane. Lines are constraints from the positron flux (separating excluded and allowed regions). Green filled contours are allowed by B/C data. The min and med benchmark models of Ref.[36] (see Tab. I) are indicated by a filled black circle. Left panel: positron contours (from PAMELA data) for a realistic modulation level of φ = 520 MV. Dashed lines correspond to the limits if the secondary positron prediction is increased by 30% to mimic a primary component. Bottom panel (very conservative): Same contours but for a modulation level larger than observed, 700 MV, and taking a very conservative analysis threshold of ∼2 GeV.

TABLE I :

 I Transport parameters (compatible with the B/C analysis) giving the max, med and min p and d dark matter flux.

We will assume K(R ≡ |p/q|) = β K 0 (R/1 GV) δ in the following -where p is the momentum, β the velocity, and q the electric charge.

ACKNOWLEDGMENTS

We would like to thank L. Derome for useful discussions about statistical analysis methods. This work was partly funded by the French ANR, Programme BLANC DMAstro-LHC, Project ANR-12-BS05-0006, and by the Investissements d'avenir, Labex ENIGMASS.