Svetlana Byrdina
email: svetlana.byrdina@univ-savoie.fr

Domingo Ramos
email: dramos@ingemmet.gob.pe

Jean Vandemeulebrouck
email: jean.vandemeulebrouck@univ-savoie.fr

Pablo Masias
email: pmasias@ingemmet.gob.pe

André Revil
email: arevil@mines.edu

Anthony Finizola

Katerine Gonzalez Zuniga

Vincentina Cruz

Yanet Antayhua

Orlando Macedo

K Gonzales Zuñiga

Influence of the regional topography on the remote emplacement of hydrothermal systems with examples of Ticsani and Ubinas volcanoes

Keywords: self-potential mapping, hydrothermal system, hot springs, temperature

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Most active volcanoes are associated with hydrothermal systems, which can be sometimes recognized by fumarolic activity and hot springs on the flanks of the edifice, high electric conductivity at depth, and self-potential anomalies on the surface (e.g., [START_REF] Aizawa | Hydrothermal system beneath Mt. Fuji volcano inferred from magnetotellurics and electric self-potential[END_REF][START_REF] Bedrosian | Hydrothermal circulation at mount st. helens determined by self-potential measurements[END_REF]Aizawa et al., 2009a,b;[START_REF] Revil | Hydrogeology of Stromboli volcano, Aeolian Islands (Italy) from the interpretation of resistivity tomograms, self-potential, soil temperature and soil CO 2 concentration measurements[END_REF]. Variations in hydrothermal circulation may reflect the temporal evolution of the volcanic activity [START_REF] Granieri | Continuous monitoring of CO 2 soil diuse degassing at Phlegraean Fields (Italy): infuence of environmental and volcanic parameters[END_REF][START_REF] Chiodini | Long-term variations of the Campi Flegrei, Italy, volcanic system as revealed by the monitoring of hydrothermal activity[END_REF][START_REF] Carapezza | Active degassing structures of Stromboli and variations in diffuse CO 2 output related to the volcanic activity[END_REF]; in some cases the hydrothermal system interacts with magma ascent leading to hydromagmatic eruptions (e.g. [START_REF] Wohletz | Explosive magma-water interactions: Thermodynamics, explosion mechanisms, and field studies[END_REF]; long term hydrothermal circulation can also contribute to destabilization of the edifice by weakening the stiffness of the rocks (e.g., [START_REF] Reid | Massive collapse of volcano edifices triggered by hydrothermal pressurization[END_REF][START_REF] Reid | Volcano collapse promoted by progressive strength reduction: New data from mount st[END_REF].

Mapping the self-potential (SP) at the surface gives insights into the geometry of the thermal fluid flow inside the edifice. A classical W-shape SP anomaly observed at strato-volcanoes is centered on the summit area and limited by the size of the edifice: as one climbs the slopes of the volcano, the self-potential first decreases several hundred millivolts to several volts. As the summit crater is approached, the self-potential rapidly recovers to the initial level measured on the base of the volcano (e.g., [START_REF] Ishido | Repeated self-potential profiling of Izu-Oshima Volcano Japan[END_REF][START_REF] Ishido | Electrokinetic mechanism for the ,W"-shaped self-potential profile on volcanoes[END_REF]Aizawa et al., 2009a). The self-potential method is sensitive to the movement of the water phase: usually, upwelling thermal waters are associated with a positive self-potential anomaly, and infiltration of cold meteoric water is associated with a negative anomaly, anticorrelating with the topography at the lower part of the edifice. Field data on hydrothermal systems sometimes agree with this simple model but in general show more variable behavior.

For example, at the basaltic shield volcano Piton de la Fournaise in Reunion island, the main hydrothermal system seems to be located within the central Dolomieu crater [START_REF] Lénat | Geolectrical structure of the central zone of Piton de la Fournaise volcano (La Réunion)[END_REF][START_REF] Lénat | The interplay between collapse structures, hydrothermal systems, and magma intrusions: the case of the central area of Piton de la Fournaise volcano[END_REF]. Despite of high precipitation, there is no clear correlation between the topography of its flanks and the self-potential, and no clear manifestation of hydrothermal activity on the flanks or on the base of the volcano. To understand the ground water flow at Piton de la Fournaise, [START_REF] Violette | Hydrologic and thermal modeling of an active volcano: the Piton de la Fournaise, Reunion[END_REF] proposed a system of interconnected perched aquifers channelizing the ground flow to a continuous basal aquifer with an outlet to the ocean.

On strato-volcanoes with steep topography the hydrothermal activity at the summit is often limited to few fumaroles; a water dominated hydrothermal system is found on the lower part of its flanks. Again, the existence of a preferential pathway for thermal water flow is often required to understand the emplacement of the hot springs (e.g., [START_REF] Ingebritsen | Vapor-dominated zones within hydrothermal systems: Evolution and natural state[END_REF]. Some hydrothermal systems extend beyond the edifice [START_REF] Aizawa | Classification of self-potential anomalies on volcanoes and possible interpretations for their subsurface structure[END_REF]Aizawa et al., 2009b;[START_REF] Onizawa | Self-potential distribution on active volcano controlled by three-dimensional resistivity structure in Izu-Oshima, Japan[END_REF]. For example, the hydrothermal system of Iwate volcano (2038 m) is displaced by several kilometers with respect to the summit and to the magmatic source inferred by seismicity and geodetic data (Aizawa et al., 2009a).

The goal of the present paper is to discuss the influence of the regional topography and permeability of the edifice on the geometry of the hydrothermal body defined by the presence of thermal waters. We build our study on field data of self-potential, soil and spring temperature at the surface, and on the results of direct numerical simulation coupling mass and heat transfer inside the volcanic edifice. The basic idea is that there is a control by the topography of the ground water flow, which in turn controls the advective heat transfer from a magmatic source and the associated self-potential distribution. Decrease of viscosity of hydrothermal fluids with temperature leads to increase of the hydraulic conductivity which might explain the existence of the preferential pathways evidenced by the field results. This idea is tested with a new set of data collected at Ticsani and Ubinas volcanoes, major dacitic-andesitic volcanoes in southern Peru. Both volcanoes are characterized by extended hydrothermal systems and are good targets for this study because of the differences in elevation on different sides of the edifice reaching almost 2500 meters at Ticsani and about 2000 meters at Ubinas.

Ticsani volcanic complex

Ticsani volcano is a 5408-m-high dacitic lava dome complex located behind the Central Volcanic Zone in Southern Peru (see position in Figure 1). Ticsani and its two neighbor volcanoes, Ubinas and Huaynaputina are often named as a single volcanic group because of the particular location outside the main arc, in the more complex environment of major regional fault-systems: NWSE strike-slip fault system, associated with the main arc, and N-S normal faults related to the Graben of Rio Tambo [START_REF] Lavallé | Structural control on volcanism at the Ubinas, Huaynaputina, and Ticsani Volcanic Group (UHTVG), southern Peru[END_REF]. The predominant fault orientations around Ticsani volcano are shown in Figure 2 b. Geochemical isotope characteristics suggest that the magmas of Ubinas, Huaynaputina and Ticsani evolved from a common reservoir localized at a depth between 20 and 30 kilometers. [START_REF] Lavallé | Structural control on volcanism at the Ubinas, Huaynaputina, and Ticsani Volcanic Group (UHTVG), southern Peru[END_REF] distinguish two main phases of volcanic activity at Ticsani. The first phase consisted of constructing the edifice of the ancient cone of Ticsani by numerous lava flows and volcaniclastic deposits. This ancient cone collapsed westward leaving a 3 km large ridge in NW-SE direction (Figure 1). In the following phase of volcanic activity, three andesitic to dacitic lava domes (D1-D3 in Figure 1) were emplaced at Ticsani accompanied by explosive eruptions dated to middle Holocene. A sub-Plinian and phreatomagmatic eruption occurred afterwards overlaying these middle Holocene deposits with brown pumice and bombs. Finally, the emplacement of the youngest, Ticsani dome, followed, and the 1600 eruption of Huaynaputina [START_REF] Mariño | Geología, historia eruptiva y evaluación de peligros del volcán Ticsani (sur del Perú)[END_REF][START_REF] Lavallé | Structural control on volcanism at the Ubinas, Huaynaputina, and Ticsani Volcanic Group (UHTVG), southern Peru[END_REF]. This complex history is responsible for the highly asymmetric edifice of Ticsani volcano with its domes sitting in an area marked by a distinct topographic gradient in the NW direction and intersected by a system of active faults.

In 2005, a seismic swarm occurred near Ticsani volcano with a normal fault main shock having a magnitude Mw=5.8, an epicenter close to summit of Ticsani volcano, and an inferred depth of 4 km below the sea level [START_REF] Holtkamp | Earthquake swarms in south america[END_REF]. Ticsani is located quite far from the megathrust system and therefore the authors suggested the regional fault system as a source of the ground deformation [START_REF] Holtkamp | Earthquake swarms in south america[END_REF]. On the other hand, the location of the seismic events below the edifice of the Ticsani complex and reactivation of the fault on its western slope close to the Putina hot springs (see Figure 1) might indicate its possible relation to the activity of the hydrothermal system.

Field Surveys

Over 4000 self-potential measurements were performed around the edifice of Ticsani volcano along six profiles, connecting the base of the volcano with one of the domes, with a spacing between the measurements of 25 m. The profiles were also connected by one closed circular profile with 50 m spacing, at the base of the volcano, in order to create self-potential loops. This approach is traditionally used to avoid the accumulation of errors along unconnected profiles. We used non-polarizable second generation Petiau Pb/PbCl 2 -NaCl electrodes [START_REF] Petiau | Second generation of lead-lead chloride electrodes for geophysical applications[END_REF] and a voltmeter with a sensitivity of 0.1 mV and an input impedance of 10 MΩ. The soil was sufficiently wet to ensure good electrical contacts. The stability of the pairs of electrodes used for mapping was regularly checked by measuring their potential difference which was always smaller than 1.5 mV. The self-potential data were combined to provide a map, which represents the distribution of the electrical potential with respect to the potential of a reference taken as zero. These data were krigged with Surfer Software and the resulting map is shown in Figure 2, the zero reference point was chosen at Toro Bravo lake. The self-potential distribution displays a classical dipolar anomaly with a peaked positive lobe (250 mV) around the youngest dome (Ti), and extended negative lobes on the upper part of the flanks and on the older domes D1 -D3. The negative anomalies reach -3400 mV south of the dome. A more variable behavior is observed in the area of the ancient collapse of the old Ticsani dome, on the NW flank of the present volcano. On this side of the volcano, the mean self-potential values are typically 1000 mV higher than on the SW side.

Soil temperature measurements were performed with a digital thermometer, every 50 m along a circular profile and along four radial profiles crossing the edifice. The measurements were performed at 0.5 m depth to minimize diurnal variations. The relationship between ground temperature and elevation is shown in Figure 3 and in a temperature map corrected of elevation effect which is presented in Figure S4 (Supplementary Material). Interestingly, even after correction of elevation effect, relatively high soil temperatures (Tcorr> 14 • C) were observed around 10 km away from the summit, in W-NW area, where most hot springs are located (see Figure 2). Most of the hot springs are located in the deep river valleys on the W-NW side of the Ticsani complex, while the flanks of the Ticsani edifice present low soil temperature (Tcorr 10 • C). Two fumaroles are located in the summit area and, for altitudes greater than 4800 m asl, the ground temperature data show a systematic thermal anomaly (seen in Figure 3 of the main text and in Figure S4) suggesting the existence of a hydrothermal system at the summit of the edifice of Ticsani.

In addition to these data, we also consider ten springs related to the Ticsani volcano [START_REF] Cruz | Caracterización geoquímica de las fuentes termales alrededor del volcán Ticsani (Moquegua)[END_REF] and located at less than 15 km from the summit. The temperatures of these springs have been measured up to four times a year since [2005][2006]. The springs are located on the W and NW slope of the Ticsani volcano and exhibit temperatures from 10 to 89 • C, i.e. close to the boiling temperature at atmospheric pressure corresponding to 3000 m asl. The temperatures have been recorded with ESPEC RT-12 probe till year 2011 and with TECPEL 318 type K sensor since 2011-2012. Hot spring locations, temperatures with its standard deviation, and values of electric conductivity used below in section 5 can be found in Table S1. Figure S5 shows the temperature versus the elevation for the hot springs together with a linear least-square fit of the linear trend shown by these data. It is of note that the hot spring temperatures anticorrelate with the elevation and distance from the summit, with the exception of the summit fumaroles.

Ubinas

Ubinas which belongs to the same volcanic group as Ticsani [START_REF] Lavallé | Structural control on volcanism at the Ubinas, Huaynaputina, and Ticsani Volcanic Group (UHTVG), southern Peru[END_REF] can be considered as a second example of a volcanic edifice with remote hydrothermal system. Geological setting and detailed description of the recent self-potential mapping of Ubinas can be found in [START_REF] Thouret | Ubinas: the evolution of the historically most active volcano in southern Peru[END_REF]; Gonzales et al. (submitted). As on Ticsani, the edifice of Ubinas is characterized by spectacular altitude differences between the opposite flanks reaching almost 2000 meters (Figure S2); the hydrothermal system is located in the valley of the river Tambo, on the lower (South-East) flank of the volcano. Numerous springs with temperatures varying between 9 and 40 • C are located in up to 20 km distance from the top of the edifice [START_REF] Cruz | Geochemical survey of hot springs around Ubinas, the most historial active volcano in southern Peru[END_REF]. In the present study, we use the data of a NW-SE self-potential profile which crosses the edifice passing through the summit and connecting the two opposite flanks with maximal altitude difference (Figure S2) as well as temperature of one hot spring on the slope monitored with HOBO probe since 2006.

Numerical modeling. Case study Ticsani volcano.

As shown in Figures 23, the existence of hot springs and the distribution of the ground temperature provide a direct evidence for an extended hydrothermal system emplaced westward from the summit of Ticsani volcano. This hydrothermal system extends to the valleys of the Putina and Carumas rivers. We used a 2-D numerical model to check whether the regional topographic gradient is sufficient to influence significantly the distribution of the thermal water and therefore the extension of the hydrothermal body. This model can be used, by trial and error, to constrain a mean permeability value consistent with the entire set of observations: location and maximal temperature of hot springs, distribution of self-potential and ground temperature. At the scale of the Ticsani edifice, ground water seems to be controlled by the gradient of the regional topography rather than by permeability distribution related to the presence of the faults. The observed linear relationships between the hot spring temperature and both the distance from the summit and the elevation (Figure S5) suggest that at the scale of the Ticsani complex, the permeability can be considered as homogeneous. Indeed, in the opposite case, hot spring temperatures would reflect local permeability variations and display more complex patterns.

A review of numerical modeling of magmatic hydrothermal systems was conducted by Ingebritsen et al. ( 2010), and we outline in the following paragraphs some points of their review relevant to our study and conclusions. In general, processes occurring in the hydrothermal system satisfy coupled partial differential equations describing fluid circulation and the corresponding equations of state. The simplest approach is to consider a single-phase liquid water, and the Boussinesq approximation (density is constant except in the term which affects the gravitational forces and responsible for buoyancy). Despite of strong assumptions, this formulation allowed to obtain several realistic conclusions. For instance, Norton and Knight [1977] showed that advective heat transport would not be significant for host rock permeabilities < 10 -16 m 2 which was confirmed with more complex models (e.g. [START_REF] Hurwitz | Groundwater flow, heat transport, and water table position within volcanic edifices: Implications for volcanic processes in the cascade range[END_REF]. The main drawback of the single-phase fluid modeling for the current study is the impossibility to take into consideration the position of the water table and vadose zone.

There exist several widely used numerical simulators taking into account multiphase flow (e.g., [START_REF] Pruess | TOUGH2 users guide, version 2.0. Rep[END_REF]. These simulators permit to calculate the position of the water table and the fluid transfer near the magma body which is not possible with full saturation models. On the other hand, the multi-phase equations include a much larger number of parameters than single phase flow equations (storage term, relative permeability, residual saturation). Some of these parameters, e.g. relative permeability, are not well constrained and give rise to potentially strong non-linear behavior of the governing equations.

Model geometry

We use COMSOL Multiphysics V3.5 to analyse a 2D finite-element model of coupled transient water and heat transport and a post-processing code resolving Laplace equation for electric potential considering partial saturation, a non-reactive flow, and temperatures not exceeding 200 • C. We believe that these assumptions will not lead to significant errors, because the main question of our study concerns the distribution of the water phase outside the summit region rather then the processes related to the degassing of the magma conduit. We postulate the existence of a magmatic body (used as a heat source) below the youngest dome in agreement with the elevated ground temperatures measured in this area (Figure 3). We place the heat source at a depth of 10 km below the surface which corresponds to the inferred depth of the major events in the 2005 seismic swarm [START_REF] Holtkamp | Earthquake swarms in south america[END_REF]. Our global groundwater flow model is considered in transient conditions times exceeding 30 000 years. When surface temperature approaches boiling temperature, simulations are stopped even if steady state is not reached. Fluid flow is due to the gradient of hydraulic heads (gravitational flow) and the free Rayleigh-Benard convection. We compare results of this model to hot spring data, electric conductivity of water, and self-potential data. Neither the permeability nor the electric resistivity structure of the Ticsani is known, therefore we make the following assumptions on their distribution. 1) The permeability decreases with depth everywhere according to Equation 8and Figure 4 b. 2) The electric conductivity is a function of temperature and saturation only (see section 5.7). Two distinct geometries will be considered: model A consists of a basal body and a homogeneous edifice, and model B has one more unit representing a magmatic conduit (Figure 4). This conduit (unit3 in Figure 4 b) is a permeable body embedded in low-permeable matrix. We consider an anisotropic permeability for unit3: vertical permeability is an order of magnitude higher than the horizontal permeability.

The numerical model is aligned along profile 1 ("Pr1-Ti-Pr1" at Figure 2). It passes through the summit of the volcano and corresponds to a high difference of elevation on opposites sides of the volcano. The profile was chosen to avoid the heterogeneous area of collapse at NW of the summit where ground temperature, self-potential data, and topography data display small scale variations related to local heterogeneities (see Figure 2).

Hydraulic problem

We use the Oberbeck-Boussinesq approximation for the the mass balance in porous medium (Nield and Bejan, 2006, e.g.). Under conditions of partial saturation, the hydraulic problem can be solved with the help of the Richards equation and van Genuchten parametrization of the capillary pressure and water retention [START_REF] Richards | Capillary conduction of liquids through porous media[END_REF].

[Ce + S e S ] ∂H ∂t + ∇ • [-K∇ (H + z)] = 0 (1)
where z is the elevation, H = p/(ρg) is the pressure head (m), p is pressure (Pa) Ce = φ ∂S w /∂H is the specific moisture capacity (in m -1 ). S e and S w are the effective and the relative water saturation (dimensionless) related to each other by S e = S w -S r w / 1 -S r w , where S r w is the residual saturation of the wetting phase. φ is porosity. S is the specific storage coefficient (m -1 ), related to matrix and water compressibilities α m and α w by S = ρ 0 g (α m + φ α w ). The hydraulic conductivity K is a function of the permeability at saturation k (in m 2 ) and the relative permeability at a partial saturation k r (dimensionless):

K = k r k ρ f g η f (2)
Here η f denotes the dynamic viscosity of the pore water (in Pa s). According to the van Genuchten model, the porous material is considered as being saturated for pressures exceeding atmospheric pressure (H ≥ 0). The relation between the water saturation and the capillary pressure is defined by:

S e =        [1 + |αH| n ] -m , for H < 0, 1, for H ≥ 0, (3)
Relative permeability is defined by:

k r =        S L e 1 -(1 -S 1/m e ) m 2 , for H < 0, 1, for H ≥ 0, (4)
Here α, n, L, m are empirical parameters describing soil properties (see Table 1).

Boundary conditions for hydraulic problem

Equation ( 1) is solved for the following boundary conditions: a hydrostatic pressure on both sides of the system. We also consider a hydrostatic fluid pressure at the bottom of the system and a constant water flux on the surface corresponding to the effective meteoric recharge. A yearly precipitation typical for the area is about 300 mm [START_REF] Pouyaud | Variabilidad pluviométrica, a escalas anual y cuatrimestral, en la vertiente peruana del Océano Pacifíco[END_REF]. We assume that the effective water input is given by 50% of the yearly precipitation (150 mm), as assumed in previous numerical studies [START_REF] Ingebritsen | Rates and patterns of groundwater flow in the Cascade Range Volcanic arc, and the effect on subsurface temperature[END_REF]. The water table position at the extremities of the Pr1 profile is given by springs (3700 m asl) and Toro Bravo lake (4600 m asl).

Heat transfer equation

The partial differential equation describing the transfer of heat by conduction and advection is given by

∇ • (λ b ∇T ) -c l ρ l u • ∇T = ρ b c b ∂T ∂t , (5)
where λ b is the bulk thermal conductivity (in W m -1 K -1 ) and c b denotes the specific heat capacity (in J kg -1 K -1 ). Bulk properties λ b , ρ b , and c b are calculated using volume fractions of solid and liquid at partial saturation:

λ b = (1 -φ)λ s + φ S w λ l (6) ρ b c b = (1 -φ)ρ s c s + φ S w ρ l c l (7)
At the ground surface, we impose as initial condition a constant temperature T = 10 • C corresponding to the average ground temperature.

Boundary conditions for heat transfer equation

Boundary conditions are convective flux everywhere (-n•(-λ b ∇T ) = 0). Other boundary conditions are convective heat flow on the sides and the bottom of the system and a constant temperature at the surface of the magmatic body. The dependent variables in the Darcy equation are the two components of the Darcy velocity, u x and u z , as well as the pore fluid pressure p, and the temperature T in the heat flow equation. Thus, the system of Equations (1-5) is used to compute the distribution of the Darcy velocity and the ground temperatures for our 2-D geometry.

Equations of state

We assume that the permeability decreases with depth as a result of compaction. Ingebritsen and Manning (1999) derived a power dependence of permeability on depth: log 10 (k) = -14 -3.2 log(z), where z is measured in km, for crustal depths down to 40 km. For the first few kilometers of depth, an exponential relation was proposed instead by [START_REF] Saar | Depth dependence of permeability in the Oregon Cascades inferred from hydrogeologic, thermal, seismic, and magmatic modeling constraints[END_REF] and [START_REF] Garibaldi | Thermal anomalies and geological structures in the provence basin: Implications for hydrothermal circulations at depth[END_REF]. Therefore we assume:

k = k 0 exp -z z c . ( 8
)
where we take z c = 4000 m. The constant k 0 is used as free parameter. For the depth of unit1 of 8500 m, equation ( 8) corresponds to a decrease of density by about one order of magnitude. The dependence of water viscosity η on temperature and pressure can be approximated according to [START_REF] Likhachev | Dependence of water viscosity on temperature and pressure[END_REF]:

η = η 0 × exp a p + E -b p R (T -T 0 -c p) , (9)
For E = 4.753 kJ/mol, η 0 = 2.4055 × 10 5 Pa s, T 0 = 139.7 K , a = 4.42 × 10 -9 Pa -1 , b = 9.565 × 10 -6 J • mol -1 Pa -1 , and c = 1.24 × 10 -6 K/Pa. The viscosity values calculated with equation ( 9) over a temperature interval of 273 -463 K and at pressures of 1 -250 bar are in good agreement with the tabulated data [START_REF] Likhachev | Dependence of water viscosity on temperature and pressure[END_REF].

The decrease of water density with temperature can be described by following equation [START_REF] Graf | Effect of viscosity, capillarity and grid spacing on thermal variable-density flow[END_REF]:

ρ f = 1000 (T -

Electrokinetic problem

Under assumption of a homogeneous distribution of the electric resistivity with spatial coordinates, the analysis of electric potential is completely dependent on the results of the groundwater flow and heat transport simulations. Therefore, the model B results presented in (Figure 6) will be used in order to constrain the electrokinetic modeling. The total electric current density j is given by [START_REF] Linde | Streaming current generation in two-phase flow conditions[END_REF][START_REF] Revil | Electrokinetic coupling in unsaturated porous media[END_REF]:

j = -σ(S w )∇ϕ + Q v S w u, (11)
where u = -K/η f ∇H, ϕ is the electrical potential, σ is the (temperature and water saturation dependent) electric conductivity, Q v -is the macroscopic excess of electrical charge per unit volume of the pore space that can be mobilized by the flow of the pore water (expressed in C m -3 ), u is the Darcy velocity. Once the Darcy velocity components u x and u z are determined and the temperatures calculated with equations (1) to ( 5), the electrical potential at the surface can be calculated using equation ( 11) together with the continuity equation ∇ • j = 0. The excess charge Q v and the electric conductivity σ are material properties and are required to compute the electrical potential distribution from the solution of the ground water flow problem. The excess of charge Q v is related to the permeability k of the medium (equation 8) and can be estimated using an empirical relation derived by [START_REF] Jardani | Tomography of the Darcy velocity from self-potential measurements[END_REF]:

log(Q v ) = -9.2349 -0.8219 log(k). ( 12
)
To date there is no data available of the electrical conductivity distribution below the Ticsani volcano, but an average value can be estimated using the conductivities σ lab of the spring waters (listed in Table S1). The electric conductivity σ lab increases with the temperature of the springs because of the higher mineralization of the pore water at higher temperatures. In addition to this mineralization effect, there is also a direct relationship between the electrical conductivity of water at a given salinity and the temperature (e.g., [START_REF] Matthess | The Properties of Groundwater[END_REF]Revil et al., 1999b):

σ(T ) = σ T 0 (1 + α(T -T 0 )), ( 13
)
where α is approximately 0.02 K -1 . This increase is due to the decrease of viscosity of the pore water with the increase of temperature and its effect on the mobility of the ions in the pore water. The σ lab values were therefore corrected using Equation 13. The resulting relationship between σ(T ) and temperature allows to determine the total temperature coefficient α tot ≈ 0.05 K -1 accounting for both ion mobility and mineralization effects (see Figure S7).

We neglect here the effect of surface conductivity along the pore water-mineral interface (Revil et al., 1999a). The (temperature dependent) conductivity of the rock σ m is related to the conductivity σ w of the pore water by:

σ m = σ w F -1 S n 1 w (14)
where F = φ -m 1 (Archie's law) is a formation factor relating the total connected porosity θ and a pore shape/grain shape coefficient called the cementation exponent m 1 . S w is the relative water saturation determined in Richards equation, and n 1 is constant close to 2. For an average spring water conductivity value at 10 • C, σ w = 0.04 S m -1 , the cementation exponent by default is m 1 = 2. With φ = 0.35, we obtain for full saturation (S w = 1) σ m ≈ 5×10 -3 S m -1 .

Finally, Equation ( 11) is solved assuming that the excess of charge per unit of volume Q v is given by expression (12) and the conductivity depends only on saturation and temperature as given by equation ( 13).

At all the external boundaries of the system, we used an insulating boundary condition for the electrical potential,

n • J = 0. ( 15
)
6. Comparison with the observations

Heat transfer and model calibration based on hot springs temperature data

The equation system (1-5) relates mass transport with heat transport which both depend on the permeability of the medium. The permeability of the unit1 together with the boundary conditions, are therefore the key parameters controlling the geometry of both ground water and heat flows. As it has been shown in numerical studies, the permeability together with precipitation impose the position of the water table [START_REF] Hurwitz | Groundwater flow, heat transport, and water table position within volcanic edifices: Implications for volcanic processes in the cascade range[END_REF]. In the following, we consider the results of the numerical model described above for different values of permeability which are consistent with yearly precipitation at Ticsani and water table position at the extremities of profile. The purpose is to understand the influence of permeability on the thermal flow behavior and to constrain its average value.

We first consider a simple homogeneous model A (Figure 4). For low surface permeability (k 0 = 4 × 10 -15 m 2 in Model A1 in Figure 5 a), we find quasi horizontal isotherms, regardless of the regional topographic gradient. We obtain similar results with higher heat source temperatures and at longer times. Heat source temperatures up to 200 • C were tested and run times up to 100 000 years. For these values of permeability and precipitation, the volcanic edifice is fully water saturated. At full saturation, the hydrostatic pressure is maximal below the summit at a given depth. The thermal flow deviates from the summit to flanks in order to escape from high hydrostatic pressure below the summit even with permeability allowing advective heat transfer.

For higher permeability (k 0 = 8 × 10 -14 m 2 , see Model A2 in Figure 5 b), the cone is de-saturated and pressure distribution does not depend on topography contrasts in the dome area (-4 < x < 2 km). In this case, the regional velocity field is defined by the altitude difference between the opposite flanks. The hydrothermal flow rises to the surface due to convection but is derived by the regional flow to the West flank where a thermal anomaly is created close to the surface. In this model realization, hydraulic conductivity at maximum exceeds the initial values at the same depth by half of the order of magnitude, creating a hydrothermal zone of high Darcy velocities and high temperature at a lateral distance of more than ten kilometers from the heat source (highlighted area in Figure 5). The model realization accounts for the position of the hydrothermal area located more than 10 km from the summit. The PA1 and PA2 hot springs are at boiling temperature and belong to a large hydrothermal area of Putina river ("hot river" in local language).This area is the most important manifestation of hydrothermalism associated with Ticsani volcano. We assume in the following that Putina valley and more exactly, PA1-PA2 hot springs, correspond to the maximum of hydraulic conductivity in model A2.

In order to estimate which mean permeability fits best the observations, we determine the hydraulic conductivity maximum for each permeability value, and compare their altitude and distance to the summit with those of PA1-2 hot springs. Figure S6 shows the result of this sensitivity analysis which estimates the optimal value of surface permeability k 0 ≈ 10 -14 m 2 minimizing the deviation between observations and location of thermal area given by the model A.

Besides of information given by temperature of the hot springs, the model can be constrained using the information on water table elevation. Water table elevation is a key parameter defining the distribution of the electric potential, especially under assumption of homogeneous electric conductivity at saturation. Therefore, it can be constrained by self-potential observed in the field (see section 6.2). The data constraints on temperature and water table elevation require low vertical and high horizontal permeabilities, which can be satisfied with a greater depth of the heat source or by anisotropic permeability. The depth of the heat source is constrained by the seismicity data, therefore we assume that the horizontal permeability of the unit1 is twice as high as the vertical one at the same depth k x = 2 k 0 . This anisotropy ratio is kept constant at all presented models. Physically, the anisotropic permeability can be explained by higher horizontal permeability of quasi parallel layers of pyroclastic deposits.

The simple geometry of the model A does not allow to explain neither elevated ground temperatures nor the selfpotential anomalies in the summit area. To account for this feature, we introduce a conduit with a larger vertical permeability in the center of the model geometry, as shown in Figure 4. We chose the permeability of the conduit in z -direction twice the permeability of the unit1 to satisfy the condition of two distinct thermal anomalies: in the summit area and on the West-flank. In x-direction, the permeability of the conduit is an order of magnitude lower than the permeability of the unit1 as discussed in section 5. Resulting simulation of thermodynamic model is represented in Figure 6 and leads to essentially the same conclusion on the possible permeability range.

In order to compare the model with temperature and locations of the hot springs, we check temperatures at several arbitrary points in the subsurface of the model geometry which we call control points in the following (Figure 4 b). The control points are chosen on the Western flank were the thermal anomaly is encountered. We place the control points below the surface, assuming a rapid adiabatic rise of the ground waters towards the hot springs. The control point Pt1 corresponds to PA1 hot spring. Therefore, the model run is stopped when the temperature at Pt1 reach the temperature observed at PA1, which is close to the boiling temperature. Figure 7 shows the temperature versus distance variation for different values of permeability as well as location and temperature of the hot springs, with control points located 400 m below the surface. For permeability value k 0 = 8 • 10 -15 m 2 , our numerical model matches well with the observed trend while higher and lower values of the mean permeability are not compatible with the observations.

Ground water flow pattern and self-potential data. Permeability constraints given by self-potential data

Once the mean permeability value is obtained, we can determine the mean value of the excess of charge Qv from the empirical equation ( 12) with permeability value from best model B2.

In order to simulate the flat central part of W-shape anomaly of self-potential, we imposed no-flow boundary condition in the interval -2000 < x < 0 m. This boundary condition imposes zero Darcy velocities of the water phase in the crater area and can have different interpretations: it can stand for a self-sealing zone where hydrothermal alteration closed the pores like it seems to be the case in Fossetta area of Stromboli volcano [START_REF] Revil | Hydrogeology of Stromboli volcano, Aeolian Islands (Italy) from the interpretation of resistivity tomograms, self-potential, soil temperature and soil CO 2 concentration measurements[END_REF]. Alternatively, zero Darcy velocity could be the consequence of high gas saturation leading to water de-saturation. Yet another reason for positive central anomaly of self-potential can be a shallow conductive zone due to intense fluid circulation and clay-rich hydrothermal alteration [START_REF] Ishido | Electrokinetic mechanism for the ,W"-shaped self-potential profile on volcanoes[END_REF]. In this case, even with non-zero Darcy velocities and streaming current, the electric potential will be close to zero. Discussion of processes occurring in the vicinity of the crater area and resulting distribution of electric potential is out of the scope of the paper and can be found in [START_REF] Ishido | Electrokinetic mechanism for the ,W"-shaped self-potential profile on volcanoes[END_REF]; [START_REF] Aizawa | Classification of self-potential anomalies on volcanoes and possible interpretations for their subsurface structure[END_REF]; Aizawa et al. (2009a).

Using the value of Q v and the values of the Darcy velocities, we solve Equation ( 11) for self-potential with σ T 0 defined by equation 13 in section 5.7. Comparing the solution of Equation ( 11) with observed self-potential data allows us to constrain the range of possible permeabilities. Figure 8 shows distributions of hydraulic conductivity, temperature, electric potential and electric conductivity of the best model B2 together with observed self-potential data. All model variations of the electric conductivity are due to saturation or temperature variations. At depth, in conditions of full saturation, the electric conductivity is determined by temperature only. On the slopes of the dome, the decrease of electric conductivity is mostly due to partial saturation. This zone corresponds probably to a zone of high resistivity on the slope of volcanic edifices, revealed typically by electric resistivity tomography (e.g. [START_REF] Finizola | Fluid circulation and strcutural discontinuities inside Misti volcano (Peru) inferred from self-potential measurements[END_REF]. Despite the simplicity of the structure of electric conductivity, this model agrees quite well with observations outside the summit area.

The most prominent signal in self-potential data, i.e. the two negative lobes on both sides of the summit, corresponds in our model to a zone of partial saturation below the cone. It is well known that the self-potential method is sensitive to the water table position (e.g [START_REF] Jackson | Regional self-potential anomalies at Kilauea volcano[END_REF][START_REF] Aubert | Internal structure of the Merapi summit from self-potential measuruments[END_REF][START_REF] Naudet | Relationship between self-potential (SP) signals and redox conditions in contaminated groundwater[END_REF]. As water table elevation depends mainly on permeability, it is possible to use self-potential data as an independent constraint. Figure 8 a in Supplementary Material displays the calculated electric potential for different permeabilities together with field data for k 0 = 4 × 10 -15 m 2 and k 0 = 8 × 10 -15 m 2 . The transition between fully and partially saturated zones corresponds to a slope change in self-potential data (compare Figure 6 a with c or 6 b with d). In B1 model, the location of two negative lobes of electric potential does not fit the anomaly of self potential observed in the field reflecting the water table position which is, at a given precipitation, defined by the permeability. In B2 model, the fit between model and observed data is better. A slight discrepancy in the East flank might be a consequence of 2-D modeling: in 2-D geometry, unit3 with its smaller horizontal conductivity represents an impermeable boundary across the regional flow which can be avoided by 3-D modeling.

Application to self-potential data collected at Ubinas Volcano

In the following, we apply the same type of model to self-potential data from Ubinas volcano. The self potential profile shows an anomaly with two asymmetric negative lobes similar to the profile across Ticsani (Figure 9). Applying the same kind of model, we obtain the best fit between the data and observation for k 0 = 6 × 10 -15 , shown in Figure 9. Temperature in a control point of the model is in agreement with the measured (by HOBO probe) temperature of the hot spring Termal Ubinas (29 • C), situated in 3300 m asl at the SE flank of the volcano.

The asymmetry of the two negative lobes correspond to the meteoric water infiltration on the opposite flanks of partially saturated edifice (Figure 9 a). This asymmetry reflects the asymmetry of the water flow system resulting from the elevation difference between the extremities of the profile. The ratio of the amplitudes of both minima depends in our model mainly on the thickness of partially saturated zone and can therefore be used to constrain the horizontal permeability of the edifice. Similar to self-potential at the surface, the distribution of the total current density (Figure 9 a) shows two asymmetric current cells: a large current cell westwards from the summit whereas on the opposite side the current cell is hardly recognizable. A similar result is obtained for Ticsani volcano (Figure 8 c). We believe that in a 3-D geometry, the fit between the model and data will be better and will allow quantitative results. This analysis implies that for certain mean permeability values, the regional topographic gradient is able to significantly divert the thermal water flow. Indeed, our model results give a rather narrow interval of permeability values consistent with the hot spring temperatures and locations (see Figure S6). The locations of fumaroles, hot springs, ground temperature and self potential anomalies which evidence the hydrothermal system and water table position, contain therefore information on the average permeability of the edifice.

Conclusion

With the purpose to understand the emplacement of remote hydrothermal systems, we considered a direct 2-D numerical model of coupled mass and heat transfer together with a post-processing code resolving Laplace equation for electric potential. We took examples of Ticsani and Ubinas volcanoes in Southern Peru, both of them characterized by remote hydrothermal systems and pronounced regional topography gradients.

Our analysis of these two volcanic examples suggests that the regional topographic gradient is able to significantly divert the thermal water flow if the mean permeability value belongs to a certain quite narrow interval. The numerical simulations show that there is a balance between the advection produced by the heat source at depth and the regional flow depending on elevation difference and permeability. If the permeability is high enough to ensure an efficient advective heat transfer and a de-saturation of the edifice, the main hydrothermal system is emplaced above the heat source even in presence of regional hydraulic flow. If the permeability is low enough that the volcanic cone is close to full saturation, a deep hydrothermal system is created with little manifestations at the surface. A remote hydrothermal system is created for intermediate permeability values, which for Ticsani volcano belong to the interval 6 × 10 -15 < k < 10 -14 . Even if initially the permeability is homogeneous, a zone of high hydraulic conductivity develops inside the thermal flow zone due to decrease of viscosity with temperature. Apart from permeability and for a given topography, the lateral distance between the heat source and the surface hydrothermal system depends on the depth of the heat source. The asymmetric negative lobes of self-potential anomaly and the ratio of their amplitudes can evidence the regional water flow system and can be used to constrain the permeability in the summit area.

The present findings have some implications for volcanic hydrothermal systems but also for hydrothermal systems in a tectonic context. An example could be the Syabru-Bensi hydrothermal system in Central Nepal which is located at the foot of the High Himalaya, in the deep valley of Trisuli river, and at a lateral distance of around 20 km from a probable heat source, a zone of high micro-seismicity and decarbonization associated with Main Central Thrust zone [START_REF] Perrier | High carbon dioxide and radon-222 gas exhalation at the Syabru-Bensi hot springs in Central Nepal[END_REF].

In the future, we plan to develop a joint approach in inverting simultaneously EM data (Controlled source magneto telluric and magneto-telluric data), geological data, and both self-potential and thermal data to constrain the geometry of the hydrothermal system and to model the pattern of ground water flow.

, [kg m -3 ] 10 3 Thermal expansion, β, [K -1 ] 2.8 × 10 -4 Thermal conductivity (water), λ l , [W m -1 K -1 ] 0.6 Thermal conductivity (rock), λ s , [W m -1 K -1 ] 2.5 Specific heat capacity (rock), cp r , [J kg -1 K -1 ] 8 × 10 2 Specific heat capacity (water), cp w , [J kg -1 K -1 ] 4.2 × 10 2 Electrical conductivity at 25 • C, σ T 0 , [S m -1 ] 6 × 10 -3 Electrical conductivity, σ(T, S w ), [S m -1 ]
Eq. 13 Excess of charge per unit volume, Q v [C m -3 ] Eq. 12 Permeability versus depth, k(z)

Eq. 8 Relative permeability, k r Eq. 4 Dynamic viscosity, η(T ), [Pa s] Eq. 9 Saturation exponent, n 1 1. Hydraulic conductivity is shown by color, isolines of temperature are given by thin gray scale lines. Arrows indicate velocity field. a) For surface permeability k 0 = 4 × 10 -15 m 2 , the cone is completely saturated, and hydrostatic pressure below the cone prevents thermal flow to rise. No temperature anomaly is encountered close to the surface, 25 • C isoline is at 3000 m depth. b) For surface permeability k 0 = 2 × 10 -14 m 2 the cone is almost de-saturated. Fluid flow in the entire model geometry is controlled by regional topographic gradient. Advective flow rises to shallow depth and is derived towards the lower West flank by this regional flow, creating temperature anomaly at the surface at 10 km distance from the summit. Significant increase of hydraulic conductivity due to viscosity decrease corresponds to a zone of intense flow of thermal water. The flow is derived towards the lower part of the flank but reaches the surface. Model B2 is the only model that agrees with following observations: the temperature anomaly at more than 10 km from the volcano summit revealed by the hot spring data, the distinct thermal anomaly at the summit confirmed by soil temperature measurements, and the position of the water table given by the self-potential data (see text). a) Model B3: for surface permeability k 0 = 8 × 10 -14 m 2 , the cone is de-saturated, similarly to A2 model. The thermal flow is directed almost vertically towards the summit.

Figure 1 :

 1 Figure 1: Location of the studied area. a and b) Location of Ticsani and Ubinas volcanoes in Central Volcanic Zone in Southern Peru. Gray dashed line indicates volcanic arc, white triangles show volcanoes. c) Simplified geological map of Ticsani volcano after Mariño et al. (2009). Red points show the self-potential data points, locations of the hot springs after Cruz (2006). Ticsani youngest dome (Ti) and older domes (D1-D3), main faults and collapse zones are shown after Lavallé et al. (2009); Mariño (2002); Gonzales et al. (2007).

Figure 2 :

 2 Figure 2: a) Digital elevation model with elevation contours superimposed on a self-potential map of Ticsani volcano. The reference value of self-potential was chosen in the area of the Toro Bravo lake. Black points show the SP data points, yellow points show the SP data points along a profile Pr1 discussed in the text. Locations of the hot springs after Cruz (2006) are indicated with white diamonds, whose numbers correspond to their order in TableS1. Interestingly, most of the hot springs are at a considerable distance (up to 12 km) from the volcano. b) Cumulative rose diagram showing primary fault and fracture directions present in the area (after[START_REF] Lavallé | Structural control on volcanism at the Ubinas, Huaynaputina, and Ticsani Volcanic Group (UHTVG), southern Peru[END_REF].

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Ground temperature data at Ticsani superimposed on a Spot Google Earth image of the studied area.

Figure 6 :

 6 Figure 6: Numerical model (geometry B shown in Figure 4 b and c) of the fluid circulation below the profile Pr1 for different values of unit1 permeability. White fat solid line shows the position of the water table, control points (black points) are shown like in Figure 4. a) Model B1:k 0 = 4 × 10 -15 m 2 ,like in A1 model, diffusive heat transfer dominates over advection even at shallow depth and the cone is almost fully saturated. b) Model B2: k 0 = 8 × 10 -15 m 2 . The flow is derived towards the lower part of the flank but reaches the surface. Model B2 is the only model that agrees with following observations: the temperature anomaly at more than 10 km from the volcano summit revealed by the hot spring data, the distinct thermal anomaly at the summit confirmed by soil temperature measurements, and the position of the water table given by the self-potential data (see text). a) Model B3: for surface permeability k 0 = 8 × 10 -14 m 2 , the cone is de-saturated, similarly to A2 model. The thermal flow is directed almost vertically towards the summit.

Figure 7 :Figure 8 :

 78 Figure 7: Hot springs temperature versus distance to the summit (full rectangles) with the results of numerical model for different permeability values. Gray diamonds show temperature values in control points and gray lines show smoothing spline interpolation. For k 0 = (8 ± 2) × 10 -15 m 2 the model agrees with observations, while at higher and lower values of permeability the model results deviate significantly from the data.

Figure 9 :

 9 Figure 9: Application of a model of type B2 to Ubinas volcano, k 0 = 6×10 -15 . Results are similar to fig. 8 and show the derivation of hydrothermal system due to gradient of regional topography.

Table 1 :

 1 Physical properties used in numerical models. Parameters in the lower part of the table are identically defined for all units.

	Parameters	unit1	unit2	unit3
	Permeability at surface, k 0 , [m 2 ]	8 × 10 -14 -4 × 10 -15		8 × 10 -14 -4 × 10 -15
	k x	2 k 0	0.1 k 0	0.1 k 0
	k z	k 0	0.1 k 0	2 k 0
	Porosity, φ	0.35	0.1	0.35
	Temperature of the heated body, T heat [ • C]			150
			all units
	Density of rock, ρ s , [kg m -3 ]		2.7 × 10 3
	Density of water, ρ l

Acknowledgments

We thank INGEMMET for financial and technical support of the field work, J. P. Metaxian for his help in the coordination of the work, F. Jouanne for fruitful discussions, S. Hurwitz and J.-F. Lenat for constructive remarks which greatly improved the modeling and the manuscript. We used Latex, JabRef, and InkScape open source software, ASTER GDEM is a product of METI and NASA. A. Revil is supported by DOE (Energy efficiency and Renewable Energy Geothermal Technologies Program) grants. Award #GO18195 and #DE-EE0005513.