Y Gaoua
email: gaoua@laplace.univ-tlse.fr

S Caux
email: caux@laplace.univ-tlse.fr

P Lopez
email: lopez@laas.fr

C Raga
email: mraga@ing.uc3m.es

A Barrado
email: barrado@ing.uc3m.es

A Lázaro
email: alazaro@ing.uc3m.es

Hybrid systems energy management using optimization method based on dynamic sources models

Keywords: Hybrid systems, Energy source model, Linear optimization, Global optimization, Fuel cell modelling, Fuel cell based vehicles I

This paper focus on offline energy management strategy based on dynamic losses computation made on accurate sources models. Hybrid energy systems (mainly Hybrid Electric Vehicle) should now be managed globally to reach the optimal operation minimizing a global cost criterion, such as hydrogen consumption. The presented operational research method is applied to a Fuel Cell based Vehicle (FCV) powertrain, and it has been adapted to the mathematical model depending on the sources characteristics (fuel cell and battery). In this paper losses and efficiency curves are issued from an accurate nonlinear model using experiments for both sources providing data representing the sources behavior to feed the optimization algorithm. Consumption comparisons on actual mission profiles demonstrate the optimal power splitting obtained.

INTRODUCTION

Due to world economic evolution driven by fossil oil market, there is an increasing interest developing hybrid systems. Mainly in transport application, it is now well known full Hybrid Electrical Vehicle should be a solution to fossil oil crisis and NOx emission reduction [START_REF] Chan | The State of the Art of Electric, Hybrid, and Fuel Cell Vehicles[END_REF][START_REF] Chau | An overview of energy sources for electric vehicles[END_REF]. In such HEV almost two energy sources are used, a main source made of a Fuel Cell Stack is used in this research paper associated to a reversible battery source (stores and delivers energy). Autonomy is also increased when consumption is minimized on a given mission. To do so, an accurate behavior of both sources should be determined and an optimization procedure should be implemented to optimally split the power demand between the two sources [START_REF] Brahma | Optimal energy management in series hybrid electric vehicles[END_REF][START_REF] Guemri | Using quasi-Newton method for energy management in electrical multi source systems[END_REF][START_REF] Koot | Energy management strategies for vehicular electric power systems[END_REF].

Depending on the nonlinear behavior of the energy sources a first step should be made in sources modeling and design. The models used for optimization loops should be as accurate as possible and not only based on static data provided by manufacturer in datasheets because the sources are solicited dynamically during an effective mission [START_REF] Adegnon | Experimental evaluation of PEM fuel cell systems efficiency[END_REF].

Considering the sources known and the mission power demand registered with a given step time from the beginning to the end, offline optimization method can be run. In fact the proposed method uses mathematical modeling of the optimization problem and presents an algorithm adapted from Operational Research field. The input data request at instant tk is provided by the models of sources running from tk-1 to tk. The cost function and constraints to respect are included in the mathematical problem and different objectives are proposed: Pb1: Minimize the total energy losses in the system (global efficiency) while respecting all the different constraints of system operating points.

Pb2: Minimize the total power provided by the FC (main source consumption) while respecting all the different constraints of the system operating points.

The paper is organized as follows: Part II presents the nonlinear electrical equivalent models of the battery and the fuel cell source. R, C elements present different values, due to various operating point used during the mission (power to be delivered), identification is based on various experiments allowing thus computing accurate losses for optimization purposes. Part III presents the mathematical formulation of the optimization problems (Pb1 and pb2) using previous losses computed data. Part III presents results obtained with actual data of a Fuel Cell based Vehicle. Conclusions and perspectives are given at the end.

II. PROPULSION SYSTEM AND ENERGY SOURCES MODELLING

A.

Powertrain chain

The powertrain analyzed, Fig. 1, shows the fuel cell as main energy source and a battery as secondary energy source, connected in parallel to the dc voltage bus. All the power and energy that the fuel cell is not able to provide will be delivered by the battery. In braking phase, energy recovery is stored in battery (no fuel cell reversibility). The nominal power of the fuel cell is 1.2kW, and the secondary source is a Lithium-ion battery characterized by its nominal voltage Vnom =72V and its capacity Qmax=4.549 Ah. The optimal sizing of the powertrain energy systems is thoroughly described in [START_REF] Raga | Optimal sizing of propulsion systems applied to fuel cell based vehicles[END_REF].

B.

Powertrain losses chain In order to deal with the mathematical formulation, the following block diagram of Fig. 2 is considered. The moto-propulsion group is resumed in power demand on the DC-bus (positive power for traction, negative in braking phase, zero at stop).

C.

Battery modeling Several works have been conducted in order to deduce the battery voltage behavior according to its operation mode: charge or discharge mode. The model used is the provided one by MATLAB ™, and the voltage equation is defined by (1):

(1)

Where is battery constant voltage (V), K is the polarization constant (V/Ah), A is magnitude in exponential zone (V), B is the inverse time constant (Ah -1 ), ∫ the instantaneous battery charge (Ah). Doing so, an equivalent nonlinear electrical model presented Fig. 3 is derived. The battery voltage depends on its capacity and the current passes through its internal resistance Rv. , K, A, B are determined (identification are presented in [START_REF] Tremblay | A Generic Battery Model for the Dynamic Simulation of Hybrid Electric Vehicles[END_REF]), using PowerSim Matlab™ toolbox, not only with data sheets, but also with several operating points coming from manufacturers and measurements and analyzed with Fig. 4.

In order to use the battery voltage equation in an optimization mathematical model, polynomial approximation is necessary to transform the nonlinear equation into polynomial equation reducing the computational complexity of the energy management algorithm and not losing accuracy. , can be also written using battery SOC after polynomial approximation (degree n = 10).

∑

(2)

D. Fuel Cell modelling

The main source in the propulsion system is a fuel cell, which produces energy via the chemical reaction of hydrogen and oxygen. The advantages of using fuel cell in electric vehicles derive from its high efficiency that exceeds 40 %, the generation of electrical energy with little fuel consumption and no NOx emission, unlike the performance of Internal Combustion Engine which varies between 25 % and 30 % emitting pollutant.

There are two possible voltage transient responses when a load current step is applied to a fuel cell: overdamped and overshooted [START_REF] Adzakpa | PEM Fuel Cells Modeling and Analysis Through Current and Voltage Transient Behaviors[END_REF], and several fuel cell models reproducing both kind of behaviors can be found in the literature [START_REF] Boscaino | Fuel cell modelling for power supply systems design[END_REF]- [START_REF] Raga | Black-box model and identification methodology for PEM fuel cell with overdamped transient response[END_REF].

A commercial Nexa Ballard 1.2kW PEM fuel cell, with overshooted transient response has been modeled, by means of the identification methodology described in [START_REF] Raga | Black-Box Model, Identification Technique and Frequency Analysis for PEM Fuel Cell with Overshooted Transient Response[END_REF]. This procedure is based on manufacturer"s datasheets and simple measurements in the time domain.

The PEM fuel cell model is shown in Fig. 6, where E oc is the open circuit voltage, resistors R s and R v (i fc ) model the static I-V curve, and the capacitor C p , the inductor L and the resistor R L (i fc ) model the dynamic behavior of the fuel cell. The diode assures the unidirectional fuel cell operation. The resistors R s and R v (i fc ) values are calculated from the quasi-static I-V curve, which can be provided by the manufacturers or by means of simple measurements on the fuel cell. The value of the resistor is constant Rs=120m, and is calculated from the fuel cell output voltage considering only the ohmic region. Due to this fact, the resistor R v (i fc ) gets variable values depending on the fuel cell current i fc . Table I shows the collected data presenting the variation. The fuel cell quasi-static model, parameterized following [START_REF] Raga | Black-Box Model, Identification Technique and Frequency Analysis for PEM Fuel Cell with Overshooted Transient Response[END_REF], is shown in Fig. 7, where E oc = 43V. And finally, the behavior of the fuel cell model is shown together the experimental response in Fig. 8. R L (i fc ) In order to reproduce the quasi-static FC behavior and run the mathematical simulations, the FC voltage equation given in (3) can be written.

R s =120m i fc + - v fc C p =330mF E oc =43V + _ R v (i fc ) L=5.54H R v (i fc ) R s i fc + - v fc E oc + _
(3) With the aim of using this model in the optimization step, the FC voltage behavior is computed with polynomial approximation (degree n = 15) to identify the R v equation (4). Allowing now, running simulations with the electrical equivalent model shown in Fig 7.

∑

To be noticed: the FC is connected to the DC bus through power electronic converter device. The battery is directly linked to this DC bus. Thus the power demand (indexed busload) is the sum of battery power and FC power on bus side (indexed bus to indicate exit of FC stack after the converter). Knowing that the efficiency of the boost converter is  generally [0.95; 1] (5) is used.

P busload =P B +P bus (power splitting) will be managed when the vehicle is in traction or when the vehicle brakes.

(5)

The losses and exact computation of all electrical values allows feeding the hereafter presented optimization algorithm with effective input data. The objective of the optimization is now to split the power demand dealing with a compromise using fuel cell or battery to minimize a criterion. These models (battery and fuel cell), can also be used in PSIM or Matlab/Simulink, or HIL environment to dynamically compute the current, the voltage and thus losses in both elements, when references are given to deliver the power demand (every 1s in our case) running from the beginning to the end of the mission.

E.

Driving profile The energy and power required by the load (dc motor and vehicle), depends the applied driving cycle [START_REF] Sadoun | Influence of the drive cycles on the sizing of hybrid storage system battery-supercapacitor supplying an electric vehicle[END_REF]- [START_REF] Yufang | Impact of driving cycles and all-electric range on plug-in hybrid vehicle component size and cost[END_REF]. This analysis is focused on a light vehicle of 800kg, propelled by a motor drive of nominal power 4.7kW which is able to provide a peak power of 18kW. An European urban driving cycle is applied, ECE-15, which maximum speed of 50km/h and the total time is 181s. Fig. 9 shows the load power requirement that derive from the driving cycle and the vehicle features. In order to evaluate the fuel consumption during a working day of 8h, the aforementioned driving cycle will be applied 161 times. Load power (kW) requirement.

III. OPTIMIZATION FORMULATION

A. Mathematical Formulation

The goal of this study is managing the energy distribution of the system in order to reduce the total power provided by the FC, and therefore the consumed hydrogen. To do this, a mathematical modeling is developed taking into account the satisfaction of the power demand on DC-bus while respecting the different constraints of system functioning and source design. The non-linear modeling developed to minimize the hydrogen consumption contains the following decision variables: 7) is used to minimize the energy quantity supplied by the FC (Pb2 = H 2 consumption is minimized) during the overall mission. ∑

P fc (
And ( 8) used for the energy losses minimization of the global system (Pb1 global efficiency).

∑ ( )

The constraints of the model are:

-Satisfaction of Power demand when the vehicle is in traction (P>0) and recovering energy from the load and the FC if necessary, when the vehicle brakes (P<0):

-The voltage of the FC is calculated using the fixed voltage, static and variable resistors (3).

-The value of the variable resistor is calculated using the FC current and polynomial approximation of degree 15 (4) in order to reduce errors.

-The FC can supply a maximum power of P max and minimum power of P min to maintain the proper functioning of its auxiliary:

-The power supplied/recovered by the battery:

-Evolution of the state of charge of the battery ( )

-Storage capacity of the battery limited to:

-The battery voltage expressed in [START_REF] Chau | An overview of energy sources for electric vehicles[END_REF].

By summarizing all the equations, the following global model (Pb2) is obtained: x. The functions f(x) and g(x) can be nonlinear and nonconvex, but should be twice continuously differentiable that is insured by the proposed polynomial formulation obtained in previous parts. The functions in this problem are convex so, the optimality of the solution found is granted.

∑ ∑ s.t. ( )
s.t.

IV. SIMULATION RESULTS

Once the power management optimization is applied, the set of voltage and current signals within the propulsion system is obtained, therefore the current reference to be applied in the converter control loop, i ref (Fig. 12) is fixed. Current control loop implemented in the powertrain First of all, the current reference signal obtained by means of both approaches (Pb1 or Pb2) are compared (see Fig. 13), and it can be observed that the obtained result from both objective functions is exactly the same mainly due to the fuel cell losses (due to FC current) which is of paramount importance in this hybrid system. Therefore, henceforth only the fuel cell power minimization data will be used. In order to verify that the developed mathematical model used in the power management optimization is well done, the current reference signal obtained in the power management algorithm is applied to the boost converter control loop. Then, the main current and voltage waveforms evolutions from the mathematical optimization and from the electric simulation are compared, Fig. 14. It can be observed that the signals evolution of both the mathematical modeling and the electrical simulation is the same. Then, it can be conclude that the mathematical algorithm has been properly developed.

Once validated the mathematical formulation, next step consist on comparing the fuel cell hydrogen consumption when the power management (PM) algorithm is applied or not. To do so, the powertrain is simulated in PSIM, and the delivered power by FC is compared when the current reference signal is directly the load demanded current i busload , or when the current reference signal is the obtained current by means of the optimization algorithm. Fig. 15 shows the fuel cell delivered power with and without power management techniques. The initial battery state of charge is 75%, and by means of the power management optimization algorithm, the final state of charge is around 10%. When the power management is not applied, the battery is discharged along the working day, and when its state of charge is lower than 10% automatically a recharge procedure is applied, until the state of charge reaches the 75%.

Table II shows the FC energy delivered, with and without power management techniques. It is concluded that within a working day, the difference in the hydrogen consumption is a 24% gain. V.

CONCLUSIONS

This paper presents a work made in collaboration between two research teams. One based in Spain working on energy sources modeling and the other one based in France focusing their research on energy management optimization method. The discrete time offline optimization is linked in this paper to an accurate and experimental-based sampled non-linear behavior of energy sources. Results obtained on actual data mission profiles allow concluding the method developed is correct, accurate and the optimal results is granted due to the effectiveness of the losses "dynamically" computed.

The objective functions used to minimize both the total energy losses in the system and the total power delivered by the fuel cell are described. The use of these functions require adequate fuel cell and battery models with the aim of computing both sources behavior, assuring the same operation than the actual FCV under consideration. Then, the optimal power splitting between both sources is achieved, and the global minimum losses in the overall system mission are reached. An important gain in terms of hydrogen energy consumed on the overall mission is obtained increasing the FCV autonomy.

On-board optimization linked to traffic information can be easy to implement and mission profile updated in a given timewindow can replace input data ECE-15. Photovoltaic Panels and/or Supercaps models can also be used in this energy management algorithm. On-line optimization algorithm using these sources model and this mathematical form, also can be now prospected, to be more robust without being too heavy in computation time, allowing to optimally react on disturbances present on an actual driving cycle or uncertainty on sources behavior (default, ageing etc.). I.

Fig. 1 .

 1 Fig. 1. Propulsion system block diagram

Fig. 2 .

 2 Fig. 2. Electrical elements and their losses in the powertrain losses chain According to the structure presented Fig. 1, a battery model and a fuel cell model are established. The modeling is oriented to well compute current and voltage behavior allowing computing instantaneous losses (at tk) in equivalent resistance whatever the battery or fuel cell operating points.

Fig. 3 .

 3 Fig. 3. Electrical model of battery source

Fig. 4 .

 4 Fig. 4. Simulation based parameters identification Simulations, using the static battery model from MATLAB, were performed by fixing the discharge (resp. charge) current to 1A (resp. -1A), with the aim of measuring both the battery open circuit voltage and the State Of Charge (SOC). After that, the battery open circuit as a function of the SOC is obtained (K=0.099, A=8.11, B=13.42), as Fig. 5 depicts.

Fig. 5 .

 5 Fig. 5. Battery behavior approximated curve

Fig. 6 .

 6 Fig. 6.

Fig. 7 .

 7 Fig. 7. Static model of the PEM fuel cell.

Fig. 8 .

 8 Fig. 8.Comparison between measured fuel cell I-V curve and simulated fuel cell model I-V curve.

Fig. 9 .

 9 Fig. 9.Load power (kW) requirement.

∑Fig 10

 10 Fig 10 presents a blockset of the problem computation, presenting first the data computation block, feeding the cost evaluation of the objective function, respecting (or not) constraints allowing step by step to search the optimal splitting at each instant tk minimizing the overall objective function.

Fig. 10 .

 10 Fig. 10. Global simulation environment The resolution of the problem is made off-line using a non linear algorithm IPOpt (Interior Point Optimizer) combined with AMPL modeler (Modeling Language for Mathematical Programming) to write the mathematical model in appropriate format. The IPOpt algorithm is an open source solver of COIN-OR platform solving a nonlinear problem (NLP) see [15-18], as shows Fig. 11.

Fig. 11 .

 11 Fig. 11. AMPL example where f(x): R^n --> R is the objective function, and g(x): R^n --> R^m are the constraint functions. The vectors g_L and g_U denote the lower and upper bounds on the constraints, and the vectors x_L and x_U are the bounds on the variablesx. The functions f(x) and g(x) can be nonlinear and nonconvex, but should be twice continuously differentiable that is insured by the proposed polynomial formulation obtained in previous parts. The functions in this problem are convex so, the optimality of the solution found is granted.

Fig. 12 .

 12 Fig. 12.Current control loop implemented in the powertrain

Fig. 13 .

 13 Fig. 13. Current reference signal obtained from fuel cell power minimization and from energy losses minimization, along one ECE15 driving cycle

Fig. 14 .

 14 Fig. 14. Mathematical optimization and electric simulation of the powertrain behavior comparison along one ECE15 driving cycle

Fig. 15 .

 15 Fig. 15. Fuel cell delivered power with and without power management technique.

  *Ibat)/3600000 INT(Vbat_(PSIM_PM_6336)*Ibat_(PSIM_PM_6336))/3600000 P FC (W)

TABLE I .

 I STATIC FUEL CELL MODEL PARAMETERS

	i fc (A)	R v (i fc ) (	i fc (A)	R v (i fc ) (
	1	1.780	25	0.312
	2	1.080	30	0.279
	3	1.013	35	0.251
	4	0.905	40	0.229
	5	0.820	42	0.225
	8	0.630	43	0.222
	10	0.550	44	0.221
	15	0.427	45	0.227
	20	0.358

  State of charge of the battery at time t. The objective function (

	I B (t): Current supplied/recovered by the battery at time t,
	V B (t): Battery voltage at time t,
	SOC(t):
	i_motor*Vbus
	2K
	1.5K
	1K
	0.5K
	0K
	-0.5K
	-1K
	50	100	150
		Time (s)
	t): Power supplied by the FC at time t,
	V fc (t): FC voltage at time t,
	I fc (t): Current supplied by the FC at time t,
	Ibus(t): Output current of the boost converter at time t,

Table II .

 II Results applying or not power management along a working day

		Applying PM	Without PM
	FC delivered energy (Wh)	967	1200


Alex Fidahoussen

Pierre Calmon

Marc Lambert

Ultrasonic phased array reconstruction methods for the localization and the characterization of defects in inspection of complex geometry configuration

Phased arrays are now widely used in ultrasonic non destructive testing, as they provide improved flexibility and adaptability to complex configurations with comparison to conventional probes. One other advantage of the phased arrays techniques is that they can give a large collection of data. To fully exploit this data, in order to localize and to characterize the defect, imaging and reconstruction methods based on a direct modelling of the inspection are required. The work presented here aims at developing such reconstruction methods which use the direct models implemented in the CIVA platform. These models allow to consider complex inspection configurations and in particular geometries with irregular surface. The proposed algorithms are derived from the so called synthetic focusing approach which is classically based on time-of-flight identification, but they also exploit CIVA modelling of the signal echo from a point-like scatterer distributed inside a region of interest used for the reconstruction area. The adopted approach is generic and can be used with different operating modes of the array: application of electronic delays at emission or/and reception, sectorial scanning (several delay laws are applied to sweep the beam at different refraction angles), independent storage of the signals received by the elements… In this communication we present the developed algorithms and discuss their performances with application examples on both simulated and experimental data.


Context of the study Theoretical formulation Validation

Context and configuration of the study

A(r) = A(r , z) u θ + A ec (r , z) = R(r )W (z)-lead to ∂ 2 ∂r 2 R(r ) + 1 r ∂ ∂r R(r ) + a 2 - 1 r 2 R(r ) = 0 (1) ∂ 2 ∂z 2 W (z) = a 2 + jωµσ(z) W (z) (2)
Following [2] the general solution given by

W (z) = CF 1 (f (z)) + BF 2 (f (z)) (3)
where F 1 and F 2 typical mathematical functions related to σ(z) Expression of

A 1c (r , z) = +∞ 0 µNII (r 1 , r 2 )e -az J 1 (ar ) e -ah1 -e -ah2
2a 3 (h 2 -h 1 )(r 2 -r 1 ) da

+ +∞ 0 C 1 e -az J 1 (ar )da, with I (r 1 , r 2 ) = ar2 ar1 xJ 1 (x)dx ∀z ∈ [0;+∞[ A 2 (r , z) = +∞ 0 [C 2 F 1 (f (z)) + B 2 F 2 (f (z))] J 1 (ar )da ∀z ∈ [-r ;0] A 3 (r , z) = +∞ 0 B 3 F 3 (g (z))J 1 (ar )da ∀z ∈] -∞; -r ] Formulation for a two-tanh profile A 1c (r , z) known, then Z deduced as Z = K +∞ 0 I (r 1 , r 2 ) 2 a 6 2 e -a(h2-h1) -1 + a(h 2 -h 1 ) + e -ah2 -e -ah1 2 φ(a) da (4) with φ(a) = C 1 K ;K = µNII (r 1 , r 2 ) e -ah1 -e -ah2 2a 3 (h 2 -h 1 )(r 2 -r 1 ) ;K = -j ωπµN 2 (h 2 -h 1 ) 2 (r 2 -r 1 ) 2 .
Continuity conditions of the quantities and/or their derivatives with respect z and/or their cancellation at ±∞ express φ(a)as

φ(a) = (aM -O)RT -(aL -N)ST + [a(LQ -MP) -NQ + OP] U (aM + O)RT -(aL + N)ST + [a(LQ -MP) + NQ -OP] U (5)
where

L = F 1 (f (z = 0)); M = F 2 (f (z = 0)); N = F ′ 1 (f (z)) z=0 ; O = F ′ 2 (f (z)) z=0 ; P = F 1 (f (z = -r )); Q = F 2 (f (z = -r )); (6) R = F ′ 1 (f (z)) z=-r ; S = F ′ 2 (f (z)) z=-r ; T = F 3 (g (z = -r )); U = F ′ 3 (g (z)) z=-r
′ means derivative with respect to z

Formulation for a two-tanh profile

Conductivity profile given by

σ(z) =          σ 12 + σ 1 -σ 12 2 1 + tanh z + c 1 2v 1 ∀z ∈ [-r , 0] σ 2 + σ 12 -σ 2 2 1 + tanh z + c 2 2v 2 ∀z ∈] -∞, -r ] (7)
Particular functions F 1 , F 2 and F 3 are

F 1 (y 2 (z)) = y µ 2 (z) [1 -y 2 (z)] ν F (µ + ν, µ + ν + 1,2µ + 1;y 2 (z)) (8)
F 2 (y 2 (z)) = y -µ 2 (z)[1 -y 2 (z)] ν F (ν -µ + 1,ν -µ, -2µ + 1;y 2 (z)) (9) F 3 (y 3 (z)) = y λ 3 (z) [1 -y 3 (z)] τ F (λ + τ, λ + τ + 1,2λ + 1;y 3 (z)) (10) with y 2 (z) = 1 + e - z+c 1 v 1 -1 , y 3 (z) = 1 + e - z+c 2 v 2 -1 µ = v 1 a 2 + jωµ 0 σ 12 , ν = v 1 a 2 + jωµ 0 σ 1 λ = v 2 a 2 + jωµ 0 σ 2 , τ = v 2 a 2 + jωµ 0 σ 12 (11)
F (α,β,γ;x) is the hypergeometric function Description of the configuration Comparison with the results given in Description of the configuration

r 1 1.3 mm r 2 3.3 mm h 1 0.5 mm h 2 7.8 mm N turn 580 σ 1 1.883 10 7 S m -1 σ 12 3.766 10 7 S m -1 c 1 0.3 mm v 1 0.1857 mm
[1] Real part of Z , (b = NR 2 ) Imaginary part of Z , (b = NR 2 ) Frequency from [1] N = 10 N = 20 from [1] N = 10 N = 20
r 1 = 2.0 mm σ 1 = 6 10 5 S m -1 r 2 = 4.0 mm σ 12 = 9 10 5 S m -1 h 1 = 0.5 mm σ 2 = 8 10 5 S m -1 h 2 = 7.3 mm N turn=200 c 1 = 0.2 mm v 1 = 0.03 mm c 2 = 0.8 mm v 2 = 0.1 mm

  Eddy-Current measurements over combustion turbine blade coatings affected by depletion of aluminium; Model taking inward and outward depletion of aluminum inside the coating into account; Conductivity profile follows a two-hyperbolic-tangent law; Analytical formulation of the variation of impedance obtained combining the approaches found in[1, 2, 3]

Frédéric Nougier, Marc Lambert, Riadh Zorgati

Frédéric Nougier, Marc Lambert, Riadh Zorgati 8/ 8


Nicolas Coppo
email: n.coppo@brgm.fr

Jean-Michel Baltassat
email: jm.baltassat@brgm.fr

Jean-François Girard
email: jf.girard@brgm.fr

Pierre Wawrzyniak
email: p.wawrzyniak@brgm.fr

Sophie Hautot
email: sophie.hautot@imagir.eu

Pascal Tarits
email: tarits@univ-brest.fr

Thomas Jacob
email: t.jacob@brgm.fr

Guillaume Martelet
email: g.martelet@brgm.fr

Francis Mathieu
email: f.mathieu@brgm.fr

Alain Gadalia
email: a.gadalia@brgm.fr

Vincent Bouchot

Hervé Traineau
email: h.traineau@cfg.brgm.fr

3-D Magnetotelluric Investigations for geothermal exploration in Martinique (Lesser Antilles). Characteristic Deep Resistivity Structures, and Shallow Resistivity Distribution Matching Heliborne TEM Results

Keywords: Magnetotelluric, 3-D inversion, Montagne Pelée, Anses d'Arlet, Martinique, TEM, Geothermal exploration

  L'archive ouverte pluridisciplinaire

INTRODUCTION

The geothermal potential of Martinique Island (Lesser Antilles, France) has been extensively investigated since more than 50 years. The first exploration started in the 60' in the Plaine du Lamentin (south of Fort-de-France) and has been extended in this area in the 80's and orientated towards the southeastern flank of the Montagne Pelée. Then, the southwestern flank of the Montagne Pelée and the Anses d'Arlet) were identified of potential interest and explored in the 2000's. Finally, three exploration wells were drilled in the Plaine du Lamentin in 2001 and confirmed the presence of a medium temperature geothermal resource. This project had for main objectives to first compile existing knowledge about the geothermal potential of Martinica Island, but also to proceed to adequate additional field exploration with the aim to complete and strengthen the existing geoscientific database and finally provide an up-to-date interpretation for the four selected areas (Anses d'Arlet, Montagne Pelée, Pitons du Carbet and Lamentin) with proposition of the most favorable zones for exploration wells. New geological (exploration of mineralogical clues related to magmatic and hydrothermal activity and structural context), hydrogeological (geochemistry and recharge assessment), geochemical (soil and water gas analysis) and geophysical data were acquired on these four sectors.

Beyond crucial geological, geochemical and hydrogeological data, geothermal exploration requires additional images of the physical properties of subsurface to understand and interpret properly the geothermal reservoir and its potentialities. This information can only be provided by adapted geophysical measurements above areas of interest. Because the amount of data collected during this two-years project is too important to be rigorously presented and described in a single communication, this paper mainly focuses on the geophysical results (new magnetotelluric and gravimetrical data, collected over the Anses d'Arlet and the Montagne Pelée province.

GEOLOGICAL AND GEOTHERMAL OVERVIEW

The Lesser Antilles island arc, which results from the westward subduction of the Atlantic Plate under the Carribbean Plate, has been active since the Early Eocene [START_REF] Nagle | K-Ar geochronology of the Limestone Caribbees and Martinique, Lesser Antilles, West Indies[END_REF]. Arc evolution has been described as occurring in three stages, representing the old, intermediate and recent arcs, corresponding to Late Oligocene-Early Miocene, Middle Miocene and Late Miocene to recent ages, respectively [START_REF] Briden | K-Ar geochronology and palaeomagnetism of volcanic rocks in the Lesser Antilles island arc[END_REF][START_REF] Westercamp | Diversité, contrôle structural et origines du volcanisme récent dans l'arc insulaire des Petites Antilles[END_REF], and having erupted from east to west. Due to the geological evolution of this subduction zone, Martinique is the northernmost island where the products of the three stages crop out from east to west without superimposition.

At the end of the Miocene, the locus of volcanism jumped to the north of the island producing the Morne Jacob shield volcano, whose oldest rocks are dated at 5.14+-0.07 Ma [START_REF] Germa | The eruptive history of Morne Jacob volcano (Martinique Island, French West Indies): geochronology, geomorphology and geochemistry of the earliest volcanism in the recent Lesser Antilles arc[END_REF], making the onset of the recent arc activity in Martinique Island. From 2.3 Ma, Trois Ilets volcanism commenced in the southwest, coeval with continued activity at Morne Jacob volcano (Germa et al., 2011a). Trois Ilets Volcanism, Carbet Complex and Conil-Pelée Complex were then active contemporaneously until 345 ka (Germa et al., 2011a). Most recently, volcanic activity has been restricted to the northern Conil-Pelée Complex, whose last eruption occurred in [1929][1930][1931][1932]. Lavas from the recent arc are tholeitic (Morne Jacob, early stage) to calc-alkaline, with basalticandesite to dacite compositions [START_REF] Labanieh | Géochimie de l'île de la Martinique aux Petites Antilles[END_REF][START_REF] Germa | The eruptive history of Morne Jacob volcano (Martinique Island, French West Indies): geochronology, geomorphology and geochemistry of the earliest volcanism in the recent Lesser Antilles arc[END_REF]Germa et al., , 2011a)). Each complex has its own geochemical signature, for major and trace elements, and for isotopic composition [START_REF] Smith | Mt Pelee, Martinique: a study of an active island-arc volcano[END_REF][START_REF] Labanieh | Géochimie de l'île de la Martinique aux Petites Antilles[END_REF]; [START_REF] Germa | The eruptive history of Morne Jacob volcano (Martinique Island, French West Indies): geochronology, geomorphology and geochemistry of the earliest volcanism in the recent Lesser Antilles arc[END_REF]Germa et al., , 2011a)).

Anses d'Arlet geothermal field (see Gadalia et al., 2015, WGC)

Recent geological observations and analysis collected during this project have improved and specified the geological context of the Anses d'Arlet prospect (Traineau et al., 2013). In particular, this geological work provides a consistent framework to the magmatic and hydrothermal activity of the southwest of the Trois-Ilets Peninsula. Detailed mapping of the fossil alteration zones connected to their tectonic controls allowed to better constrain the geothermal reservoir context and its evolution up to now. A three steps continuous evolution is proposed in agreement with previous works by Germa et al. (2011) and Westercamp et al. (1990).

The initial stage is characterized by the development of a hydro-magmatic system thought to be related with the polyphased activity of the Roches Genty eruptive center (not dated). Many volcanic edifices grew up from north-west to south-est (Pointe Burgos-Le Diamant volcanic axis date between 1.5 and 0.35 Ma). Geological mixing between deep basaltic lavas and more superficial dacitic lavas with quartz suggests that magmatic differentiation occurred at relatively shallow depth (6 km from mineralogy). It is not obvious whether a single or multiple magmatic chambers fed all of these the volcanic edifices. The local tectonic setting (NW-SE and NE-SW structures) of the Trois-Ilets peninsula may have favored successive ascents of a deep basaltic magma during at least 1.1 Ma. These magma chamber(s) and associated intrusions would have been the heat source of an old acidic hydrothermalmagmatic system (early hydrothermalism) related to Roches Genty eruptive center.

The intermediate stage is characterized by the decline of the magmatic component and the evolution of the hydrothermal-magmatic acidic system toward a large high temperature neutral geothermal system. A clay cap develops together with a surface fumarolic activity along a corridor connecting the Anses d'Arlet and Petite Anse, reflecting the evidence of boiling related to thishigh temperature fossil system. Structural analysis demonstrates that faults reworked many time in stress tension creating favorable conditions for fluid circulation. This hydrothermalism can be described by three alteration products, spatially distinct: evidence of a kaolinite-alunite alteration (relief), the presence of sinter (silica deposits) and the existence of clay mineral with smectite (dominant) and trace of interstratifed illite/smectite. During the last 0.5 Ma, the Morne Jacqueline experienced a lateral collapse that makes this fossil caprock outcropping together with intrusions.

The final stage represents the current geothermal activity expressed by the Eaux-Ferrées thermo-mineral springs in the Petite Anse bay where degassing and small travertine deposits can be observed. The small flow (0.03 l/s) of the springs associated with the fossil character, fumarolic and of high temperature, characteristic of the alunite-kaolinite alteration suggest a reduction both in extension and intensity of the surface manifestations. Deeper this could indicate a contraction (or a migration) of the HT geothermal system southwards.

Montagne Pelée geothermal field

Recent geological observations carried out during this project and reported by Traineau et al. (2013) completed the geological knowledge of the Montagne Pelée. Particularly, the pelean volcanic system, its formations, its magmatic heat source, its volcanotectonic structures and the context of the surface hydrothermal manifestations disseminated on the southwestern flank of the stratovolcano. The Montagne Pelée is one of the most active volcano of the Antilles with 28 magmatic eruptions during the last 16 ka (Germa et al, 2011). Edifice growth and its plumbing system started about 127 ka ago [START_REF] Germa | The eruptive history of Morne Jacob volcano (Martinique Island, French West Indies): geochronology, geomorphology and geochemistry of the earliest volcanism in the recent Lesser Antilles arc[END_REF][START_REF] Boudon | Role of large flank-collapse events on magma evolution of volcanoes. Insights from the Lesser Antilles Arc[END_REF]. Lavas geochemistry shows that the magmatic system nowadays feeding the Montagne Pelée is the same that fed Mont Conil in the past, both being located on the same NW-SE volcanic axis, the major tectonic direction of Martinica.

The main surface hydrothermal manifestations testifying a deeper geothermal activity are mainly located on the southwestern flank of the volcano. They are sulfuric springs of the uphill Rivière Claire (nowadays disappeared), fumaroles of the Etang Sec (also disappeared), hot springs from rivers Chaude Mitan and Picodo (still existing). The geochemical components of the two latters indicate the existence of another geothermal reservoir that could be separated from the main one. Hot spings and well located downhill on the coast are interpreted as a lateral extension of the river Chaude spring system. Other "cold" springs have a geochemistry indicating a deep component but very diluted.

No recent tectonic structure, except the three lateral collapses scars proposed by Le Friand et al. (2003) were identified on the southwestern flank. The first flank collapse occurred about 100 ka ago with a volume estimated to 25 km 3 , the second (13 km 3 ), about 32 ka ago. Both of them considerably modeled the topography of the southwestern flank of the Montagne Pelée. We consider they had major implications on the geologic volcanic deposits pile accumulated around the volcano, causing deeper altered part of the edifice to reach the surface with a probable thinning of the caprock. Precise triggering factors of such catastrophic event are still a matter of debate but can be related to volcanic activity, caldera collapse, strong rainfall, climatic change [START_REF] Coppo | A deep scar in the flank of Tenerife (Canary Islands): Geophysical contribution to tsunami hazard assessment[END_REF][START_REF] Quidelleur | Causal link between Quaternary paleoclimatic changes and volcanic islands evolution[END_REF], seismic activity, fault activity. A third smaller (2 km 3 ) flank collapse occurred about 9 ky ago and had a reduced impact of the stratovolcano structure.

To the summit of the Montagne Pelée, two calderas and eruptive conduits have been identified by [START_REF] Westercamp | Carte géologique de la Montagne Pelée à 1/20 000 (département de la Martinique) et notice explicative simplifiée[END_REF]. A larger old one has been under discussion since the evidence of lateral collapse.

MAGNETOTELLURIC METHOD

Different rocks, sediments, geological structures and fluids have a wide range of different electrical conductivities. Therefore, measuring electrical resistivity allows different materials and structures and fluids to be distinguished -sometimes identified -from one another and can improve knowledge about type of fluid content and its distribution, clay content, tectonic and geological structures.

The magnetotelluric method is a passive geophysical technique for imaging the electrical conductivity and structure of the Earth from the near surface down to several tens of kilometers [START_REF] Vozoff | The magnetotelluric method, Electromagnetic Methods in Geophysics[END_REF]. Its principle consists in simultaneous measurements and recordings of natural (internal) variations of electric and magnetic fields at the earth surface whose strength and geometry are governed by the sub-surface resistivity distribution. The primary external natural signal comes mainly from both the solar wind that induces pulsations in the magnetosphere and the lightning activity which both cause natural variations in the earth's magnetic field, inducing electric currents flowing under the Earth'surface. These naturally varying electromagnetic fields are measured over a wide range of magnetotelluric frequencies (from 10 kHz to 1000 s in our case). Low frequencies (< 1 Hz) are generated by the solar activity and highest (> 1 Hz) by worldwide thunderstorm activity. Combined, these natural phenomena create MT source signals over the entire frequency spectrum. Therefore, data quality is strongly dependent on external source activity and can be severely reduced -especially in the low frequency rangeduring phases of poor solar activity.

Magnetotelluric data comprises time series of the horizontal electric components (Ex and Ey measured in orthogonal directions) with their two associated magnetic components (Hx and Hy). When possible the vertical magnetic component (Hz) is also recorded. The ratio of the electric field to the magnetic field (E/H) at selected frequency, known as the "impendance tensor" embodies the full information about sub-surface conductivity structure. The ratio E/H is usually displayed as both apparent resistivity and phase as a function of frequency. Because of the "skin effect" phenomenon (explaining the attenuation of electromagnetic fields with depth), the depth of penetration of EM waves increases with lower frequency and higher resistivity. At most sites, collocated transient electromagnetic data were also acquired to correct the magnetotelluric data for the amplitude of galvanic distortion [START_REF] Sternberg | Correction for the static shift in magnetotellurics using transient electromagnetic soundings[END_REF] caused by small-scale, shallow conductors.

A variety of techniques are available to assess the dimensionality of the causative structures from the impedance tensor, and, in the case of a two-dimensional (2D) Earth, determine the geoelectrical strike direction [START_REF] Simpson | Practical Magnetotellurics[END_REF][START_REF] Caldwell | The magnetotelluric phase tensor[END_REF]. In the present work, 3-D inversion has been directly applied to integrate tensors multi-dimensionality and sea-effect (see below).

Data acquisition

MT data have been acquired using Metronix MT stations (four ADU06 and three ADU07), MFS06e, MFS07 and MFS07e induction coils (for the magnetic field) and unpolarizable electrodes for electric field (Petiau, 2000). All MT stations were GPS synchronized. A remote reference has been permanently deployed in the Pitons du Carbet massif, in the center part of Martinique Island during the survey. Every day three new MT sites were occupied according to the survey design and objectives, and to local constraints (accessibility, authorization, proximity of noise sources). In order to increase data quality in the low frequency band (< 1Hz), our acquisition program has been extended to two days (three during the week-end). A total of 36.5 hours of low frequencies (sampling frequency = 256 Hz), 2 hours of intermediate frequency (sampling frequency = 2048 Hz) and 9 minutes of high frequencies (sampling frequency = 65536 Hz) were recorded. The remote reference was programmed to run on a daily basis in order to be precisely synchronized with the 6 MT stations deployed simultaneously on the field. Magnetic sensors were systematically buried at about 20-30 cm depth and electric and magnetic cables fixed or buried in the ground. A bentonite mud was prepared to insure good and long contact resistance of unpolarizable electrodes. Data quality revealed to be during better November than in December 2012 because of a lower solar activity during the last period as shown by comparison with geomagnetic observatories.

Processing

Before processing, data format was first homogenized and data arranged in a database dedicated for processing. Due to the large amount of data collected, a daily quality check was performed to insure both acquisition parameter accuracy, correct synchronization with the remote reference and data quality. We then applied a robust processing with remote reference for frequencies above 0.1 Hz based on a robust estimation of the impedance tensor errors described in Wawrzyniak et al. (2012). This statistical approach, based on BIRRP code (Chave and Thompson, 2004), allows to improve impedance tensor estimates providing that time series are long enough. Consecutively, low frequencies were processed using BIRRP. MT tensors and interpretational quantities were then plotted as function of frequency for analysis.

3-D MT inversion

Because of the complex volcanic island context, merging proximity of the sea with complex shoreline morphology, combined conical structure (with radial and circular geological objects) of a volcanic edifice and the linear main constitutive tectonic structures, 3-D inversion has been run with sea effect taken into account, to avoid potential distortion of the MT tensors. This task has been performed using the minim3D code [START_REF] Hautot | The deep structure of the Baringo Rift basin (central Kenya) from 3-D magneto-telluric imaging: Implications for rift evolution[END_REF][START_REF] Hautot | 3-D magnetotelluric inversion and model validation with gravity data for the investigation of large igneous provinces[END_REF]. A first solution to solve this effect is to consider the shoreline and the bathymetry around Martinique island as model parameters. However, a compromise between model resolution and time computation had to be chosen for the inversion. Furthermore, taking a coarse bathymetry and shoreline may result in numerical instabilities or imprecision. Therefore, correction coefficients for sea effects have been computed prior to the inversion, using digital elevation model of Martinique and existing detailed bathymetry (NOAA). Electric field E and magnetic field B computed at the surface of a "normal" medium (without sea, respectively E n , B n ), i.e. our 3-D resistivity model during inversion, are related to the "true" fields (with sea, respectively E t , B t ) by distortion tensors K and R using:

E t =K -1 E n (1) B t =R -1 B n (2)
K and R distortion matrices have been computed using bathymetry until a distance greater than induction length of the maximum period acquired (1000 s), i.e. a distance about 200 km around each prospected area. The maximum resolution of the model in the investigation zone is 350 m. Then, during 3-D inversion process, at each iteration during forward model computation, E and B fields are corrected of the distortion factor before comparison to real data. The eight components of the impedance tensor are tried to be fitted during the 3-D inversion process. Tipper has not been inverted.

RESULTS & INTERPRETATION OF ANSES D'ARLET GEOTHERMAL FIELD

32 magnetotelluric soundings distributed in the area of interest (the uppermost three by helicopter owing to the rugged topography) were measured during the end of November 2012. After processing (see section 3.2), data analysis showed that most of the data exhibit a 1-D behavior above 1 to 10 Hz, while the longest periods reflect more complicated structures with sometimes diagonal components one order of magnitude lower than off-diagonal ones.

We also observed a lack of coherency around 10 s which affected data quality in this dead band (Figure 1). Although surficial resistivity is generally below 50 Ωm in the investigated area, static shift effect was also observed at some stations. Before 3-D inversion, model parametrization has been designed to optimize the number of data versus the number of parameters, especially cells are merged together as a function of depth to reduce model parameters. A model of 18x16x12 cells (in the north, east and to depth direction, respectively is then adjusted to the grid data. The smaller cell of the model is 425x425m wide in the center of the model and the largest measures 2000 m at the boundary. Resistivity distribution of layer 6 (621-1081 m depth) illustrates the main geoelectric structure of the investicated area. The coastal resistive body identified on layer 5 is better defined and extends form Petite Anse (south of Morne Jacqueline) towards the Morne Réduit flanks. Its wide extent, 6 km along the NNO-SSE axis and 2.5 km perpendicularly suggests the presence of a relatively massive body. Only two conductive spots remain, one below Morne Larcher and the other below Morne La Plaine (northwards) delimited by two faults approximately orientated NE-SO. Although no MT sounding was recorded at the Anses d'Arlet, we observe a more conductive area always existing close to the hydrothermally altered zone, and disconnecting clearly the Morne Jacqueline. Note the southern well-defined boundary of the resistive body which concentrates most of the Petite-Anse springs having geothermal characteristics.

Resistivity distribution of layer 7 (1081-1851 m depth) shows the continuity with the overlaying layer 6. Northwestwards a homogeneous zone of 10-20 Ωm appears and contrasts with the massive resistive body (10-200 Ωm) along a N-S axis. The latter shows two poles because of the absence of MT data in the Anses d'Arlet area to constraint more precisely its extent. The conductive body of the Morne Larcher is still present but they are no more evidence between the hydrothermally altered zone and resistivity, suggesting that it is not so deeply rooted or that size may not be big enough to be highlighted.- The last layer 8 (1851-3051 m depth) is globally similar to layer 7. Because of increasing depth, resolution decreases, but structuration remains the same. The main axis, now N-S is generated by the association of mesh cells at depth. We note a resistivity increase for the massive body and the terrain of intermediate resistivity. Geological observations around Morne Jacqueline, together with geochemical and previous geophysical results lead to an integrated interpretation of the area, illustrated on figure 3. The coast N-S resistive body is interpreted as a shallow intrusion deeply buried, which could be related to the observed superficial dioritic dikes in the Eaux Ferrées sector of Petite Anse (Traineau et al., 2013[START_REF] Gadalia | Compléments d'exploration géothermique en Martinique : conclusions et recommandations pour les zones de la Montagne Pelée, des Anses d'Arlet, des Pitons du Carbet et du Lamentin[END_REF]. This resistive anomaly (> 50 Ωm) is strongly correlated with the density model (contour with density > 2.75 g/cm3, Figure 3). We suggest that it could act as the heat source for the current geothermal system characterized by a probable decreasing activity [START_REF] Gadalia | Compléments d'exploration géothermique en Martinique : conclusions et recommandations pour les zones de la Montagne Pelée, des Anses d'Arlet, des Pitons du Carbet et du Lamentin[END_REF]. This interpretation does not discard a potential heat source related to the Morne Larcher (346 ka). Conductive spots are associated for some of them to the internal altered structure of volcanic edifices that could have developed a local hydrothermal system during eruptions phase. Following the classical model of [START_REF] Johnston | Evaluation of electromagnetic methods for geothermal reservoir detection[END_REF], the uppermost conductive layer (thicker in the southern part of the area) in interpreted as a cap rock of the geothermal system The geophysical synthesis proposes the area of Petite Anse to be the most promising for exploration (red square of circle on Figures 2 and3). Processing (see section 3.2) and data analysis showed that most of the data exhibit a 1-D behavior above 10 to 1 Hz, while the longest periods reflect more complicated structures. Signal induced by solar activity in the longest periods was significantly lower during December 2012 resulting in noisy components of the impedance tensor in the MT dead band (3 to 32 seconds) affecting resolution of the model at depth. Figure 4 shows typical results collected during the survey. A 3-D inversion has been run [START_REF] Hautot | The deep structure of the Baringo Rift basin (central Kenya) from 3-D magneto-telluric imaging: Implications for rift evolution[END_REF][START_REF] Hautot | 3-D magnetotelluric inversion and model validation with gravity data for the investigation of large igneous provinces[END_REF] taking into account coast and sea effects. In order to resolve at best the complex structure of Montagne Pelée, model parametrization has been designed to optimize the number of data versus the number of parameters, especially cells are merged together as a function of depth to reduce model parameters. A model of 22x23x12 squared cells (in the north, east and to depth direction, respectively is then adjusted to the grid data. The smaller cell is 500m wide in the center of the model and the largest measures 750 m at the boundary. MT investigations on the Montagne Pelée reveal a first surficial resistive (20 <  < 2000 Ωm) layer 200-500 m thick attributed, at least for a part, to volcanic deposits related to the last eruption cycles of the Montagne Pelée. This unit overlays a thick (500-1000m) and massive conductive ( < 15 Ωm) which follows approximately the topography. Few major discontinuities apparently affect this layer that could be related with the long term volcanic and tectonic history of the Montagne Pelée. The conductive layer is mainly attributed to older deeply altered volcanic series (pre-Pelée?). Its high conductivity indicates high clay content which is unfavorable (at large scale) to fluid transfer and may play a caprock role above potential geothermal system(s). Despite the masking effect of this layer, below the conductive layer, MT data reveal a third and thick more resistive layer that can be interpreted in some parts as the potential geothermal reservoir. Deeper, the 1 km thick resistive layer is poorly constrained by the data although refraction seismic data (re-processed from [START_REF] Eschenbrenner | Coupes sismiques des structures superficielles dans les Petites Antilles -II : Martinique. Pageoph[END_REF] show a good morphological correlation with it, resulting of increasing velocity. Slightly southwestwards from the summit, we delimited a major discontinuity (D on figure 5, -200 m) based on geophysical data (magnetotelluric and gravity data). Because two different geochemical signatures characterize the springs of Mitan-Picodo Rivières (Cl-Na) and Rivière Chaude (HCO 3 -Cl-Na), we infer the potential existence of two separatedat least partlygeothermal reservoirs. This discontinuity could act at depth either as a potential barrier that would favor separated geochemical evolution or as drain for preferential fluid circulation.

Resistivity layers from -400 down to -600 m below sea level point out a highly conductive axis NW-SE orientated, running parallel to the Caribbean sea, and matching gravity lineaments detected on the gradient of the Bouguer anomaly (corrected for a density of 2.0). We may also note on layer -600 b.s.l. that the scars of the latest lateral collapse LC3 (located on figure 5, 600m layer, 9 ka, Le [START_REF] Friant | Large-scale flank collapse events during the activity of Mount Pelée, Martinique, Lesser Antilles[END_REF] circumscribes precisely a 5 to 10 Ωm conductive area with higher resistivity on both sides belonging to the NW-SE main axis. This observation is thought to reflect lateral flank collapse processes that excavate one part of the flank of the volcano, including a substantial thickness of the characteristic and general conductive layer found in many volcanic edifices [START_REF] Descloitres | Electrical structure beneath the eastern collapsed flan of Piton de la Fournaise volcano, Reunion Island: Implications for the quest for groundwater[END_REF][START_REF] Nurhasan | Two electrical conductors beneath Kusatsu-Shirane volcano, Japan, imaged by audiomagnetotellurics, and their implications for the hydrothermal system[END_REF][START_REF] Monteiro Santos | An audio-magnetotelluric investigation in Terceira Island (Azores)[END_REF][START_REF] Coppo | A deep scar in the flank of Tenerife (Canary Islands): Geophysical contribution to tsunami hazard assessment[END_REF]. These catastrophic events deeply affect the stability of the edifice providing potential new preferential pathways for fluids migration and re-formation of the conductive layer which may later acts as cap rock. Although quasi ubiquitous, the conductive layer has a complex morphology that is difficult to connect to a simple a global model with geological and known geothermal manifestations. Particularly, it is difficult to define precisely areas where the conductive layer could act with certainty as cap rock as expected above a well-developed geothermal system. In the classical model of conductive caprock overlying a geothermal, more resistive system [START_REF] Johnston | Evaluation of electromagnetic methods for geothermal reservoir detection[END_REF], the upflow zone is located at the top of the reservoir [START_REF] Crosby | Bulls-eye! -Simple resistivity imaging to reliably locate the geothermal reservoir[END_REF]Usher et al., 2000). In the case of the Montagne Pelée, such typical zone is well highlighted on the -200 m map (Figure 5) where the resistive is surrounded by the conductive layer. The apex of this resistive body is approximately located at altitude 0 m and below the southeastern flank of the summit of the edifice.

To better understand the Montagne Pelée structure, we integrated available SP (Zlotnicki et al., 1998) and refraction seismic (reprocessed form [START_REF] Eschenbrenner | Coupes sismiques des structures superficielles dans les Petites Antilles -II : Martinique. Pageoph[END_REF] data from previous works (Figure 6) along profile B-B'. Some attempts were made to correct the topographical electrokinetic effect on SP in order to better evidence positive anomalies. . Re-interpretation of seismic refraction data show a high velocity anomaly (Vp = 3000-4000 m/s) between 1000 and 3000 m depth below the Montagne Pelée summit (Figure 6). We thus observe on the SW flank of the volcano, a correlation between, SP positive anomalies (A, B, C), resistive anomalies within the main condcutive layer (R1, R2, R3), dense anomalies of the gravity field (G0, G1 and G2) and a high velocity anomaly at relatively shallow depth as illustrated in figure 6. This configuration of dense and rapid structures may indicate the presence of shallow solidified intrusion possibly related to the present volcano as indicated by the recent lava dome outcropping close to G2 anomaly.

Location of both G1 and G2 dense structures coincident with the uphill scar of the last lateral collapse (LC3 on Figure 5, 600 m), indicate that both structures may be linked. The massive nature of these body and their recent age, both resulting in weak alteration and low clay content can explain the high resistivity anomalies observed and may manage a potential permeability which favours upward fluid flow and associated positive SP anomalies.

This whole set of observations more or less converges towards a general schema of which make it possible the presence of a geothermal system heated by quite superficial and recent intrusion, cap by a general and thick conductive layer which is cut by fault and intrusion on the volcano southwest flank which leave hydrothermal fluid weakly flowing upward and outcrop as the known scarce surface manifestation

SURFICIAL RESISTIVITY

During the MT survey, several soundings were affected by static shift effect, a physical galvanic effect causing apparent resistivity curves to be shifted relative to each other from an unknown scalar factor, while the phase remain unaffected. The time-independent nature of static shift means that there is no impedance phase associated with the phenomenon [START_REF] Simpson | Practical Magnetotellurics[END_REF]. Many publications deals about static shift, techniques to recover "true" resistivity (for ex : Bahr, 1988;[START_REF] Sternberg | Correction for the static shift in magnetotellurics using transient electromagnetic soundings[END_REF][START_REF] Jones | Static-shift of magnetotelluric data and its removal in a sedimentary basin environment[END_REF] and the need or not to correct it before 3-D inversion [START_REF] Meju | Joint inversion of TEM and distorted MT soundings: some effective practical consideration[END_REF].

At the beginning of 2013, within the framework of the MarTEM project, a heliborne TEM survey carried out by SkyTEM covered the whole Martinique. With the aim to illustrate the ability of the 3-D inversion scheme [START_REF] Hautot | The deep structure of the Baringo Rift basin (central Kenya) from 3-D magneto-telluric imaging: Implications for rift evolution[END_REF][START_REF] Hautot | 3-D magnetotelluric inversion and model validation with gravity data for the investigation of large igneous provinces[END_REF] Right plots show variation between the two models (log scale).

CONCLUSIONS

This paper deals with a part of the geophysical results collected during the last multi-disciplinary geothermal exploration phase of Martinique (2012Martinique ( -2013)). We present both the 3-D resistivity models of the Anses d'Arlet and Montagne Pelée geothermal fields resulting of separate 3-D magnetotelluric inversions including sea effect. We point out a potential geothermal resource on each field/area characterized by local specific structures. Complementary geological and geochemical information is provided to help the reader understanding each local context.

On the Anses d'Arlet sector, MT and gravity data reveal a massive dense and resistive body, running parallel to the Caribbean sea and elongated in the NNW-SSE direction. It is interpreted as a N-S trending deep intrusion related to volcanic activity that affected the area during the last My, and could act as heat source of the high temperature geothermal system. A cap rock structure is interpreted at shallow depth. Although many geological and geochemical elements seem to indicate that the geothermal system is currently becoming less active and extensive in space, its southern boundary, which concentrates most of the surface hydrothermal manifestations is proposed to be the best place for further exploration wells.

The geothermal field related to the Montagne Pelée is characterized by a quasi ubiquitous bell-shaped conductive layer and with quite monotonic resistivity distribution making interpretation difficult in terms of geothermal targets. According to geological and geochemical data, we were able to distinguish a potential main discontinuity that could be related to two separated geothermal reservoirs approximately located below the summit of the edifice and southwestwards. Resistivity distribution of the conductive layer seems to be closely connected with edifice evolution and flank collapse events.

Beside these interesting deep structures, we demonstrate, after analyzing the results of the recent heliborne TEM survey covering the whole Martinique island, that surface resistivity distribution obtained from 3-D inversion reproduce faithfully the resistivity distribution observed by TEM. In spite of a very different sampling scale, this comparison illustrates the ability of 3-D MT inversion to take into account and reproduce static shift effects in the sub-surface resistivity distribution.

Figure 1 :

 1 Figure 1: MT soundings 05 (left) and 26 (right) showing data (off-diagonal Zxy and Zyx components in red and blue, respectively) and model responses in green. Upper plots display resistivity and lower plots phase for the same components.

Figure 2

 2 Figure2shows four layers of the resistivity model of the Anses d'Arlet obtained through 3-D inversion (very surficial resistivity distribution is presented on Figure7together with TEM data for comparison). Resistivity distribution of layer 5 (351-621 m depth) displays the emergence of a resistive costal body, globally orientated NNO-SSE, more evidenced around Morne Champagne (MCh, NNO). Inland, conductive areas mostly trail off excepted below Morne Larcher where a thickening is observed. Note the persistence of the conductive body centered on the hydrothermally altered zone between the Anses d'Arlet and Petite Anse (MT soundings 22 and 23), and elongated in the main tectonic direction of Martinique. Further north, a conductive body appears below Morne Réduit, and could be interpreted as for Morne Larcher to the altered core of the volcano.

Figure 2 :

 2 Figure 2 : a) Location map of the investigated area on geological background (Westercamp et al. 1900) and resistivity models at 351 to 3051m deep(b, c, d, e). Superimposed elements are indicated in the legend.

Figure 3 :

 3 Figure 3 : 2-D profiles crossing the area of main interest for further exploration drilling. Resistivity profile (above) is extracted from the 3-D resistivity model and the density profile (below) from the 3-D density model.

Figure 4 :

 4 Figure 4 : MT soundings 19 (left) and 47 (right) showing data (off-diagonal components in red and blue) and model responses in green. Upper plots display the resistivity and lower plots the phase for the same components.

Figure 5

 5 Figure5presents the results of 10 layers extracted from the 3-D resistivity model by altitude to illustrate characteristic changing resistivity with depth. From 800 m down to 200, a typical sub-circular resistivity distribution with radially increasing resistivity from the center highlights the conductive core of the volcano. The central conductive (< 5 Ωm) structure is interpreted as altered core thought to be related with the development of hydrothermal alteration during and post eruptive phases. Going down (200 m), the conductive layer extends forming a bell-like geometry whose top follows approximately the topography of the edifice. The presence of the conductive layer at higher altitude on the northeastern part of the volcano is strongly related to the topography and its historical evolution. The three lateral collapses that affected the southwestern flank of the Montagne Pelée during the last 100 ka (Le[START_REF] Friant | Large-scale flank collapse events during the activity of Mount Pelée, Martinique, Lesser Antilles[END_REF] deeply excavated the structure of the volcano, thought to be initially more homogeneous. At sea level, we observe the development of a conductive body NNW-SSE just to the south of the summit, rotating to NW-SE and growing 200m below sea level. Note the changes between 0 and 400 m below sea level where the conductive layer has moved southwestwards leaving the upper part of the southwestern flank and the eastern part of the volcano with intermediate resistivity (10-50 Ωm).

Figure 5 :

 5 Figure 5 : Ten horizontal slices extracted from the 3D resistivity model of the Montagne Pelée stratovolcano. Superimposed information is detailed in the legend.

Figure 6 :

 6 Figure 6 : Section B-B' (indicated on Figure 5, -400 m map). A) SP original data (Zlotnicki et al., 1998) and corrected for topographical variations, and vertical gradient data (correction density = 2.0). b) Resistivity section extracted from the 3-D resistivity model ( R0, R1, R2 identify more resistive bodies within the main conductive layer) with superimposed velocity model (re-processed from Eschenbrenner et al., 1980). Magenta lines display the top of the denser bodies identified on c). Dashed vertical black line corresponds to the main discontinuity showed on Figure 5. Dashed vertical white lines illustrate discontinuities in the conductive layer. c) Density section extracted from 3-D density model with velocity model. G0, G1, G2 identified denser bodies.

  to take into account and reproduce static shift effects.to reproduce surficial resistivity distribution, we compared results of the 3-D resistivity model with the results of the TEM resistivity (1-D stitching). To homogenize TEM distribution, TEM resistivity data are first krigged by layer over the investigated area. The median resistivity of this new grid is then extracted above each MT cells and compared with MT resistivity (Figure7, right plots). Cells that no contain raw TEM information or that are located too far away from fly lines are removed from the statistics. Results show that 82 % of the remaining cells have a maximum resistivity variation below +/-400% for the Anses d'Arlet prospect, 79 % for the Montagne Pelée.These results demonstrate that in spite of a fully opposite sampling scale (53 MT soundings and more than 20000 TEM for the Montagne Pelée and 32 MT soundings and more than 6200 TEM), surficial resistivity distribution is quite faithfully reproduced by 3-D inversion.

Figure 7 :

 7 Figure 7 : MT versus TEM resistivity for superficial layers. Above are represented the data of the Anses d'Arlet (AA) geothermal province and below from the Montagne Pelée (MP) geothermal province. Left plots display resistivity from the 3-D MT models for both sector at the indicated depth. Center plots display resistivity of the TEM data (see text for details).Right plots show variation between the two models (log scale).

ACKNOWLEDGEMENTS

The authors wish to thanks the FEDER and ADEME institutions, the Conseil Régional and the Syndicat Mixte d'Electricité de la Martinique (SMEM) for funding this project.


Hana Krichene
email: hana.krichene@inria.fr

Mouna Baklouti
email: mouna.baklouti@enis.rnu.tnmohamed.abid@enis.rnu.tn

Mohamed Abid

Philippe Marquet
email: philippe.marquet@lifl.frjean-luc.dekeyser@lifl.fr

Jean-Luc Dekeyser

Communication-Computation overlap in massively parallel System on Chip

This paper presents the concept and the implementation of a communication-computation overlap in massively parallel System on Chip (mpSoC). This paradigm allows to decrease the execution time of parallel programs using specific strategies in the programming level and partially decoupled system control in the hardware level. The related experiment in VHDL language for system design and in assembly language for case study implementation are described.

I. INTRODUCTION

Many modern application domains are concerned by the conjunction of parallel algorithms and high computing resources. They include signal and image processing applications such as software radio receiver, sonar beam forming, or image encoding/decoding. Furthermore, the implementation of the system on a single chip will be of prime interest for those applications that also require some degree of embeddedness.

Nowadays, SoCs, used for massively parallel computing, consist of many processing units connected by communication networks where running the data to be processed and the results of various intermediate calculations. When it comes to ten thousand of processing units, the network must be efficient to avoid slowing down the whole system. The main complaint of these systems against communication networks is that they consume computation time: they must cut messages into packets in the sender side and reassemble it in order in the receiver side, without error on path. A communicating program on the network will therefore spend time to make these communications, instead of using that time to calculate. Given that the cost of power consumption and area occupation, this mess is unsatisfactory.

Several models of communication-computation overlap, e.g., Message Passing Interface (MPI) [START_REF] Gropp | Mpi: The complete reference[END_REF], Parallel Virtual Machine (PVM) [START_REF] Geist | Pvm parallel virtual machine: a users guide and tutorial for networked parallel computing[END_REF], multi-threaded MPI [START_REF] Jiayin | Overlapping communication and computation in mpi by multithreading[END_REF] and Hybrid MPI/SMPS [START_REF] Marjanovic | Overlapping communication and computation by using a hybrid mpi/smpss approach[END_REF] were implemented to reduce execution times and improve speedups. MPI and PVM rely on non-blocking elementary communications, asynchronous send and receive primitives as provided by well-known message-passing libraries while the multi-threaded MPI model proposes MPI point-to-point operations in order to overlap communication and computation. Compared to single-threaded active polling communication used by most of the implementations, multithreaded provides more parallelism. Hybrid MPI/SMPS uses shared memory programming model to introduce asynchrony necessary to overlap communication and computation. If the architecture allows the integration of these high level program software solutions, the speed of parallel applications could double. But they still be limited by the hardware complexity to adopt these approaches.

In this work, we define a new mechanism to amortize the time spent by the communication phase. In fact, processing unit should indicate that it wants to send or receive something, then return directly to his computations while the network would transfer requested. This is called overlapped communications by computation.

The objective of this paper is to reconsider the interest and the feasibility of communication-computation overlap into massively parallel System on Chip. This model allows for communications operations in the background while the program continues to run. It is thus possible to hide a portion of any communication costs.

The next section presents our proposed overlapping communication-computation model. Section 3 discusses the structure implementation and the performance evaluation through the matrix multiplication case study.

II. COMMUNICATION-COMPUTATION OVERLAP MODEL

Distinguish calculation stages of those communication needs the separation of these two stages in different blocks. This repartition should be provided by the designer at programming level. Then, the overlapped execution of these blocks will be done in parallel according to the program description. Separate communication from computation provides an opportunity to ensure the autonomy execution in mpSoC. To ensure this feature while avoiding the centralized control, we define second hierarchical control level, over than the master controller, achieved by the slaves controllers. We will discuss this structure in the next section.

A. Hardware design: Master-Slave control structure

A novel control structure was proposed for the massively parallel System-on-Chip, referred to as master-slave control structure [START_REF] Krichene | Master-slave control structure for massively parallel system on chip[END_REF]. Its concept departs from the centralized configuration and instead of a uni-processor master controlling a set of parallel Processing Elements (PEs), the master cooperates with a grid of parallel slave controllers which supervises the activities of cluster of PEs. We define, as shown in fig. 1, the hardware implementation of this configuration in massively parallel system: The hardware architecture is composed of a single MCU and multiple Slave controllers (SCUs) combined with local processing element (PE) (or a cluster of 16 PEs), known collectively as Nodes. The MCU and SCU array are connected through single level hierarchical bus and the SCUs are connected together through X-net interconnection network. This network is clocked synchronously with the SCUs and respectively with the PEs. SCU controllers in the grid care for the instruction execution activities that involve a large degree of parallelism and the communication activities that need to coordinate all the PEs in the grid. Note that the SCU controllers do not limit performance; they do not become a sequential bottleneck. They participate only in the controlling and scheduling of very large groups of PEs execution at a time.

Using master-slave control structure, the PEs in massively parallel system can execute independently and then can communicate synchronously. Separate communication phase from computation one, not only allows the programmer to optimize distinct processors for their intended tasks, but also gives the ability to overlap these two phases which reduces the system execution time.

B. Software program: Master-Slave program management

Based to the decoupled control design presented in the previous section, the programmer can easily manage the masterslave program to overlap communication by computation stage. Therefore, the basic idea to implement this paradigm is to divide the principal program into small blocks of parallel instructions, called Slave Program (SP), and send these blocks to the activated PEs of the system. Then, according to a predefined mask, the SCU sends the order of begin execution. In parallel to computation, the master manages a synchronous inter-node communication, and so on until the end of the program.

To implement this master-slave program we have defined some instructions as shown in fig. 2: Master side

• initial.mask: master controller informs the activated node, using broadcast with mask technique [START_REF]Broadcast with mask on a massively parallel processing on a chip[END_REF], that PEs must start the parallel execution of SP blocs. • SP.execution: slave controller orders the parallel execution of SP blocks. The program in SP blocks is the same for all PEs, but the control flow may be different in each node according to the mask. Each SP block starts with Begin-instruction and finishes with End-instruction.

When SP blocks are executed, the communication is established. GO signal is activated. • wait.ortree: master controller waits for the SP execution end of all the active PEs to start a new parallel cycle. This end execution is controlled by the ortree signal [START_REF] Krichene | Master-slave control structure for massively parallel system on chip[END_REF].

Slave Side

• wait.go: slave PE waits for go signal to start SP execution. • end.SP: slave PE informs the master that the SP execution is finish by sending the end signal.

We require also that the SP block code contains only memory reference of local PE or local register in the relative node and no reference to the global memory or the local memory of the neighbor PE.

Referring to these requirements and specifications, it is possible to perform communication operations in the background while the slave program continues to run, as it is shown in fig. 2. However, This feature is allowed only if the execution time of SP block is higher than time spent for communication phase. One of the basic computational kernels in many data parallel codes is the multiplication of two matrices (C=AxB). For this application we have implemented two architectures: the first one is a 64 PEs and the second one is a cluster of a 64 (PE+ muladd IP). These computation units are arranged in 8x8 node-grid with each PE has a 64bytes data stack . To perform multiplication, all-to-all row and column broadcasts are performed by the nodes. The following codes are executed by the master and all slave-PEs simultaneously (Fig. 3): Fig. 4 shows the Virtex6 ML605 FPGA resource occupation and the execution time results. We validated that the architecture based on PE+muladd IP is more efficient but occupies a large area on the chip. In fact, the operation performed in FC16 [START_REF] Haskell | A VHDL forth core for FPGAs[END_REF] PE needs 19 clock cycles but only one clock cycle in muladd IP. However, the space on a chip is increased because the additionally area occupied by muladd IP other than the PE area. On the other hand, we see that the execution time when using the overlapping communication-computation paradigm is decreased in the both implementations (only PE and PE+muladd IP). This is due to the fact that communications are done when the processing units execute their locally slave program, so that we introduce a second parallelism level in addition to the massively parallel execution.

IV. CONCLUSION

This paper presents the paradigm of communicationcomputation overlap, not only in the programming level but also in the architectural level. The decoupled master-slave control structure allows the independent execution of parallel slave programs and the globally synchronous communication. Separating the computation control from the communication one allows the designer to perform easily programming strategies that overlap these two phases. Using this paradigm, we save around 5% of the globally execution time. It is thus possible to hide a portion of any communication costs.

Fig. 1 .

 1 Fig. 1. Hardware prototype of Master-Slave control structure for mpSoC

Fig. 3 .Fig. 4 .

 34 Fig. 3. Overview of matrix multiplication program