Jean-Baptiste Jorcin

Christine Blanc

Nadine Pébère
email: nadine.pebere@ensiacet.fr

Bernard Tribollet

Vincent Vivier

Galvanic Coupling Between Pure Copper and Pure Aluminum Experimental Approach and Mathematical Model

The corrosion behavior of a pure aluminum/pure copper couple in a weakly conductive sulfate solution was investigated. Potential and current distributions on the surface of the model couple at the beginning of immersion were obtained by solving the Laplace equation using a finite element method ͑FEM͒ algorithm. The potential distribution predicted by the calculations was checked using a Ag/AgCl microreference electrode. A good agreement was found between experimental and theoretical results. It was shown that the reaction occurring at the copper electrode was oxygen reduction, while aluminum remote from the Al/Cu interface remained in the passive state. Moreover, calculations predicted a large cathodic current, related to an increase in oxygen reduction, restricted to copper at the Al/Cu interface. This led to a local pH increase reaching values higher than 9, allowing the dissolution of aluminum to occur close to the interface. Combining these data with optical and scanning electron microscope observations after 24 h of immersion in the sodium sulfate solution allowed a three-step mechanism to be proposed to explain the corrosion damage, and particularly the presence of a copper deposit on the aluminum surface, some distance from the Al/Cu interface, a phenomenon currently observed in commercial copper-rich aluminum alloys. DOI: 10.1149/1.2803506

There is significant interest in the corrosion behavior of Cucontaining Al alloys, such as the 2024 aluminum alloy ͑AA2024͒, which remain of importance for aerospace applications due to their high strength-to-weight ratio. Consequently, the corrosion susceptibility of 2XXX Al-alloy series has been studied for a long time. [1][2][3][4][5][6][7][8][9][10][11][12] The corrosion of AA2024 is mainly caused by the heterogeneous microstructure of the alloy, which has been intentionally developed to improve the mechanical properties of the material.

Among the different types of intermetallic particles present in AA2024 alloy, both coarse Al 2 CuMg particles, also called S-phase or Al-Cu-Mn-Fe particles ͑size: 1-30 m͒ and smaller Al 2 CuMg particles ͑about 200 nm͒ present in the grain boundaries have been widely studied due to their reactivity that constitutes initiation sites for corrosion. [13][14][15][16][17][18][19][20] Even though the influence of copper-rich intermetallics on the corrosion resistance of aluminum alloys was clearly shown, the mechanisms explaining the dissolution of the intermetallics, the copper enrichment, and the pit nucleation at these sites are still not clearly explained. [13][14][15][16][17][18][19][20][21][22][23] Authors usually refer to a galvanic coupling phenomenon between the particles and the surrounding matrix. Zhu and Van Ooij 24 showed that anodic Al 2 CuMg particles dealloyed both Al and Mg ͑with a more severe dealloying of the Mg͒ when immersed in a neutral 0.6 M sodium chloride solution. Simultaneously, a significant dissolution of the Al matrix surrounding the coarse particles was also observed. This mechanism is similar to that proposed by other authors. 15,18,25 In previous studies, Blanc et al. 13,14 discussed a slightly different mechanism, in which the copper-rich intermetallic reactivity was studied in a 0.1 M sulfate solution and described by a three-step process consisting of homogeneous dissolution of the particle, copper redeposition, followed by local dissolution of the matrix surrounding the particle. These phenomena have been widely investigated by electrochemical measurements in combination with scanning electron microscopy ͑SEM͒ and electron dispersive spectroscopy ͑EDS͒. However, conventional electrochemical methods lack spatial resolution, while microscopy techniques do not allow the in situ chemical composition to be determined; these are crucial for the understanding of the role of the intermetallic particles on the corrosion process. Consequently, many authors have used atomic force microscopy ͑AFM͒, which provides high lateral resolution to study dissolution phenomena occurring in the S-phase. 18,[25][26][27] To understand the corrosion behavior of copper-rich aluminum alloys, another alternative is to study the electrochemical behavior of model alloys that are representative of the different metallurgical phases. 22,28 For example, thin films of Al-Cu alloys can be easily obtained using magnetron sputtering and are useful for such studies. 29,30 However, these alloys are mostly multiphase alloys and thus their corrosion behavior remains intricate.

In the present study, a simple system was considered consisting of a pure aluminum/pure copper ͑Al/Cu͒ couple. Sample immersion for 24 h in a 10 -3 M Na 2 SO 4 solution was followed by observations with an optical and a scanning electron microscope to reveal the corrosion morphology. Particular attention was paid to the Al/Cu interface. Measurements of local potentials were carried out and compared with results obtained from calculations based on the resolution of the Laplace equation to describe the first moments of immersion. On the basis of the experimental approach and the mathematical model, the interpretation of the corrosion phenomena occurring at the Al/Cu interface is proposed.

Experimental

Samples.-The sample consisted of a pure aluminum/pure copper couple. A cylinder of pure aluminum ͑99.999 wt %͒, provided by Alfa/Aesar, was drilled in its center and a cylinder of pure copper ͑99.9 wt %͒, provided by Goodfellow, was then introduced by force into the hole. The assembly of the two materials gave a perfectly joined interface, avoiding any crevice corrosion due to surface defects. The diameters of the two cylinders were chosen to obtain a ratio aluminum surface area/copper surface area of 10 ͑external diameters were equal to 10 and 3.15 mm for aluminum and copper bars, respectively͒. The electrode was then embedded in an epoxy resin so that a disk electrode was obtained at the extremity.

Before immersion in the electrolyte, the disk electrode was mechanically polished with SiC papers up to 4000 grade and ultrasonically cleaned with ethanol, then with distilled water. The electrolyte was a 10 -3 M Na 2 SO 4 solution prepared with analytical grade chemicals in contact with an air atmosphere and at room temperature.

Local electrochemical measurements.-Local electrochemical measurements were performed using a homemade device. The local potential variations were measured using a 100 m diameter Ag/AgCl microelectrode. It consisted of a Ag wire ͑100 m in diameter, Goodfellow͒ laterally insulated with cataphoretic paint and sealed in a glass capillary with epoxy resin. The AgCl was deposited in a 2 M KCl solution under potentiostatic oxidation of the Ag electrode at 0.4 V/SCE. The microreference electrode was moved with a three-axis positioning system ͑UTM25, Newport͒ driven by a motion encoder ͑MM4005, Newport͒ allowing a spatial resolution of 0.2 m in the three directions. A homemade analog differential amplifier with both variable gain and high input impedance was used to record the local potential. All measurements were performed at the corrosion potential E corr using a Keithley 2000 digital multimeter. The experimental setup was controlled with homemade dataacquisition software developed under a LabView environment.

Surface characterizations.-The samples were observed before and after 24 h of immersion in a 10 -3 M Na 2 SO 4 solution by optical microscopy with an Olympus PMG3 microscope. Observations with a scanning electron microscope were also performed using a Leo 435VP apparatus to obtain a better description of the corrosion morphology, particularly at the Al/Cu interface.

Theoretical Description of the Current and Potential Distributions

Mathematical model.-Numerical simulations were performed to obtain a description of the potential and current distributions on the surface of the disk electrode and in the surrounding electrolyte solution at the beginning immersion. Figure 1 gives a schematic representation of the disk electrode used in this study.

The potential ⌽ in the solution surrounding the electrode is governed by the Laplace equation

ٌ 2 ⌽ = 0 ͓1͔
Using cylindrical coordinates ͑r,,z͒, Eq. 1 can be expressed as

1 r ץ ץr ͩr ץ⌽ ץr ͪ + 1 r 2 ץ 2 ⌽ ץ 2 + ץ 2 ⌽ ץz 2 = 0 ͓2͔
where z is the normal distance to the electrode surface, r is the radial coordinate, and is the azimuth. The cylindrical symmetry condition requires the geometry to be invariant under rotation about the y axis, i.e.,

ץ⌽ ץ = 0 ͓3͔

The combination of Eq. 2 and 3 yields the governing equation in a two-dimensional ͑2D͒ domain as

ץ 2 ⌽ ץr 2 + 1 r ץ⌽ ץr + ץ 2 ⌽ ץz 2 = 0 ͓4͔
On the surrounding insulator and far from the electrode surface, the boundary conditions are given by

ͯ ץ⌽ ץz ͯ z=0 = 0 at r Ͼ r 0 ͓5͔
and

⌽ → 0 as r 2 + z 2 → ϱ ͓6͔
Under the assumption of a kinetic regime, the current density at the electrode surface can be expressed as

I = -ͯ ץ⌽ ץz ͯ z=0 ͓7͔
where is the electrolyte conductivity. The boundary conditions were determined from experimental measurements performed individually on each material. On the aluminum part of the electrode, the boundary condition corresponds to an anodic current of 1 A cm -2 . On the copper electrode, both oxygen reaction and copper dissolution were taken into account. The global current can be expressed as The simulations were performed using a finite element package Femlab in a 2D axial symmetry. The mesh size was refined to obtain a numerical error lower than 0.1% evaluated from the net current of the system, which is the sum of the current passing through the copper and aluminum electrodes.

I Cu = k Cu expͩ 2.303 b Cu ϫ ⌽ͪ -k O 2 exp ͩ -2.303 b O 2 ϫ ⌽ ͪ ͓8͔
Results of the simulations.-Figure 2 shows the potential distribution calculated on the electrode surface along the electrode radius with the electrolyte conductivity as a parameter. Independently of the electrolyte conductivity, the general shape of the potential distribution remains similar. The potential is seen to be constant over the copper electrode; it strongly increases immediately after the Cu/Al interface and reaches a maximum value over the aluminum electrode near the Al/insulator interface just before a slight decrease. Such behavior is fully consistent with the boundary conditions determined for the calculation. Figure 2 also shows that the variations of the potential along the disk-electrode radius are higher for a low conductivity of the electrolyte which, from a practical point of view, made these variations easier to detect. Figure 3 shows the influence of Ag/AgCl probe position ͑i.e., readings at two distances from the electrode surface͒. Calculations were performed in the case of a weakly conductive electrolyte to highlight the influence of the position of the probe on the potential distribution. No significant difference was observed between the two measurements; when the probe was withdrawn from the electrode surface, the sudden variation of the potential at the Al/Cu interface was barely reduced by comparison to measurements performed on the disk electrode surface itself. Moreover, the amplitude of the potential variations was smaller when the electrode was far from the substrate.

The resolution of the Laplace equation also allows the distribution of the current on the disk electrode surface to be determined: both normal and radial current distributions were calculated as shown in Fig. 4. The normal current ͑Fig. 4a͒ on aluminum is constant and corresponds to the boundary condition chosen, while for the copper electrode, the calculations show a cathodic current distributed along the electrode radius. Thus, calculation led to results in good agreement with the initial hypothesis, a passive behavior for aluminum and a cathodic current related to oxygen reduction on copper. Furthermore, potential and current distributions revealed a particular behavior at the Al/Cu interface. Figure 4a shows a high cathodic current on copper at the Al/Cu interface corresponding to an increase of oxygen reduction on copper

O 2 + 2H 2 O + 4e -→ 4 OH - ͓9͔
This suggests that local variations of the chemical composition of the electrolyte can occur from the beginning of immersion with a local alkalinization of the electrolyte at the Al/Cu interface. Figure 4b also shows an increase of the radial component of the current over the aluminum in the vicinity of the interface between the two metals. The variations of the radial contribution are linked to the potential distribution over the aluminum, whereas for a constant value of the potential at the electrode surface ͑for instance over the copper electrode as shown in Fig. 3͒ the normal contribution of the current varies from the electrode center to the Cu interface ͑Fig. 4a͒ and the radial current contribution is zero. These results suggested a particular evolution of the interface morphology of the disk elec-trode after immersion in the electrolyte related to corrosion phenomena restricted to the Al/Cu interface.

Results and Discussion

Morphology of the Al/Cu interface after 24 h of immersion in a 10 Ϫ3 M Na 2 SO 4 solution.-The Al/Cu interface was observed using optical microscopy and SEM before and after immersion in the sodium sulfate solution. Figure 5 shows an optical micrograph of the interface before immersion ͑copper on the right and aluminum on the left of the micrograph͒. The two materials are perfectly joined without any defects observable at the interface. After 24 h of immersion in the sulfate solution ͑Fig. 6a͒, two main differences were observed in the vicinity of the interface by comparison with the observations performed before immersion.

1. A deep crevice was formed at the Al/Cu interface. Close to the crevice, aluminum presents a bright color, while, except in this zone, it seems to be covered by an oxide film over its whole surface, which is in good agreement with previous hypotheses concerning a passive state for aluminum in the couple. The bright color suggests that the crevice was at least partially related to the dissolution of the aluminum close to the interface.

2. An orange-brown ring of about 50 m thick occurred all around the Al/Cu interface. It was located on the aluminum material at a distance of about 150 m from the interface. The ring was analyzed by EDS as being copper and is attributed to the redeposition of copper coming from the dissolution of the Al/Cu interface. 31,32 The observations also show that the corrosion processes were restricted to the Al/Cu interface and only concerned the interface or the zone close to the interface ͑copper ring at 150 m from the interface͒. No other domain of the sample exhibited corrosion damage, whether on aluminum or copper. Moreover, Fig. 6b clearly shows that the corrosion phenomenon is independent of the size of the interface involved in the reaction process. Similar phenomena were observed around a 10 m particle of copper which became embedded in the aluminum part of the electrode during mechanical polishing. The particle and its surrounding aluminum also behave as a microcouple. The formation of a crevice is seen at the Al/Cu interface and copper redeposition occurred on the aluminum part not far from the copper particle. Furthermore, these results are in good agreement with those obtained on AA2024, in which Al 2 CuMg particles act as copper sources after immersion in sulfate solutions with or without chlorides. 13 Further SEM observations allow the zone close to the interface to be more accurately described. In Fig. 7a andb, the gray zone corresponds to aluminum while the white zone corresponds to copper. The deep crevice previously observed using optical microscopy was more accurately seen on the SEM micrographs. Moreover, aluminum was found to be corroded close to the Al/Cu interface, with a dissolution depth decreasing from the interface to some micrometers from the interface ͑Fig. 7a͒. Except in this zone, aluminum was not corroded. Similar observations were performed for copper, for which no corrosion damage on the copper disk except at the Al/Cu interface was seen ͑Fig. 7b͒. Accurate analysis of the topography of the copper at the interface gave striking results for the corrosion process. At the Al/Cu interface, copper formed a kind of ledge, with no corrosion damage on the copper disk surface but a strong dissolution of copper at the Al/Cu interface, suggesting the occurrence of crevice corrosion of copper. 33 Comparison between theoretical calculations and experimental observations.-Observations of the morphology of the Al/Cu interface after immersion in the electrolyte suggest that corrosion occurs in different steps. The strong dissolution of aluminum at the Al/Cu interface partially leads to the formation of the crevice, as observed. Moreover, it can be assumed that, in the crevice, the chemical composition of the electrolyte varies significantly ͑occluded zone͒, leading to local chemical conditions which favor crevice corrosion of copper. The localized corrosion of copper explains the formation of the copper ring observed 150 m from the interface. Thus, the initial step of the corrosion process appears to be the dissolution of aluminum at the Al/Cu interface, while aluminum remains passive on the remaining surface of the electrode. This is in good agreement with the current distribution previously shown from the numerical calculations ͑Fig. 4͒. A strong cathodic current density, related to increased reduction of the oxygen, was predicted on copper at the Al/Cu interface. Such a reaction should lead to a local alkalinization of the electrolyte. 32 To confirm this assumption, bromothymol blue was added to the 10 -3 M Na 2 SO 4 solution ͑pH 5.7͒ before immersion. This chemical turns blue when the pH exceeds a value of 7.6. Figure 8 shows a blue ring all around the disk electrode at the Al/Cu interface. The blue ring was formed only a few minutes after immersion of the disk electrode, indicating a fast and local increase of the pH of the electrolyte. This result shows that the kinetics of the cathodic reduction of oxygen was strong enough at the Al/Cu interface to induce a significant variation of the pH in the electrolyte.

Figure 9 shows the potential distribution measured in the 10 -3 M Na 2 SO 4 solution at the beginning of immersion at the surface of the disk electrode using a Ag/AgCl microreference electrode. The conductivity of the electrolyte was 10 -6 S/m, but this value can be locally modified with time due to the corrosion reactions occurring on the disk electrode. Two consecutive measurements were performed and both gave the same results. The curves were also in good agreement with the potential distribution deduced from the theoretical calculations ͑z = 0.315 mm͒ ͑Fig. 3͒, showing that the hypothesis of the calculations was relevant. Thus, from these results, a three-step mechanism is proposed to explain the corrosion phenomena observed on the aluminum/copper model couple.

1.

Step 1: During the immersion of the disk electrode in the electrolyte, aluminum is the anode of the system and is in the passive state, while copper is cathodically polarized. The current distribution shows that a significant oxygen reduction occurs on the copper electrode close to the Al/Cu interface, whereas copper dissolution can be ignored during the first moments of immersion. Thus, due to the strong cathodic reduction, alkalinization of the electrolyte occurs locally.

2.

Step 2: When the pH reaches a value of 9, the aluminum depassivates and dissolves, leading to the formation of a crevice. The current distribution shows that the high cathodic current is restricted to the Al/Cu interface, which explains why the dissolution of aluminum only occurs at the Al/Cu interface, leading to the specific topography observed with SEM. 3.

Step 3: In the bottom of the crevice formed by aluminum dissolution, local variation of the chemical composition of the electrolyte can occur because this is an occluded zone. Then, local chemical conditions enable crevice corrosion of the copper. Finally, dissolved copper can form a deposit on the disk electrode. The place on which copper forms the deposit depends on the potential distribution. Additional experiments were performed with copper sulfate; the aluminum/copper model couple was immersed in sulfate solution, copper sulfate was added, and a copper ring rapidly formed on the disk electrode, showing that when copper ions are present in the electrolyte, they can easily form deposits on the disk electrode at the right potential.

Such phenomena are not dependent on the scale of the events ͑Fig. 6͒, and the three-step mechanism could explain the corrosion mechanisms in commercial alloys such as AA2024 and why the dissolution of copper-rich intermetallics occurs independently of the nature of the electrolyte. 34 Indeed, the results showed that the most significant feature was the potential distribution, which provokes the increase of oxygen reduction at the interface and thus induces a strong variation of the pH of the electrolyte. Our study shows that a simple model ͑pure aluminum/pure copper͒ is relevant for the study of what happens in copper-rich aluminum alloys and also underlines the specific and important role of the Al/Cu interface. This feature cannot be observed if the galvanic coupling is studied in the usual way with two independent electrodes with electrical contact in the same solution.

Conclusion

The results show that studying such a simple system as the pure aluminum/pure copper couple is relevant to understand the corrosion phenomena associated with copper-rich intermetallics in aluminum alloys. The potential and current distributions on the surface of the model couple were deduced from theoretical calculations and the potential distribution was checked using a Ag/AgCl microreference electrode. Increased cathodic activity related to oxygen reduction was observed at the Al/Cu interface from the beginning of immersion. Combination of the measurements and calculations with optical and SEM observations of the model couple after a longer immersion time ͑24 h͒ allowed a three-step mechanism to be proposed to explain the corrosion damage. In this mechanism, the formation of a copper deposit, always observed on commercial copper-rich aluminum alloys, does not appear to be the main step of the process but only a consequence of the corrosion damage. Our work underlines that the local alkalinization of the electrolyte at the Al/Cu interface plays a key role in the degradation of the system.

  where the constants k Cu and k O 2 are determined from experimental measurements of both the corrosion potential of the Cu electrode and the anodic and cathodic Tafel slopes of copper dissolution and oxygen reduction ͑b Cu , b O 2 ͒, respectively. The values determined experimentally were: E corr = -0.610 V/SSE, b Cu = 55 mV/dec, and b O 2 = 220 mV/dec, allowing calculation of the constants k Cu = 9.4 A cm -2 , and k O 2 = 7.6 10 -6 A cm -2 .

Figure 1 .

 1 Figure 1. ͑Color online͒ Schematic representation of the pure aluminum/ pure copper model couple.

Figure 2 .

 2 Figure 2. ͑Color online͒ Potential distributions on the surface of the Al/Cu model couple deduced from theoretical calculations for different values of electrolyte conductivity.

Figure 3 .

 3 Figure 3. ͑Color online͒ Potential distributions on the surface of the Al/Cu model couple deduced from theoretical calculations for two values of the probe position ͑electrolyte conductivity equal to 5 ϫ 10 -5 S m -1 ͒.

Figure 4 .

 4 Figure 4. Current distributions on the surface of the Al/Cu model couple deduced from theoretical calculations: ͑a͒ normal current and ͑b͒ radial current.

Figure 5 .

 5 Figure 5. ͑Color online͒ Optical observation of the Al/Cu interface of the disk electrode before immersion.

Figure 6 .

 6 Figure 6. ͑Color online͒ Optical observations after 24 h of immersion in a 10 -3 M Na 2 SO 4 solution of ͑a͒ the Al/Cu interface and ͑b͒ a copper particle embedded in pure aluminum on the disk electrode.

Figure 7 .

 7 Figure 7. SEM observations of the Al/Cu interface of the disk electrode after 24 h of immersion in a 10 -3 M Na 2 SO 4 solution: ͑a͒ and ͑b͒ show two different views of the crevice.

Figure 8 .

 8 Figure 8. ͑Color online͒ Observation of the disk electrode during immersion in a 10 -3 M Na 2 SO 4 solution. Bromothymol blue was added to the electrolyte.

Figure 9 .

 9 Figure 9. ͑Color online͒ Potential distributions measured in a 10 -3 M Na 2 SO 4 solution on the surface of the pure aluminum/pure copper model couple using a Ag/AgCl microreference electrode. Two consecutive measurements were performed.

Acknowledgments

The authors gratefully acknowledge Yannick Thebault for his valuable assistance in SEM-EDS.

Centre National de la Recherche Scientifique assisted in meeting the publication costs of this article.


Monitoring seasonal snow dynamics using ground based high resolution photography (Austre Lovénbreen, Svalbard, 79 o N)

Keywords: Ground based camera, snow cover dynamics, in-situ sensing Ny Ålesund Austre Lovénbreen

Arctic glaciers are reliable indicators of global climate changes. However, monitoring snow and ice dynamics in Arctic regions is challenging: some fast but key events can be missed since they are short in time but significant in the hydrological budget.

In the context of long term monitoring with high temporal and spatial resolutions of the snow cover dynamics, automated digital cameras were installed around the Austre Lovénbreen glacier basin (Spitsbergen, Norway, 79 o N). Despite data losses due to rough weather conditions and control electronics failure, a dataset of 2411 pictures (out of an expected 3294) was gathered over a one year hydrological period to assess the snow coverage of the glacier as a function of time with daily resolution. 73% of the total number of expected

Introduction

The results described in this paper are considered within the context of global climate change. Alpine-type polar glaciers are considered reliable indicators of global warming [START_REF] Moholdt | Recent elevation changes of svalbard glaciers derived from ICESat laser altimetry[END_REF][START_REF] Nuth | Glacier geometry and elevation changes on svalbard (1936-90) : a baseline dataset[END_REF][START_REF] Hagen | Geometry changes on Svalbard glaciers: mass-balance or dynamic response ?[END_REF][START_REF] Lefauconnier | Glaciers and climate in svalbard: Statistical analysis and reconstruction of the Brøggerbreen mass balance for the last 77 years[END_REF][START_REF] Dowdeswell | The mass balance of circum-arctic glaciers and recent climate change[END_REF]. Indeed, small polar hydro-systems exhibit significant responses to climatic anomalies. One original aspect of the present study is the very local scale of the analysis in order to understand in great details the phenomena under investigation. While classical analysis of global processes is usually performed on a wide glacier range [START_REF] Humlum | Mapping snow cover duration, avalanches and other geomorphic processes by automatic digital cameras[END_REF][START_REF] Farinotti | Snow accumulation distribution inferred from time-lapse photography and simple modelling[END_REF], such detailed analysis is more unusual to analyze a polar hydro-system. This originality is even more significant in an Arctic environment, where extreme weather conditions yield additional challenges for the continuous monitoring. Thus, the basic conclusion is that satellite imagery is insufficient to investigate the snow and ice dynamics of a glacier basin at daily scale, both due to cost and availability of materials as well as insufficient spatial and time resolutions.

Hence, in the framework of the 2006-2010 International Polar Year program Hydro-Sensor-FlOWS, a sensor network has been setup over the 10 km 2 large Austre Lovénbreen glacier basin (Fig. 1) at altitudes ranging from sea level to 450 m Above Sea Level. Dedicated automated digital cameras were installed at strategic spots of the glacier basin providing image coverage of the whole glacier (Fig. 2). While designed to provide continuous informations on the hydroglaciological dynamics, this setup also provides a continuous monitoring of the snow cover dynamics [START_REF] Hinkler | Detection of spatial, temporal, and spectral surface changes in the Ny-Ålesund area 79 o N, Svalbard, using a low cost multispectral camera in combination with spectroradiometer measurements[END_REF][START_REF] Newbery | An automated camera system for remote monitoring in polar environments[END_REF][START_REF] Buus-Hinkler | Snow-vegetation relations in a High Arctic ecosystem: Inter-annual variability inferred from new monitoring and modeling concepts[END_REF][START_REF] Corripio | Snow surface albedo estimation using terrestrial photography[END_REF][START_REF] Hinkler | Automatic snow cover monitoring at high temporal and spatial resolution, using images taken by a standard digital camera[END_REF][START_REF] Keller | Mountainview -precision image sensing on high-alpine locations[END_REF]. Unlike previous studies [START_REF] Hinkler | Detection of spatial, temporal, and spectral surface changes in the Ny-Ålesund area 79 o N, Svalbard, using a low cost multispectral camera in combination with spectroradiometer measurements[END_REF], this project focuses on the use of high resolution digital camera for snow cover monitoring, and hence requires more detailed signal processing than vegetation cover estimate. High resolution digital camera provide sub-meter resolution over the whole field of view, a result only accessible with sensor resolutions above the 10 Mpixel range [START_REF] Hinkler | Automatic snow cover monitoring at high temporal and spatial resolution, using images taken by a standard digital camera[END_REF]. Furthermore, unlike previous studies using oblique view [START_REF] Humlum | Mapping snow cover duration, avalanches and other geomorphic processes by automatic digital cameras[END_REF], geometric correction is used for a quantitative analysis of this dataset, with a high fraction of the area of the glacier covered by overlapping pictures taken by different cameras in a mosaique of projected pictures.

Our purpose here is to observe snow and ice melt dynamics and assess its contributions to the global budget of this characteristic small polar hydro-system. The signal processing sequence is based on qualitative datasets from which meaningful quantitative informations are extracted. By applying geometrical corrections to the ground-based oblique view pictures, mosaics of zenithal projected maps are obtained. After describing the image processing steps, we will exhibit the first results of the application of this dataset to snow cover monitoring: the assessment is performed on the 2009 hydrological season, whose short visible melt period extends from the beginning of July to the end of August.

Hydrological dynamics of a drainage basin are defined by several factors. We focus on the upstream aspect of the snow and ice contributions to these processes. Although satellite imagery is available, the phenomena under investigation require improved time and space resolutions. While informations gathered through remote sensing is rich by itself, it must be complemented with our ground-based pictures, better suited to our problem due to their improved time and spatial resolutions. Indeed, the available sensor network and the small basin area are consistent with this detailed monitoring, and justify the use of ground based high resolution pictures, assuming that the appropriate digital image processing is applied as will be described here.

Materials and methods

Instrument setup

Six automated instruments located around the Austre Lovénbreen glacier provide the data needed for this study. The instruments are based on a commercial off the shelf Leica D-Lux 3 or D-Lux 4 camera with 10 Mpixel resolution. These compact cameras were selected for their excellent lens, high resolution at the time the program was started (2007), and availability of a SD-card non-volatile storage medium (requiring only nine wires between the dry airtight enclosure in which the cameras are located and the user-accessible memory card). The wide angle lens, selected as the default setting at power up, is suitable for imaging the nearby small glacier basin: the furthest topic of interest is located 1,700 m away, as opposed to other similar references aimed at observing from more than 4 km away [START_REF] Farinotti | Snow accumulation distribution inferred from time-lapse photography and simple modelling[END_REF]. Arctic weather conditions, with wide temperature variations between shadow and sunlit conditions, or electrostatic buildup during snow storms, appear as significant electronics reliability limiting factors.

With this setup, oblique view pictures are gathered with high time and spatial resolution (the pixel size is at worst 0.46 m × 0.46 m considering the lens aperture), covering 96% of the glacier surface. Beyond the appropriate sensor resolution, the camera locations selected with safety and ground stability in mind provide a nearly full coverage of the glacier (Fig. 2).

We have selected, as a compromise between limited power consumption and improving the chances of appropriate weather and lighting conditions, the acquisition of 3 images every day at a fixed hour set to 7 h, 11 h and 15 h Coordinated Universal Time (UTC), corresponding to 8/12/16 h mean solar time at this longitude). Hence, weather permitting, a daily coverage of the glacier condition is guaranteed by this instrument network, providing a complement to space-borne satellite imagery whose availability is strongly dependent on weather conditions and satellite programming over the region of interest.

After selecting the 6,500 usable pictures for a given year, the resulting images must be converted from oblique view to zenithal projection to comply with the requirements of a quantitative analysis [START_REF] Laffly | High temporal resolution monitoring of snow cover using oblique view groundbased pictures[END_REF].

A detailed analysis of the success rate of this instrument network is provided in Fig. 3: recorded pictures out of an expected 3294: the solar pannels used to power the camera real time clock are only functional when the sun rises over the horizon, hence the poorer rate achieved when including the winter season. The winter season is anyhow characterized with dark pictures hardly usable for snow cover analysis, and is not considered during the image processing step.

out
A camera is considered successful if a picture has been grabbed during a given date, and if a correct date (embedded real time clock onboard the camera still powered) is associated to the EXIF header of these images. Considering that weather conditions might be the most significant factor affecting the proper working of the embedded electronics, the bottom chart (Fig. 3) displays the variation of the air temperature and precipitation in Ny-Ålesund during the same period: temperatures fluctuating around 0 o C, associated with condensation, appear as a significant factor of erratic behaviour of the electronics (June 2009 and October 2009). One parameter not displayed here is wind velocity: strong winds are associated with electrostatic charge buildup which might yield to microprocessor reset.

Digital image processing is performed in three steps. First, a geometric correction aims at projecting the images to a zenithal view and yield a mosaic of pictures covering the region of interest of the glacier. Secondly, with the objective of measuring the snow coverage in mind, specific image sets for which sudden snow cover evolution is observed are selected: these pictures are assembled as mosaics in order to compute the snow covered and iceexposed area fractions. Finally, the snow covered area delimitation on the mosaics with daily resolution over a whole melt-season is performed manually.

Geometric correction

The geometric correction strategy [START_REF] Laffly | High temporal resolution monitoring of snow cover using oblique view groundbased pictures[END_REF] is based on the use of geolocated markers positioned on the glacier surface. The GPS position of a hundred 2 m×2 m flags located on the whole glacier surface is measured (Fig. 4 a-c) in order to compensate for the lack of Ground Control Points (GCP), an issue often met when processing images in Arctic areas [START_REF] Korona | SPIRIT. SPOT 5 stereoscopic survey of polar ice: Reference images and topographies during the fourth International Polar Year (2007-2009)[END_REF][START_REF] Berthier | Recent rapid thinning of the "Mer de Glace" glacier derived from satellite optical images[END_REF].

The most challenging aspect of these processing steps is the identification of which flags are seen on each picture, requiring manual processing and a good knowledge of the field.

One reference picture from each oblique-view camera on which all flags are clearly visible is used for identifying the pixel position of each flag (Fig. 4 d). In addition, a satellite image including the geographic coordinates of each flag is mapped over the DEM in order to simulate similar viewing conditions under the ArcScene module of ArcGIS (Fig. 4 e).

The identification is achieved by comparing the oblique view picture with the 3D simulation.

These flags complement identifiable natural features on the slopes used as GCP in order to generate a robust geometric correction model based on Delaunay triangulation (Fig. 4 f).

The resulting geometrically corrected photo is mapped on the satellite image to assess the quality of the processing (Fig. 4 g).

Snow melt event selection

Following a classification step aimed at identifying the images with poor visibility (fog, clouds, water condensation on the camera lens, total darkness) in order to select the most interesting pictures, the photo dataset is available as contact prints on which significant snow melt events are readily visually identified. This selection step also allows for the removal of unusable pictures due to poor weather conditions, and identify missing pictures within the dataset due to electronics or battery failures. Nevertheless, unusable pictures in terms of snow cover analysis due to poor visibility provide valuable qualitative informations since they are often indicators of sudden weather events, either strong winds or precipitations.

Each melt event is characterized by a significant retreat of the ice-snow boundary as seen on the oblique view pictures. For a given day, depending on the availability of pictures from the cameras and on the lighting conditions, the best combination of images gathered during a single day, possibly at different times, is manually selected in order to improve the visibility of the phenomena under investigation. While some images include significant informations on the snow melt, quantitative processing is impossible due to water and frost condensation on the outer lens of the camera enclosure. Following this detailed analysis of the contact prints, the relevant dataset is selected for further processing.

Picture mosaics

For a given date selected because it exhibits a significant melt event, we wish to overlap some of the pictures representative of the phenomena under investigation in order to generate a zenithal projected mosaic compatible with a quantitative analysis and comparison with satellite images. The selection of the overlap order is performed for best visual results.

The mosaic quality is compared to reference satellite orthorectified images for assessing the positioning as well as the glacier boundary: once the mosaic is performed by overlapping the geometrically corrected images, a visual control assesses the result quality in order to avoid any chance of poor interpretation.

We have decided not to automate this mosaic generation step based on the Erdas software: although time consuming, the visual quality assessment is mandatory for the best quantitative exploitation. For instance, when unsuitable deformation is observed on the mosaics with respect to the satellite images, additional GCP can be added to compensate for the divergence of the geometric correction algorithm, as seen for example on bedieres flowing late during the melt season which are not considered as initial GCP. Using one image from each ground-based cameras, mosaics from 5 or 6 overlapping pictures are obtained (Fig. 5)

Results

In 2009, the melt season under investigation spans from the beginning of July to the end of August, during 43 days. From the first to the last visually identified melt events, one might theoretically look for 43 daily mosaics and 129 mosaics with 4 hour time resolution.

Considering the unusable or missing pictures, the available dataset allows for the generation of mosaics of the necessary picture subset to define the ice-snow boundary for 39 days.

This dataset allows for the regular monitoring with daily resolution of the snow cover. The quantitative analysis of the picture mosaic is performed by manually defining the snowice boundary, visible as a significant color change in the zenithal view. From the glacier boundary and this snow-ice limit, we compute a snow cover fraction of the glacier (Fig. 6).

Repeating this process for all significant events during the melt season, a snow cover fraction calendar is generated (Fig. 7), which provides the basic data for further snow cover dynamics observation.

A detailed analysis of the melt season shows that the phenomenon is discontinuous, with successive steps. Based on the observations on the raw photos gathered during the 43 days, only 12 events, during which the snow cover significantly varies, are observed (Fig. 7) and deserve the mosaic assembly.

During the 2009 melt season under consideration here, we 4 steps during which the glacier variation starts from a 100% snow cover condition (July 10 th ) to reach at the end of the melt season a 37% snow cover (September 12 th , Fig. 7). Following this date, the first snow fall again yields a 100% snow coverage. These steps are correlated with the location at which the process occurs: the first step is concerned with the glacier snout. The second and third phases are located where ice flows from the upstream cirques converge, close to the central area of the glacier: in this flattest region of the glacier, water storage yields slow snow cover variations with multiple small events. The last phase is characterized by a stabilization of the snow cover downstream of the accumulation cirques.

During these sudden melt events, the ice-snow boundary always brutally rises before stabilizing. The snow cover is stable between the melt events. High resolution images allow for the qualitative observation of rare events, such as the snow cover creeping downstream (Fig. 8). This phenomenon, recorded 3 times in 2009, occurs when the capillary cohesion of the snow cover reaches saturation, and is thus mainly associated with water concentration in the snow cover [START_REF] Dewalle | Principles of snow hydrology[END_REF].

Discussion

Image processing

From a technical point of view, we have considered the automated classification and generation of the mosaics as relevant due to the huge number of pictures [START_REF] Shukla | Estimation of debris cover and its temporal variation using optical satellite sensor data: a case study in Chenab basin, Himalaya[END_REF][START_REF] Racoviteanu | Optical remote sensing of glacier characteristics: A review with focus on the Himalaya[END_REF][START_REF] Casey | Characterization of glacier debris cover via in situ and optical remote sensing methods: a case study in the Khumbu Himalaya, Nepal[END_REF]. The results we obtained were not conclusive when attempting to classify the whole set of images from a given camera due to the large differences in snow coverage and exposed rocks and boulders between images obtained in the beginning and end of the season. Although automating the classification on a subset of images saves time by removing unusable pictures in this context and allowing the user to focus on the significant contact prints, the lengthy manual processing provides the opportunity of visually identifying interesting phenomena. Working on a dedicated region of interest rather than on the whole picture might improve the classification scheme by reducing the influence of weather conditions.

The automation of mosaic generation appears challenging due to the better capability of the human eye to select which picture to put on top of the stack in order to build the mosaic, and thus improve the resulting mosaic quality by enhancing the phenomena under investigation. Both shadows in good weather conditions and low cloud cover appear as possible factors for moving a picture towards the bottom of the stack while rising the contribution of the pictures least disturbed.

At the end, the lack of automation was also validated by the few significant events yielding a reasonable number of pictures to be processed: this conclusion would certainly differ if 43 mosaics had to be generated and analyzed.

Considering our initial goal of observing and quantifying the snow cover melt dynamics, we conclude that not all the images are needed at a given date. While obtaining a mosaic including contributions from all 6 cameras provide satisfactory technical and visual results at a given time, the quantitative analysis of ice-snow boundary identification and snow cover ratio over the glacier surface only requires 2 or 3 pictures under most conditions. In this particular example, once the snow cover has completely disappeared from the field of view of the downstream camera (camera # 2, Fig. 2), processing the images from this instrument is no longer mandatory. Similarly, identifying the ice-snow boundary on the downstream looking cameras means that the images from the upstream looking cameras are no longer needed to compute the snow cover fraction over the glacier surface.

Consistent with field trip observations, only a few melt events are relevant. Indeed, only 12 such events were identified during the 2009 melt season, justifying the removal of some of the melt models. As an example, linear ice-snow boundary motion as used in [START_REF] Li | Simulation of snow distribution and melt under cloudy conditions in an alpine watershed[END_REF] is inconsistent with the observations.

Snow cover

The melt process starts earlier than the first bare ice surface observation, while the icesnow boundary evolution provides information on the snow cover disappearance dynamics.

Camera photos only provide a binary information of ice or snow presence onto the glacier.

Since they do not provide any snow thickness information, the snow cover dynamics cannot be analyzed in terms of melt volume. The motion of the snow-ice boundary deduced from the picture set must be complemented with field measurements (snow thickness and density, air temperature) to quantify the melt rate on each point of the glacier. Indeed, these binary ice/snow glacier coverage maps provide the required data to feed a degree-day melting model in complement to the daily evolution of the glacier thermal state [START_REF] Braithwaite | Calculation of degree-days for glacier-climate research[END_REF][START_REF] Johannesson | Degree-day glacier mass balance modelling with applications to glaciers in iceland, norway and greenland[END_REF][START_REF] Zhang | Application of a degree-day model for the determination of contributions to glacier meltwater and runoff near keqicar baqi glacier, southwestern tien shan[END_REF]. In such a model, the melt coefficient is dependent on the nature of the glacier surface which includes such effects as albedo and heat capacity. Snow cover acts as a protective layer which delays ice melt [START_REF] Dumont | Retrieving surface glacier albedo using spatial and terrestrial remote sensing data. assimilation of albedo into snow model CROCUS to simulate distributed mass balance on Saint Sorlin glacier (French Alps)[END_REF].

Hence, a common melt process is applied to the whole glacier: water generated by the melt of the upper part of the snow cover accumulates in the snow until capillary cohesion ruptures and water flows at the snow-ice interface. Snow cover dynamics is thus interpreted in terms of glacier topography:

• the observation of a fast change of the ice-snow boundary on the snout is not only associated with a thinner snow cover and increased temperature but also a steeper slope inducing faster runoff,

• the central part of the glacier exhibits a similar process, but with a water-filled snow swamp: the six central steps observed in Fig. 7 exhibit a slow snow-ice boundary motion on this central area while we interpret this slow process as a consequence of the lack of slope which prevents the swamp from evacuating.

The stream network on the ice surface only organizes once all the snow has been removed: before this process takes place, water circulation is diffuse at the ice-snow interface and takes the snow away.

The last phase is characterized by a stabilization of the snow-ice boundary below the cirques, associated with a reduction of the water flow and snow melt. Since the snow melt rate is related to the snow stock, the thinner the snow cover, the lower the melt potential.

However, water circulation processes throughout the whole snow cover make this simple statement much more complex [START_REF] Dewalle | Principles of snow hydrology[END_REF].

Furthermore, we provide a calendar of the snow melt dynamics (Fig. 7) related in this paper to the weather chronicles acquired at the weather station of Ny-Ålesund and exhibiting significant melt events: these dates are either associated with unusually high temperatures and/or liquid precipitations.

1. The most intuitive cause of melt is temperature rise: the melt events dated July 14th and 22nd are solely associated with temperature increase while no precipitation occurs at this time, 2. some melt events follow significant rainfall, often associated with temperature increase:

the events dated July 11 th , August 11 th and August 23 rd follow rain by 1 to 2 days.

This delay indicates that melt is not associated with energy brought by the liquid water inducing snow fusion, but rather with capillary diffusion of the liquid water in the snow cover, yielding loss of cohesion.

3. some melt event are immediately associated with rain as well as temperature increase:

most significantly, the long rain period from July 28th to August 4th as well as August

Conclusion and perspectives

The novel approach of generating mosaics of geometrically corrected ground-based digital camera pictures provides the means for a precise monitoring of the snow cover dynamics over a small Arctic glacier basin. Its contribution to the hydrological processes is thus assessed.

The processing algorithm follows three main steps:

1. identifying and geolocating on one representative picture from each camera ground control points, followed by geometrical corrections of the image, 2. assembling mosaics of the relevant events during one snow melt season, 3. manually drawing the ice-snow boundary over the glacier surface, and computing the relative area of the snow-covered fraction of the glacier.

We observe, from the resulting dataset covering one melt season, the discontinuous melt rate of the snow cover. A few significant snow melt events are observed, associated either with temperature increase or liquid water diffusion in the snow. Liquid water, either provided as rain or surface runoff on the glacier, appears as a significant accelerating factor defining melt rate beyond the influence of air temperature.

Further study aim at exploiting these binary ice/snow cover maps in a degree-day melting model, crossing the information of snow-ice presence with daily thermal state of the glacier.

From these computations, we will provide the potential daily water volume produced by melt, and compare these with hydrological measurements performed downstream in the moraine.

The expected result aims at a better understanding of the hydrological processes in a polar glacier region, and most significantly the snow or ice contribution to the runoff during the melt season.

Figure 1 :

 1 Figure 1: Located 79 o N, the 10 km 2 Austre Lovénbreen glacier basin is the subject of these investigations.

Figure 2 :

 2 Figure 2: The network of automated digital cameras observing the Austre Lovénbreen glacier basin. Bottomright: a double enclosure keeps the digital camera in a dry atmosphere while a dedicated low-power microcontroller based electronics controls the power to the camera at preset time intervals.

Figure 3 :

 3 Figure 3: Chart sumarizing the success rate of the cameras located around the Austre Lovénbreen basin during the 2009 hydrological season, and associated weather conditions.

Figure 4 :

 4 Figure 4: Signal processing steps requiring the preliminary measurement of the position of ground control points over the glacier surface. A significant issue is the interpretation of the pictures in order to locate these features and identify which flag is located at known positions.

Figure 5 :Figure 6 :

 56 photo from cam 3

Figure 7 :

 7 Figure7: Top: snow melt over one season is separated in four main phases, with a snow cover evolution from 100% to the final value of 37% over the glacier area. Bottom: the generation of mosaics allows for a calendar of the melt rate in which the snow cover variation is compared to weather conditions (precipitations and temperatures in Ny-Ålesund, 6 km from the Austre Lovénbreen; data provided by Meteo Norsk).

Figure 8 :

 8 Figure 8: Creeping of the snow cover has been recorded at the snow-ice interface three times in 2009: here a zoom on the set of photos from camera 3 captured July 29 th .

7

  th exhibit multiple snow cover retreat, 4. finally, a single event occurring August 16 th remains unexplained since it is neither associated with precipitation nor temperature rise. Mechanical dynamics of the snowcover could explain it by drainage of the top snow swamps and by a break in the water-supersaturated snow.

  of an expected 6570 pictures to be gathered during one year (April 1st 2009-March 30th

2010) by the 6 cameras, we have recorded 3921 pictures, or 60%, including the full date and time informations in the EXIF header. When focusing on the hydrological year spanning from April 1st 2009 to September 30th 2009, the success rate increases to 73%, with 2411

Acknowledgements

Funding for this research program was provided by the French National Research Agency (ANR) Hydro-Sensors-FlOWS program, IPY #16 program, and IPEV.


J Dudouet

D Juliani

M Labalme

J C Angélique

B Braunn

J Colin

Guérin

D Cussol

Ch Finck

J M Fontbonne

H Guérin

P Henriquet

J Krimmer

M Rousseau

M G Saint-Laurent

Comparison of two analysis methods for nuclear reaction measurements of 12 C + 12 C interactions at 95 MeV/u for hadrontherapy

Keywords: fragmentation, ∆E-E, cross section, KaliVeda PACS: 25.70.Mn, 29.85.Fj

published or not. The documents may come

Comparison of two analysis methods for nuclear reaction measurements of 12C +12C interactions at 95 MeV/u for hadrontherapy

Introduction

Carbon ions with energies ranging from 80 to 400 MeV/u are used in hadrontherapy to treat cancerous tumors. Their advantage compared to photons lies in the great accuracy of ion beams to target the tumor while sparing the surrounding healthy tissues (due to the high dose deposition in the Bragg peak and the small angular scattering of ions). The physical dose deposition is affected by the fragmentation of the ions along their penetration path in the patient tissues [1]. This process leads to a consumption of the primary beam (up to 70% for 400 MeV/u 12 C in water) and to a delocalisation of the dose in the healthy tissues [2].

Simulation codes are used to compute the transport of ions in matter, but the constraints on nuclear models and fragmentation cross sections in the energy range used in hadrontherapy (80 to 400 MeV/u) are not yet sufficient to reproduce the fragmentation processes with the required accuracy for clinical treatments. To improve the knowledge on the 12 C fragmentation process, experiments have been performed in Japan and in Europe for more than 15 years. Measurements of light charged fragment production in water and PMMA (C 5 H 8 O 2 ) have been started by the Japanese treatment centers (Chiba and Hyogo) and by the Gesellschaft für SchwerIonenforschung (GSI) biophysics department in the energy range 200-400 MeV/u [3] [4] [5] [6] [7] [8].

To extend these data to the lowest energies, a first integral experiment on thick water equivalent targets has been performed by our collaboration in May 2008 at the Grand Accélérateur National d'Ions Lourds (GANIL, France). Energy and angular distributions of the fragments produced by nuclear reactions of 95 MeV/u 12 C with thick PMMA targets have been obtained [9]. Comparisons between experimental data and simulations, using different physics processes of GEANT4 [10] for modeling the passage of particles through matter, have shown discrepancies. Nevertheless, it is difficult to constrain the nuclear reaction models by a direct comparison to thick targets experiments.

To improve the models and reach the accuracy required for a reference simulation code for hadrontherapy, a second experiment has been performed on thin targets on May 2011 at GANIL to study C-C, C-H, C-O, C-Al and C-nat Ti reactions at 95 MeV/u. The particle identification has been performed by using telescopes (stack of detectors) and the ∆E-E technique. The experimental set-up will be presented in the next section.

In such experiments, the usual procedure to extract the particle natures and the particle energies from raw data consists in two separate steps: the particle identificaton and the energy calibration of the detectors. The particle identification is done by using graphical cuts on so called ∆E-E maps. The energy calibration is done by using well defined reference energies (elastic energies, punch through energies). We propose in this paper an alternative procedure which performs the particle identification and the detector calibrations in one single step. The two different analysis methods will be described, the first one called "graphical method" uses graphical cuts, the second one called "KaliVeda method" uses the KaliVeda framework developed by the INDRA collaboration [11][12]. The main advantage of the KaliVeda method is that once the identification grids are created, they allow to do both the detector calibration and the particle identification at the same time thanks to automated procedures which saves time in data analysis. A comparison of the results obtained with the two methods for the C-C reaction will be presented in the last section.

Experimental setup

Charged particles detection

The experiment has been performed using the ECLAN reaction chamber at the GANIL G22 beam line. The energy of the impinging carbon ions was 95 MeV/u on different thin targets (C, CH 2 , Al, Al 2 O 3 , nat Ti) with area densities of about ∼100 mg/cm 2 . In this section we will report only about carbon target data. The setup consists of four ∆E thin -∆E thick -E telescopes mounted two by two on two rotating stages that allows rotation inside the chamber from 7 • to 43 • at a fixed distance of 22 cm behind the target. Thus the solid angle covered by each telescope was 6.5 msr. A fifth telescope was at fixed angle at 4 • , downstream, covering a solid angle of 0.4 msr as shown in Fig. 1.

The five ∆E thin -∆E thick -E telescopes consist in two silicon detectors followed by a 10 cm thick CsI scintillator. The thicknesses of the silicon detectors were 150 µm and 1000 µm respectively.

Beam Monitor

The beam monitoring is based on the measurement of fluorescence Xrays emitted by a thin Ag foil (7 µm thick) set in the beam after the reaction chamber. X-rays are detected by means of two Si(Li) detectors located at 90 • with respect to the beam direction. These two detectors are calibrated at different beam intensities ranging from 10 3 to 10 6 ions/s by means of a plastic scintillator intercepting the beam. The beam intensity is determined with the two detectors (Si(Li)). The detected fluorescence X-rays, emitted by the thin Ag foil, are represented by two peaks in the energy spectra of the Si(Li) (Fig. 2(a)). These peaks have been integrated while taking into account the background. Different background subtraction methods have been tested: no background, constant background, linear background, polynomial background. They are all compatible within 2%. The energy spectrum of the plastic scintillator presents several peaks corresponding to the number of carbon ions detected in coincidence (Fig. 2(b)). These peaks have been fitted with gaussian functions and integrated hereafter.

Data Analysis using graphical cuts

Particle identification

The ∆E-E technique is used to perform the identification of the charged particles. This method is based on the fact that the energy loss of charged particles in matter depends on their charge and mass. Figs. 4 and5 show examples of identification maps for the 250 µm thin Carbon target. Each hyperbole corresponds to a particle characterized by its charge and mass numbers (Z,A). Figure 4 represents the energy deposited in the thin silicon detector relative to the energy deposited in the thick one. The hyperbolas back bending are due to the particles which cross the second stage of the telescope. Figure 5 shows the energy deposited in the thick Si detector relative to the energy deposited in the CsI scintillator. This map is complementary to the previous one, thanks to the identification of the particles passing through the first two stages. The identification is performed by a set of contours for each telescope in order to separate each particle in charge and mass, from Z = 1 to Z = 6. These contours are represented by black lines on the maps of the figures 4 and 5. Since the data acquisition is triggered when the second stage is hit, the identification threshold at low energy is determined by the thickness of the first stage. The minimum energy requisite to cross the first stage (150 µm) and thus to hit the second one depends on the nuclear charge of the particle, these threshold values are given in Table 1. The CsI crystal thickness prevents most of the charged particles from passing through the scintillator. Indeed, the lightest particles namely protons, deuterons and tritons, need respectively energies greater than 198, 130 and 103 MeV/u to cross the CsI crystal.

Energy calibration

The energy calibration of the two silicon detectors is based on the maximal deposited energy in the second stage for a given particle without crossing it. Thanks to the known thicknesses of the Si detectors, this energy threshold is computed, for each particle, by using SRIM [13] simulations and is linked to the corresponding channel number of the thick Si detector. The energy loss in the thin Si detector is deduced from the energy deposited in the thick Si detector. This allows to establish the correspondence between the energy deposited in the thin Si detector and the channel number in the corresponding charge to digital converter (QDC) module. Fig. 6 shows an example of the calibration points obtained for one thick Si owing to this method. As the response of a silicon detector does not depend on the nature of the particle, the correlation between the energy loss and the channel number is the same for all particles. The calibration parameters are obtained by a linear fit of the points as shown on Fig. 6.

Contrary to silicon detectors, the produced light in the scintillator depends on the nature of the particle and on its energy. It is necessary to use a calibration function with a distinct set of parameters for each isotope. The chosen function, based on the one proposed by Stracener et al [14], is described by the following formula:

E MeV = α + βE channel + γln(1 + δE channel ) + ηE 2 channel (1)
This equation takes into account the quenching of the scintillator at small energy as well as the linear response at higher energy. The term, ηE 2 channel , has been added to take into account a possible saturation at very high energy. The energy loss in the CsI detector is deduced from the energy deposited in the thick Si detector using SRIM simulations. This method propagates the errors made on the energy calibration of the thick Si detector to the energy calibration of the CsI detector. The lower the deposited energy in the thick Silicon detector is, the larger the error on the residual energy released in the CsI is. With our experimental setup, this error has been estimated, from almost 30% for 95 MeV protons to 4% for 95 MeV/u carbon ions. An example of the calibration functions obtained is shown on Fig. 7. The five parameters (α, β, γ, δ and η) are determined for each isotope and for each telescope.

Data Analysis using KaliVeda

An new analysis method based on the KaliVeda toolkit [12] has been used for the analysis of the experiment. KaliVeda is an object-oriented software, based on the ROOT [15] framework for data processing, with many useful classes for data analysis and energy loss calculation. It has been developed to analyse data obtained by the ∆E-E technique. It was firstly developed for the analysis of the INDRA multi-detector.

KaliVeda tools

In this part, the KaliVeda tools which are used for the analysis will be presented.

The Energy Loss functional

KaliVeda is using a functional describing accurately the energy loss ∆E in a detector as a function of the residual energy E deposited in a second detector in which the particle has stopped. This functional is the one proposed by L.Tassan-Got in which the charges Z and the masses A of particles are parameters [16]. For detectors delivering a linear response (as silicon detectors), it reads:

∆E = [(gE) µ+ν+1 + (λZ α A β ) µ+ν+1 + ξZ 2 A µ (gE) ν ] 1 µ+ν+1 -gE (2)
where g, µ, ν, λ, α, β, ξ are fit parameters.

For detector delivering a non linear response versus the deposited energy this functional needs to be corrected. For scintillators, it allows to take into account the quenching effect of scintillation leading to a non linearity of the emitted light versus the deposited energy. The energy E in equation 2 must now be expressed as a function of the light h emitted by the scintillator. The light response of CsI(Tl) crystals is deduced from the Birks formula [START_REF] Birks | Theory and Practice of Scintillation Counting[END_REF] and allows to express the energy released in the CsI as a function of the emitted light :

E = h 2 + 2ρh 1 + ln 1 + h ρ (3)
where ρ = ηZ 2 A is a new fit parameter.

The KaliVeda identification grid

As with the previous analysis method, the particle identification is performed using the ∆E-E method. In such a representation, the particles are grouped around quasi-hyperbolas characterized by their charge and mass. Then, for each detected particle, the closest line has to be identified and the corresponding charge and mass have to be assigned to the particle. It can be done by using graphical contours as presented in the previous section or by using the KaliVeda identification tools. KaliVeda is using identification grids which are graphs that contain a set of identification lines as described by equation 2 (cf Fig. 8 and Fig. 9). Each line corresponds to an isotope and reproduces the shape of its quasi-hyperbola on the ∆E-E map. In a first step, to construct the grids, the user has to do a "handmade identification grid" for a few lines on the ∆E-E map. To do so, one has to give the coordinates of about ten points of the relevant lines and to specify charge and mass of the corresponding particle as shown on Fig. 8. Only few lines are needed here (those with the highest statistics), whereas the other one will be computed afterwards.

In a setup with three stages telescopes, two identification grids are needed per telescope (one for the particles stopped in the thick Silicon detector and one for those stopped in the CsI). Fig. 8 represents the "handmade identification grid" for the particles stopped in the thick Silicon detector (superimposed to the ∆E thin -E thick map).

In a second step, the functional of equation 2 is fitted to the handmade ∆E thin -E thick identification grids. The same method is used for the ∆E thick -E CsI identification grids but the energy deposited in the CsI is here expressed as a function of the emitted light as described in equation 3. The fit result allows to obtain the analytical expression of the identification lines for all isotopes. New identification grids including the lines of isotopes with low statistics can now be created, as shown on Fig 9.

By this identification method, two sets of accurate ∆E-E identification grids (∆E thin -E thick and ∆E thick -E CsI ) have been determined.

The KaliVeda method

In this part, the method we have developed will be presented. This method is using the KaliVeda tools described above.

Particle identification

To identify the detected particles, the procedure works as follows:

• For each experimental measurement, the associated identification grid is loaded and the algorithm computes the distance on the ∆E-E map between the coordinates of the detected particle ((∆E thin ;E thick ) or (∆E thick ;E CsI )) and the different isotope's lines.

• The charge and mass assigned to the particle are the one corresponding to the closest line of the identification grid.

• A new variable named A real is determined as a function of the distance to the isotope's line with A real ∈ [A-0.5;A+0.5]. The value of A real is close to the mass number A when the experimental point is close to the identification line. The minimal or maximal values of A real correspond to points located at the same distance of two lines with consecutive mass numbers.

• This new variable is used to determine the Particle Identifier (PID) expressed as PID = Z + 0.1

×(A real -2×Z).
As examples, the PID for α particles may range from 1.95 to 2.05, the PID of 3 He isotopes may ranges from 1.85 to 1.95 and the PID for 6 He may range from 2.15 to 2.25. Fig. 10 illustrates an example of PID calculation for alpha particles. The result of the PID calculation for all the detected isotopes is shown on Fig. 11. The quality of the identification is very good and the different isotopes are clearly separated, even those with low statistics.

Energy calibration

Once the particle identification has been achieved, the energy calibration of the detectors can be performed.

For energy calibration of the Silicon detectors, the analytical expressions (eq. 2) of the identification grids for the ∆E thin -E thick map are used.

First, a theoretical identification grid is generated. The points of this grid are obtained with the energy loss calculation tool of KaliVeda, taking into account the Silicon detectors thicknesses. The energy loss functional is then fitted to the theoretical identification grid. It gives an analytical expression of the theoretical lines. Then, the analytical expression of the experimental lines are obtained as described in part 4.1.2.

Two sets of parameters have thus been obtained, one for the theoretical grid expressed in MeV, and one for the experimental one expressed in channels. As Silicon detectors are delivering a linear response, their energy calibration consists in determining the affine relation between channels and MeV. Four calibration parameters have to be determined, two for the thin Silicon detector and two for the thick Silicon detector. By applying the correct transformations respectively to ∆E and E in the experimental energy loss functional, the two functionals (experimental and theoretical) will be equal. This procedure is done by using an iterative algorithm which adjusts the four calibration parameters to make the two functionals as close as possible. Fig. 12 represents a ∆E-E calibrated map on which the theoretical grid is superimposed. Both are in perfect agreement. The CsI energy calibration can not be done with the same method because of the non linear response of the emitted light versus the deposited energy (cf eq. 3). The thick Silicon detector energy calibration is used to determine the deposited energy in the CsI as follows. Once a particle is identified in the ∆E thick -E CsI map, its charge, mass and energy loss in the thick Silicon detector are known. Then, a KaliVeda tool computes the corresponding incident energy of the particle thanks to the known energy loss in the Silicon detector. The residual energy which will be deposited in the CsI can thus be deduced. The disadvantage of such a method, as the one used in section 3.2, is that the error made on the thick Silicon detector calibration parameters is propagated to the energy deposited in the CsI. This error has been estimated, in this case, from around 20% for 95 MeV protons to 4% for 95 MeV/u 12 C.

The main advantage of using the KaliVeda toolkit is that once the identification grids are created, they allow the detector calibration and the particle identification at the same time. With such a method, the energy calibration is done not only on a few points but for all the points of the identification lines. Moreover, these procedures are automated which saves time in data analysis. The time required, for the user, to construct the necessary elements for the calibration and identification is about one hour per telescope.

Comparison of the two analysis methods

The obtained data have been analyzed with the two methods described above. In this part, some results of fragmentation cross section measurements for the carbon target will be presented and the two analysis methods will be compared.

Firstly, the angular distributions, obtained for angles from 4 • to 43 • , will be compared. The 12 C fragmentation cross sections for a A Z X fragment are obtained as follows :

dσ dΩ ( A Z X) = NA Z X × A target N12 C × Ω × ρ × th × N A (4)
where NA Z X is the number of A Z X fragments detected, A target is the target mass, N12 C is the number of incident carbon nuclei, Ω is the solid angle of the detector, ρ and th are the target density and thickness respectively and N A is the Avogadro's number. Fig. 13, Fig. 14 and Fig. 15 display the angular distributions obtained with the two analysis methods for different emitted fragments. Statistical errors only are represented, systematic errors are still under study. As expected, the fragment emission is dominated by H and He isotopes with a predominance of α at low angles (below ∼ 10 • ) which is compatible with the alpha cluster structure of the 12 C. The results also show an angular emission most forward peaked for heavier fragments. The two analysis methods give consistent results. The discrepancies between the two methods are of ± 5% for most of the isotopes and angles. However, the discrepancies can reach 15% at low angles for light isotopes or for isotopes with low statistics. These differences are due to the widths of the particle selection cuts which are not the same for the two methods (width of the contours for the first method and width of the identification lines for the second one). The differences are more important for protons because of the increasing "background" on the ∆E-E maps at low energy. As this "background" is mainly due to protons escaping the CsI by scattering, a large identification line width has been used in the KaliVeda method in order to not lose proton statistics.

Secondly, the energy distributions will be compared. Fig. 16 represents the energy distributions of protons and α fragments at different angles obtained with the two analysis methods. These distributions are dominated at low angles by a peak centered close to the beam energy (95 MeV/u). The energy and amplitude of this peak are decreasing with the emission angle. This observation confirms that most of the detected fragments are coming from the projectile fragmentation. The energy distributions are also in good agreement. The discrepancies between the two methods are due on one hand to the propagation of the energy calibration of the thick Silicon detector. This effect can be seen on Fig. 16(b) in which the differences between the two analysis method are increasing with the energy. It should be noticed that the discrepancies stay under 10% even for high energy protons. On the other hand, the discrepancies are also coming from the particle selection width as explained for the angular distribution. This effect is clearly visible for the higher energies of Fig. 16(a) in which the largest width of the KaliVeda methods leads to the selection of background (due to the presence of 2 H and 3 H). As the background of the experiment is low, this effect is statistically low even with the choice of a selection width larger in the KaliVeda method than in the classical one. The comparisons between the energy distributions obtained with the two analysis methods show the coherence and the robustness of the two identification and calibration methods.

Conclusion

Two different methods for the data analysis of impinging 12 C ions at 95 MeV/u on thin carbon target have been presented in this paper. Measurements for twenty different angles led to production rates, angular distributions and energy spectra for fragmented particles, ranging from protons to carbon ions. Two different analysis methods have been used, the first one is based on graphical cuts while the second one uses the KaliVeda framework. Both of them are in good agreement, within a discrepancy close to 5% for most angular distributions, the discrepancy is higher concerning energy distributions in some specific cases. Angular distributions are dominated by the emission of light fragments (Z < 3), more especially by α particles at forward angles, that is consistent with the α cluster structure of the 12 C ion. Energy distributions dominated at low angles by a peak close to the beam energy indicate that most of the emitted particles result from the projectile fragmentation during nuclear reactions. As mentioned above, other targets have been used in this experiment using the same beam at the same energy, which will be the subject of an upcoming paper. The analysis of these data will help to determine the double differential cross section for the different couples C-H, C-O, C-Al and C-Ti.

Figure 1 :

 1 Figure 1: Representation of the experimental setup.

Figure 2 :

 2 Figure 2: Left-hand side: Energy spectrum of the X-rays emitted by the Ag foil obtained with a Si(Li) detector. Right-hand side: Energy spectrum obtained with the plastic scintillator for the beam monitor calibration.

Figure 3 :

 3 Figure 3: Intensity calibration of the SiLi detector thanks to the plastic scintillator intercepting the beam at low energy (I plast = p0×I Si(Li) + p1).

Figure 4 :Figure 5 :

 45 Figure 4: Identification map between the first stage and the second stage of a telescope. The black lines represent the different graphical contours used for the particles identification.

1 Figure 6 : 1 Figure 7 :

 1617 Figure 6: Energy calibration for the thick silicon detector of a telescope (linear fit).

Figure 8 :

 8 Figure 8: Handmade identification grid built on a ∆E-E map for the most produced isotopes.

Figure 9 :

 9 Figure 9: New identification grid obtained by fitting the Tassan-Got functional to the handmade identification grid.

12 Figure 11 :

 1211 Figure 11: PID calculation for all detected isotopes.

Figure 12 :

 12 Figure 12: Result of the Silicon detectors calibration. The theoretical grid is superimposed to the ∆E-E calibrated map (the right-hand side is a zoom of the red square). Only statistical errors are represented.

3 Figure 13 :

 313 Figure 13: Comparisons of the Hydrogen angular distributions obtained with the two analysis methods. Only statistical errors are represented.

4 Figure 14 :Figure 15 :

 41415 Figure 14: Comparisons of the Helium angular distributions obtained with the two analysis methods. Only statistical errors are represented.

Figure 16 :

 16 Figure 16: Comparisons of the energy distributions obtained with the two analysis methods for protons and α at different angles.

Table 1 :

 1 Minimum energy for detection as a function of the charge of the particles.

	Nuclear charge Threshold (MeV)
	1	4
	2	14
	3	30
	4	46
	5	61
	6	81

  Example of PID calculation for 4 He fragments.

	2000 Thin silicon (Channel) 3000	200 300 400 500 Counts 600
		100
		0
		1.96	1.98	2	2.02	2.04
				PID
	1000
	1000	2000		Thick silicon (Channel) 3000 4000
	Figure 10: