Frances Laëtitia

Grivet Manuel

Renault Jean-Philippe
email: jean-philippe.renault@cea.fr

Ducret Didier
email: didier.ducret@cea.fr

L Frances
email: lfrances@univ-fcomte.fr

Hydrogen Radiolytic Release from Zeolite 4A / Water Systems Under γ Irradiations

Keywords:

Although the radiolysis of bulk water is well known, some questions remain in the case of adsorbed or confined water, especially in the case of zeolites 4A, which are used to store tritiated water. An enhancement of the production of hydrogen is described in the literature for higher porous structures, but the phenomenon stays unexplained. We have studied the radiolysis of zeolites 4A containing different quantities of water under 137 Cs gamma radiation. We focused on the influence of the water loading ratio. The enhancement of hydrogen production compared with bulk water radiolysis has 1 Abbreviations : WLR : Water Loading Ratio ZWS : Zeolite 4A / Water System been attributed to the energy transfer from the zeolite to the water, and to the influence of the water structure organization in the zeolite. Both were observed separately, with a maximum efficiency for energy transfer at a loading ratio of about 13%, and a maximum impact of structuration of water at a loading ratio of about 4%.

Introduction:

Zeolites are widely used nanoporous materials. They are encountered in a broad field of applications, such as catalysis, filtration and isotopic separation [START_REF] Kotoh | Multi-component behaviors of hydrogen isotopes adsorbed on synthetic zeolites 4A and 5A at 77.4K and 87[END_REF][START_REF] Kotoh | Molecular-sieving effect of zeolite 3A in adsorption of H2, HD, D2[END_REF][START_REF] Montanari | On the mechanism of adsorption and separation of CO2 on LTA zeolites: An IR investigation[END_REF][START_REF] Zhu | Water vapour separation from permanent gases by a zeolite-4A membrane[END_REF]. Zeolites A, which are hydrophilic, are also used to store tritiated water, generated and required by nuclear applications like ITER. This experimental reactor is supposed to produce, amongst other tritiated wastes, high quantities of pure tritiated water during its working period [START_REF] Pamela | ITER tritiated waste management by the Host state and first lessons learned for fusion development[END_REF]. Its development requires safe storage solutions.

Synthetic Zeolites A are built with corner sharing TO4 tetrahedrons, where T corresponds to silicon or aluminum [START_REF] Breck | Critical Review of Rate Constants for Reactions of Hydrated Electrons Hydrogen Atoms and Hydroxyl Radicals (OH/O -) in Aqueous Solution[END_REF]. At the nanometer scale, their crystalline structure is composed of two kinds of cages. Sodalite cages, also called β cages, are a truncated octahedral shape, with a diameter of 6.6 Å, disposed at each corner of a cubic arrangement. The second kind of cages, called supercages or α cages, are located in the center of this array, and are characterized by a diameter of 11.4 Å. Sodalites cages are linked together by prisms, called double-four-membered rings (D4R).

The framework negative charge (induced by the aluminum valence) is compensated by cations, included in the structure. The nature and location of those cations affects the properties of the zeolite.

In zeolites 4A, the positive charge is provided by monovalent sodium cations, leading to the following chemical composition for a unit cell: [Na12(SiO4)12(AlO4)12]. The zeolite 4A is strongly hydrophilic. A unit cell can adsorb 27 molecules of water. Among those 27 molecules, some are commonly localized as molecules adsorbed in the crystalline structure: 4 or 5 in each sodalite cage and 20 to 23 in each supercage [START_REF] Crupi | T dependence of vibrational dynamics of water in ion-exchanged zeolites A: A detailed Fourier transform infrared attenuated total reflection study[END_REF][START_REF] Demontis | Dynamical Properties of Confined Water Nanoclusters: Simulation Study of Hydrated Zeolite NaA: Structural and Vibrational Properties[END_REF].

The quantity of water is expressed as Water Loading Ratio (WLR) given in percentage, equal to the mass of water divided by the mass of zeolites. The maximum loading ratio for zeolites 4A depends on the synthesis method and the binder quantity, but is about 20 %. It has been shown, considering the adsorption heat [START_REF] Moïse | Adsorption of water on X and Y zeolites exchanged with barium[END_REF] or FTIR results [START_REF] Crupi | A FT-IR absorption analysis of vibrational properties of water encaged in NaA zeolites: evidence of a structure maker role of zeolitic surface[END_REF], that the water adsorption is heterogeneous, and occurs in three-steps, leading to different kinds of confined water. The first adsorbed water molecules solvate the cations which then migrate, opening the access to the β cages.

The α cages are then filled.

Even if bulk water radiolysis is well known, even in the case of self-radiolysis of tritiated water [START_REF] Breck | Critical Review of Rate Constants for Reactions of Hydrated Electrons Hydrogen Atoms and Hydroxyl Radicals (OH/O -) in Aqueous Solution[END_REF][START_REF] Ershov | A model for radiolysis of water and aqueous solutions of H2, H2O2 and O2[END_REF], [START_REF] Stolz | Self radiolysis of tritiated water[END_REF], many questions remain about the radiolysis of adsorbed and confined water. An enhancement of hydrogen production, in some porous materials has been observed, for different irradiation conditions [START_REF] Cecal | Hydrogen yield from water radiolysis in the presence of zeolites[END_REF][START_REF] Le Caër | Water Radiolysis: Influence of Oxide Surfaces on H2 Production under Ionizing Radiation[END_REF][START_REF] Le Caër | Radiolysis of Confined Water: Hydrogen Production at a High Dose Rate[END_REF][START_REF] Rotureau | Radiolysis of Confined Water: Molecular Hydrogen Formation[END_REF]. The same phenomenon is observed in the precise case of zeolites A [START_REF] Nakashima | Radiolytic Gas Production from Tritiated Water Adsorbed on Molecular Sieve 5A[END_REF][START_REF] Nakashima | Radiolytic Gas Production from Tritiated Water Adsorbed on Molecular Sieve 5A[END_REF], or structurallyclosed zeolites: faujasites [START_REF] Nakashima | Radiolytic Gas Production from Tritiated Water Adsorbed on Molecular Sieve 5A[END_REF], but the literature is quite limited in this area, even if the impact of confinement on radiolysis has been shown [START_REF] Foley | Radiolysis of Confined Water: Production and Reactivity of Hydroxyl Radicals[END_REF].

Since the interaction and behavior of water strongly depends on the loading ratio of zeolites, we followed gas production depending on the level of water filling of the zeolites. We chose gamma irradiations, which are close to β -irradiations induced by tritiated water in terms of the dose rate. We took special care with respect to the dosimetry: the primary energy deposition in water is commonly estimated to be equal to the total energy deposited multiplied by its mass ratio [START_REF] Laverne | H2 Production in the Radiolysis of Water on CeO2 and ZrO2[END_REF], [START_REF] Laverne | H2 production in the Radiolysis of Aqueous SiO2 Suspensions and Slurries[END_REF]. We proposed an improvement of this consideration using Monte-Carlo simulation. We focused on the dihydrogen production since it may create, with dioxygen, an explosive atmosphere over a large range of concentration. We adopted two different approaches, regarding both primary energy deposition in the water and energy deposition in the whole system containing zeolites 4A and water.

Materials and Methods

2.1.Sample conditioning

Synthetic hydrophilic zeolites 4A were supplied by Molsiv Adsorbents. Their maximum water capacity, expressed as the ratio of water mass under dried zeolite mass, given by the supplier is 19.6 %. To improve the mechanical properties of zeolites 4A, a chemical inert non-porous binder, presenting the same atomic composition as the zeolites, is included during the synthesis. We estimated the proportion of binder by comparing the experimental water quantity adsorbed at saturation for our samples and that for a pure 4A sample, provided by CECA. The water saturation obtained by exposure to ambient water vapor, is investigated by thermogravimetry and is obtained as a mass percentage of 25 % for pure zeolites, against 19 % for our storage samples. The binder proportion is so estimated to be 24 %. Nevertheless, as the atomic composition of binder is close to the one of pure zeolites, no significant influence on interaction sections, and thus on irradiation effects, is expected. 4 g of zeolites 4A were introduced in 10 cm 3 air tight glass ampoules. Samples were then conditioned with a unique bench, comprising a pump, a pure water tank, a manometer and a junction to connect ampoules containing samples (figure 1). The use of a unique bench and of a unique ampoule for the different steps of samples conditioning prevents from the exposition of zeolites 4A to atmospheric water. As they are hydrophilic, such an exposition would have led to residual water adsorption by zeolites. Each ampoule was degassed under vacuum at 623 K, during at least 12 h, to ensure the elimination of water from zeolites initially equilibrated with water vapor from the atmosphere. The complete dehydration has been checked by weighing, after isolation and removal of ampoules from the bench.

Adsorption of controlled quantity of water, from about 3 % to saturation, was realized by manometry. Zeolites were exposed to the chosen water vapor pressure, released in the line whose volume has been measured to 2.58 dm 3 . The quantity of water adsorbed was thus pre-determined using the perfect gas equation, and checked by the mass increase. Details of the quantities of water adsorbed are given in Table 1.

Water loading ratio (%) Ewater (%) [b] (from mass ratio) 3.0 3.9 5.3 6.6 11.5 16.0 Ewater (%) [b] (improved with MCNPX) 4.9 6.1 8.1 10.3 17.4 23.6

[ a] Correction coefficient used to calculate total energy deposited in zeolite/water systems from energy deposited in the Fricke dosimeter [b] Percentage of total energy deposited in water Table 1. Characteristics of irradiated zeolite 4A / water systems Ampoules were then connected to the analysis bench, and filled with argon, used as gas carrier during irradiations and analysis, with a pressure of 1.46 bar. Each sample underwent several irradiations. The analysis line is also used to replace radiolysis gases after each irradiation and each measure, by pure argon. We then proceeded with the following irradiation.

Irradiations

γ irradiations by a 137 Cs source were carried out with a dose rate of 5.4 Gy.min -1 . The experimental dose was determined by irradiation using the Fricke dosimeter, irradiated under the same conditions as the zeolites samples (equivalent mass, equivalent ampoule and equivalent volume).

After the initial conditioning, consisting in water adsorption at controlled ratios, each ampoule underwent several cycles of irradiations. One cycle includes a filling with argon, used as the carrier gas, exposition to the 137 Cs source, analysis of the gas released, elimination of radiolysis gases, and weighing. This last step is carried out successfully to ensure that the quick elimination of radiolysis gases under primary vacuum did not induce significant modifications of the amount of water adsorbed. Samples have been exposed to radiations during periods from 14 hours to 63 hours. The apparent quantity of hydrogen produced was determined by gas chromatography (Varian model CP 2003), using argon as the carrier gas. Chromatogram peaks were integrated with the software Soprane. Calibration of the chromatograph was performed using a dedicated standard gas, containing argon, and H2, in a known concentration of 99.6 ppm. The error on the H2 concentration measure, in ppm, was calculated to be 2 %.

The absence of any significant hydrogen release under irradiation has been checked for empty glass ampoules and for an ampoule containing dried zeolites 4A.

Results and discussion:

Six systems containing the same quantity of zeolites (about 3.5 g) but with different quantities of water underwent cumulative γ irradiations. The quantity of hydrogen released was measured after each irradiation, and eliminated before the next one.

Dosimetry

The energy deposition is commonly deduced from the Fricke dosimetry [START_REF] Laverne | H2 Production in the Radiolysis of Water on CeO2 and ZrO2[END_REF][START_REF] Laverne | H2 production in the Radiolysis of Aqueous SiO2 Suspensions and Slurries[END_REF][START_REF] Le Caër | Radiolysis of Confined Water: Hydrogen Production at a High Dose Rate[END_REF]. However, if the Fricke dosimeter is relevant for homogenous systems, its precision for heterogeneous systems is debatable. Discussion about energy transfers also requires a precise idea of how primary energy deposition is divided between the zeolites and water. This repartition is frequently evaluated from the weight ratio (LaVerne and [START_REF] Laverne | H2 Production in the Radiolysis of Water on CeO2 and ZrO2[END_REF][START_REF] Laverne | H2 production in the Radiolysis of Aqueous SiO2 Suspensions and Slurries[END_REF], which does not take into account the different composition between the adsorbent and the adsorbed molecules. The different atomic composition between zeolites and water implies, for example, different interaction cross sections, physical and electronic densities. Therefore, we conducted Monte-Carlo simulations, to improve the evaluation of the dose deposition in the systems, in zeolites and in water, rather than considering untreated experimental dosimetry results and mass ratios.

The experimental dose deposed in the samples, is experimentally estimated from the dose deposited in the Fricke dosimeter, exposed to gamma radiation, in similar conditions (equivalent mass, equivalent geometry, and equivalent volume) to the Zeolite / Water Systems (ZWS). The experimental dose rate is thus estimated to be 5.4 Gy.min -1 .

To calculate the energy deposited in the samples, as a function of water loading, we carried out simulations with the MCNPX code, version 2.7.0, based on Monte-Carlo considerations. MCNPX is able to follow gamma particles and secondary electrons generated. It evaluates the statistical energy deposition and location until a threshold of 1 keV. Under this energy, the program consider that the energy carried by the particles is deposited where it is locate. A model of two spheres was used (figure 2). The first sphere contains the water and the second sphere has the interaction properties of zeolites, that is to say, the atomic composition and density (1.57) of zeolites 4A. For each system, the diameter of the water sphere has been calculated to provide a volume ratio between the zeolite and water corresponding to their mass ratio, considering their respective density. If rw, on figure 2, represents the radius of the water sphere and is directly linked to the water mass and volume, R represents the radius of the bigger sphere, but is not directly linked to zeolite mass and volume, as water replaces zeolite material in the geometry. The ratio between rw and R have been calculated for the different quantities of water adsorbed.

We indicated a density of 1 for water and 1.57 for zeolites and atomic composition of zeolites and water is given as an input parameter. Results are given for 10 8 gamma followed. Using the Fricke dosimeter, MCNPX gave a correction coefficient called F (table 1) for energy deposition in the whole sample, slightly evolving with the water loading. Simulation also provides, for each one of the six systems studied, an indication about the energy distribution between the zeolite and the water. This distribution concerns first the energy deposition events and does not take into account the energy transfer occurring after them. The primary energy deposition in water (Ewater)

is given as a percentage in table 1. It does not imply a drastic change compared to primary energy deposition calculated from mass ratios. Note that the energy deposition in percentage calculated from the water mass ratio is different from the Water Loading Ratio. Actually, the latter is conventionally defined as the mass of water divided by the mass of dried zeolites. Energy deposition has to be calculated from the ratio between the mass of water and the total mass of the system.

Hydrogen production

Using these simulation results, we extracted two kinds of information from the hydrogen release, one according to the energy deposition in the whole system, and a second according to the primary energy deposited in the water. The first one consists in considering the extreme case of a complete energy transfer from the zeolite to the adsorbed water. The second one consists in neglecting this energy transfer. Following this second strategy, we obtained a hydrogen release which can be compared to hydrogen measured in free bulk water. This quantity of hydrogen is calculated from the primary radiolytic yield encountered in the literature: 0.045 µmol.J -1 [START_REF] Rotureau | Radiolysis of Confined Water: Molecular Hydrogen Formation[END_REF]. Our measurements first confirmed the enhancement of hydrogen production that can be multiplied by three in the presence of the zeolites, previously described as a catalyst [START_REF] Cecal | Hydrogen yield from water radiolysis in the presence of zeolites[END_REF].

As observed by Nakashima et al. [START_REF] Nakashima | Radiolytic Gas Production from Tritiated Water Adsorbed on Molecular Sieve 5A[END_REF][START_REF] Nakashima | Radiolytic Gas Production from Tritiated Water Adsorbed on Molecular Sieve 5A[END_REF], hydrogen release seems to be linearly dependent on the energy deposited for the higher water loading ratios (higher than 12~13 %) (figure 3). This behavior is conserved for lower loading (figure 4), contrary to the results available for 5A zeolites [START_REF] Nakashima | Radiolytic Gas Production from Tritiated Water Adsorbed on Molecular Sieve 5A[END_REF]) that showed a slow-down of the hydrogen release rate with respect to the dose deposition. In the case of 5A zeolites, irradiations were carried out in one-step, with different irradiations stages. In our case we exposed ZWS repeatedly to irradiations, eliminating hydrogen between the two steps. Otherwise, after our long-step irradiations, for the lower loading ratios, a slight increase is observed for the hydrogen release ratio (figure 4). Similar results are obtained if the energy deposition in the whole ZWS is taken into account.

Hydrogen radiolysis yield

For each water loading ratio, the hydrogen radiolytic yields are extracted from the slopes between two consecutive points. Values obtained for similar slopes are joined to calculate the average hydrogen radiolytic yield. When they are not compatible, two radiolytic yields are calculated and extracted. The calculated hydrogen radiolytic yields, using the two strategies previously described are given in figure 5 and 6. The dashed line in figure 5 shows the reference value for irradiated bulk water irradiated under similar conditions (primary hydrogen radiolytic yield: 0.045 µmol.J -1 ). The dashed line in figure 6 represents the apparent G value calculated from the bulk water radiolytic yield, for the irradiation of a ZWS system, with the total energy deposition taken into account, but without energy transfer from the zeolite.

In both cases, whatever the energy taken into account, hydrogen release seems to be favored in ZWS systems compared to free bulk water. Nevertheless, the yield does not increase linearly with the WLR. A maximum appears and its position depends on whether the radiolytic yield is calculated from the primary energy deposition in the water or in the whole system.

The hydrogen radiolytic yield, calculated from the primary energy deposition in water, shows a maximum for a loading ratio close to 4 %, and then decreases, approaching the yield obtained in water (figure 5). This loading ratio corresponds to a particular hydration condition of the zeolite. At a loading ratio of 4 %, β cages are filled and the water adsorbed is the water which presents the strongest interaction with the zeolite. This is highlighted by the conditions required to evacuate this water. Dehydration of the zeolite from saturation to 4 or 5 % is possible by pumping under secondary vacuum for several hours. Then, to eliminate the residual water, additional heating over 423 K is required [START_REF] Demontis | Dynamical Properties of Confined Water Nanoclusters: Simulation Study of Hydrated Zeolite NaA: Structural and Vibrational Properties[END_REF]. Besides, the entropy of the water adsorbed in β -cages is close to the entropy of solid water at absolute zero temperature [START_REF] Mizota | Entropy of zeolitic water[END_REF]. This illustrates the confinement effect on the organization of water. The strong confinement effect and interaction between zeolite and water at such a loading rate seems to favor energy transfer and dihydrogen production.

The hydrogen radiolytic yields calculated from the energy deposited in the whole ZWS shows a different maximum, at a WLR of 13 % (figure 6). If the transfer of energy from zeolite to water is complete, the hydrogen radiolytic yield should vary in the same way than the radiolytic yield calculated for the primary energy deposition in water. The maximum obtained at 13% highlights a maximum efficiency of all cumulated interactions between zeolite and water. It is coherent with the hydration of the zeolite. Actually, above 13 %, water adsorbed in zeolites, is sometimes called "bulklike" water, as it is located in the center of α cages, and only undergoes a low indirect influence from the zeolite [START_REF] Demontis | Dynamical Properties of Confined Water Nanoclusters: Simulation Study of Hydrated Zeolite NaA: Structural and Vibrational Properties[END_REF]. At about 13 %, the maximum internal surface coverage is reached, leading to a configuration in which the higher proportion of water interacts with the zeolite.

It seems that the more the internal surface of the zeolite is covered by water, the more efficient the energy transfer is. Then, for higher WLR, the hydrogen radiolytic yield decreases tending to join that for bulk water.

If the confinement effects are neglected, we can evaluate the percentage of energy transferred from the zeolite to water at the end of the irradiation stages, by assuming that the hydrogen radiolytic yield of adsorbed water would be equal to the one in bulk water (0.045 µmol.J -1 ).

Results are given in table 2. "p" is the percentage of energy transferred from the zeolite to water, compared to the quantity of energy deposited in the zeolite. "nH2" is the quantity of hydrogen released, "GH2" the radiolytic yield of dihydrogen in water, "Ewater" the primary energy deposition in water and "Ezeolite" the primary energy deposition in zeolite. [a] 6.2 х 10 -3 9.9 х 10 -3 1.2 х 10 -2 1.7 х 10 -2 2.2 х 10 -2 1.7 х 10 -2 [a] : calculated according to the energy deposition in the whole system and excluding the long step irradiations which slow down the hydrogen production Table 2. Evaluation of the energy transfer from the zeolite to the water according to water loading ratios in our systems.

The maximum percentage of energy transferred, evaluated this way, is about 34 %, for the WLR of 13.0 %. This value is high compared to the energy transferred from CeO2 or ZrO2 to adsorbed water, which are of 4 and 18 % respectively, and which are calculated with the same reasoning [START_REF] Laverne | H2 Production in the Radiolysis of Water on CeO2 and ZrO2[END_REF]. This comparison tends to confirm that energy transfers from zeolite to adsorbed water are efficient.

Moreover, the confinement effect has been shown to favor the radiolysis in controlled pore glasses [START_REF] Foley | Radiolysis of Confined Water: Production and Reactivity of Hydroxyl Radicals[END_REF]. Our result, obtained for smaller cavities, is consistent with these observations, as GH2 calculated from the energy deposition in water shows a maximum for a complete filling of the smallest cavity, where the water is in the most organized state. For this loading ratio, according to the strength of interaction between zeolite and adsorbed water, confinement could result in a quasiimmobilization of water. This interpretation is confirmed by a decrease tending to the value of the hydrogen radiolytic yield obtained for bulk water, when the loading ratio increases, and so the strength of the interaction between zeolite and water, and the confinement effects decrease.

Nevertheless, the specific state of confined/fixed water does not simply correspond to a high organization of water, if we compare the radiolysis of liquid water and radiolysis of ice water.

Actually, the radiolytic yield of hydrogen in ice under gamma radiation is much lower than in liquid water, with values such as 0.007 µmol.J -1 at -196°C, 0.010 µmol.J -1 at -100°C and 0.026 µmol.J -1 at -15°C [START_REF] Ghormley | Effects of γ-Radiation on Ice[END_REF][START_REF] Siegel | Isotope Effects in the 77°K γ Irradiation of Ice[END_REF]. Freezing and increasing the organization of water seems to decrease the hydrogen radiolytic yield rather than favoring it. A modification of the migration properties of the intermediate species of radiolysis in confined water could also be implied. It would lead to a perturbation of the recombination mechanisms, resulting in the excessive production of hydrogen compared to free liquid water.

Conclusion:

Hydrogen released has been quantified for different Zeolites 4A/ Water Systems (ZWS), containing different water quantities. An enhancement of hydrogen production in ZWS compared to bulk water has been first confirmed. The influence of the water loading ratio (WLR) on the hydrogen released by radiolysis has been studied with two approaches.

A maximum hydrogen radiolytic yield is observed at a WLR of about 4 %, if the primary energy deposition in water is used for the calculations. In this configuration, water is confined in the β cages, in the most organized configuration, and undergoes the strongest interaction with the zeolites.

It seems that the strongest interaction leads to the highest hydrogen radiolytic yield. β cages water shows the most efficient hydrogen production under radiolysis by gamma radiation.

If the energy deposition in the whole system is used, the maximum hydrogen radiolytic yield is displaced to 13 % WLR. This results reveals the WLR that shows the most efficient hydrogen production, when cumulated interactions are taken into account. The 13 % WLR value, corresponding to a maximum coverage of the inner surfaces of zeolite, is consistent with the hydration of zeolites 4A.

The WLR of 13 % also corresponds to the situation when the gas phase contains the most important concentration of hydrogen, in absolute value, during the radiolysis under gamma radiations.

Figure 1 .

 1 Figure 1. Schematic representation of the bench used to adsorb controlled amounts of water.

Figure 2 .

 2 Figure 2. Geometry used with MCNPX with radius details, to estimate energy depositions in water and in zeolites, from Fricke dosimeter.

Figure 3 .

 3 Figure 3. Cumulated hydrogen released from ZWS, depending on the cumulated primary energy deposition in water, for water loading ratios between 10 and 20%. The dashed line shows the quantity of hydrogen which would have been released for free water irradiated under similar conditions. (WLR: Water Loading Ratio)

Figure 4 .

 4 Figure 4. Hydrogen released from ZWS, depending on the primary energy deposition in water, for the water loading ratio of 5.6%. The dashed line shows the quantity of hydrogen which would have been released for free water irradiated under similar conditions. (WLR: Water Loading Ratio)

Figure 5 .

 5 Figure 5. Apparent hydrogen radiolytic yield as a function of the water loading ratio of zeolites, according to the primary energy deposition in water.

Figure 6 .

 6 Figure 6. Apparent hydrogen radiolytic yield as a function of the water loading ratio of zeolites, according to the energy absorbed by the whole system (zeolites 4A and water)

7.Figure 1 .

 1 Figure 1. Schematic representation of the bench used to adsorb controlled amounts of water.

Figure 2 .

 2 Figure 2. Geometry used with MCNPX with radius details, to estimate energy depositions in water and in zeolites, from Fricke dosimeter.

Figure 3 .Figure 5 .

 35 Figure 3. Cumulated hydrogen released from ZWS, depending on the cumulated primary energy deposition in water, for water loading ratios between 10 and 20%. The dashed line shows the quantity of hydrogen which would have been released for free water irradiated under similar conditions. (WLR: Water Loading Ratio)

Figure 6 .

 6 Figure 6. Apparent hydrogen radiolytic yield as a function of the water loading ratio of zeolites,according to the energy absorbed by the whole system(zeolites 4A and water)

Table 1 .

 1 Characteristics of irradiated zeolite 4A / water systems

Table 2 .

 2 Evaluation of the energy transfer from the zeolite to the water according to water loading ratios in our systems.


J Pontabry
email: pontabry@unistra.fr

F Rousseau
email: rousseau@unistra.fr

C Studholme

M Koob

J.-L Dietemann

A discriminative feature selection approach for shape analysis: application to fetal brain cortical folding study

Keywords: Structural MRI, Fetal imaging, Feature selection, Brain development

The development of post-processing reconstruction techniques has opened new possibilities for the study of in utero fetal brain MRI data. Recent cortical surface analysis have led to the computation of quantitative maps characterizing brain folding of the developing brain. In this paper, we describe a novel feature selection-based approach that can be used for studying the most discriminative brain folding patterns using in utero reconstructed fetal MR data. Results show that this temporal process related to brain maturation can be characterized by a small set of points, located in the main sulci.

Introduction

Development of the central nervous system (CNS) is governed by a complex set of several spatio-temporal mechanisms. Understanding this sequence is fundamental for understanding perinatal neurology [START_REF] Volpe | Neurology of the Newborn[END_REF]. Many processes are involved in brain development at fetal and neonatal stages. These are asynchronous processes taking place at different locations in the brain. For instance, the main period of neural migration is 5 to 25 week postmenstrual age (PMA). Once neurons have been generated, they migrate through two different mechanisms (passive cell displacement and active cell migration). Axon and dendrite sprouting starts around 25 weeks PMA and declines around 1 year. The brain is also the scene of other major events such as synapse formation, glial cell proliferation, myelination, etc. (see de [START_REF] De Graaf-Peters | Ontogeny of 94 the human central nervous system: What is happening when? 95[END_REF] for a recent review of the ontogeny of the human CNS).

As highlighted by de Graaf-Peters and [START_REF] De Graaf-Peters | Ontogeny of 94 the human central nervous system: What is happening when? 95[END_REF], the knowledge on the exact timeline of ontogenetic events occurring during human brain development will provide new insights on the influence of injuries appearing at a specific point in time during this important period of brain building. The modeling of this process timeline from in vivo data is of great importance for improving, for instance, neonatal care services. It is clear that fetal and neonatal periods are key steps of brain devel-expansion [START_REF] Rajagopalan | Mapping directionality specific volume changes using tensor based morphometry: An application to the study of gyrogenesis and lateralization of the human fetal brain[END_REF]. The availability of 1 high resolution 3D images of the fetal brain with tissue 2 labelling enables also studies of tissue boundary shape 3 changes. This has led to the computation of quantita-4 tive maps of brain folding further used for statistical de-5 tection of sulci and brain asymmetry emergence (Habas 6 et al., 2012). Cortical folding analysis has been also used 7 to predict physiological age [START_REF] Wright | Automatic quantifi-109 cation of normal cortical folding patterns from fetal brain MRI[END_REF]. Such [START_REF] Caldairou | Segmentation of the cortex 82 in fetal MRI using a topological model[END_REF]. Let be P = {I 1 , I 2 , . . . , I N } a set of N images. Theses images are normalized in an average space I computed as described previously in section 2.3. The transformation T i that maps the average space I to the image I i is a composition of an affine transform and a non-linear deformation field. In order to capture local shape changes, we only consider the non linear components of the mappings

{T i } i=1,...,N .
Each of the N non linear deformation fields that maps voxels from average space to population space consists of M vectors of R 3 (where M is the number of voxels considered in the region of interest). These 3D vectors are further arranged in a matrix Y ∈ M M×N (R 3 ) where the rows and columns correspond respectively to the features and the samples:

Y =    p 1,1 • • • p 1,N . . . . . . . . . p M,1 • • • p M,N    , (1)
where p i,j is the 3D displacement vector of the voxel i in the deformation field of subject j. Since we focus in this study on brain folding, the analysis is restricted to the displacement vectors belonging to the cortical plate. The matrix Y is the full set of features from which a sparse representation of brain folding is estimated. Instead of adopting a voxel-based approach, that is to say applying a statistical testing method to detect the most significant folding patterns, we propose a data driven approach aiming at extracting from the matrix Y a small set of discriminative features. [START_REF] Tenenbaum | A global geo-100 metric framework for nonlinear dimensionality reduction[END_REF] (see [START_REF] Gerber | 112 Manifold modeling for brain population analysis[END_REF] for an example of Isomap application to manifold modeling for MRI dataset analysis). However, the lower-dimensional space may not be easily interpretable, meaning that in our context, the extracted patterns do not correspond to displacement fields anymore. In this work, we investigate an alternate method by relying on a feature selection approach, aiming at computing a subset of the original variables.

Let Y j be the j th column of the matrix Y and X j a subset vector of Y j (X j ⊂ Y j ) of size M ′ ≤ M . The key assumption here is that the entire set of non linear deformations can be reconstructed using only a very small set of 3D displacement vectors. Within such a sparse framework, Y j and X j are linked by a reconstruction function f :

Y j = f (X j ).
(2)

In this work, we chose to formulate the reconstruction function f using a Nadarya-Watson kernel regression function [START_REF] Nadaraya | On estimating regression[END_REF]:

f (X, h) = N j=1 w j (X, h)Y j (3)
where the weighting function w j is written as

w j (X, h) = K h (X -X j ) N k=1 K h (X -X k ) . ( 4
)
The parameter h is the bandwidth of the kernel K h acting on the subset X. The proposed sparsity-based approach relies on the assumption of redundant information in the original sample Y . Therefore, the relevant information is carried by a subset of parameters. Such a subset of M ′ parameters should minimize the following criterion:

{γ, ĥ} = arg min γ,h J(γ, h) s.t. γ 0 ≤ M ′ , (5)
where γ ∈ M p×1 ({0, 1}) is a binary vector used to define the activated parameters and J is the reconstruction error defined as follows:

J(γ, h) = N j=1 Y j -f (diag(γ) • Y j , h) 2 . ( 6
)
Optimization problems such as defined by equation ( 5) error threshold E.

2

The error criterion is the mean reconstruction error of parameters and is then expressed in mm:

Ē(γ, h) = 1 N p N j=1 p i=1 Y i,i -f (diag(γ) • Y j , h) i 2 , (7)
where f (diag(γ) • Y j , h) i is the reconstruction of the i th

γ := (0 • • • 0) T Estimate parameter h J 0 := J(γ, ĥ) k := 1 while k ≤ M ′ and J(γ, ĥ) ≤ E do J k := J k-1 foreach feature i such that γ(i) = 0 do γ(i) := 1 if J(γ, ĥ) < J k then J k := J(γ, ĥ) i (+) k := i end γ(i) := 0 end if one parameter to add is found then γ i (+) k := 1 Estimate parameter h k := k + 1 else Stop the loop end end X := diag(γ) • Y

Results

5

Experiments have been conducted on both synthetic In each experiment, a Gaussian kernel is used within The synthetic dataset is composed of 20 two-dimensional 19 ellipses (rasterized to images of size 64 × 64 pixels) whose 20 shape varies according to the vertical axis (see Figure 1).

21

The first image of this dataset is then an elongated ellipse 22 while the last image is a circle. This dataset has been gen-23 erated using only one degree of freedom (i.e. the length of 24 the vertical axis). On the first image, the vertical axis is 25 four times smaller than the horizontal one.

26

The results of the proposed algorithm on synthetic of 20 different shapes (reflecting shape variability), or 2) as a set of an evolving shape over time (i.e. temporal evolution of a shape). In later case, it is interesting to study time varying features that could reflect the growth rate for instance, such as length of the vertical axis.

Given the half-size d of the ellipse on the vertical axis and at t = 1, the half-size of the ellipse on the same axis and at t = 20 should be equal to 4dd = 3d, because the vertical axis on the first image is four time smaller than the same axis on the last image. Since there are 20 images sampled over time, the growth rate should be equal to 3d 20 = 0.15d. In our case, the half-size of the vertical axis on the first image is 5.5 pixels. Therefore, the theoretical growth rate of the simulated dataset is 0.825 pixels per time unit.

Figure 4 shows the mean magnitude of vectors over time in the reference space of the first image (t = 1). The linear function fitted to the data has a slope of 0.812, which is close to the theoretical growth rate.

Application on fetal brain MRI dataset

Since our first objective is to study the geometrical changes of the cortical folding through in utero brain maturation, we define as region of interest the cortical gray matter. Moreover, average deformation fields (computed in the orthogonal direction of brain surface) along the cortical gray matter have been used in order to reduce the initial parameter set to a smaller set of 3D displacement vectors [START_REF] Rajagopalan | Mapping directionality specific volume changes using tensor based morphometry: An application to the study of gyrogenesis and lateralization of the human fetal brain[END_REF]. Figure 5 shows five examples of deformation fields considered in this study. Time (a.u.) Mean magnitude (pixels)

Figure 4: Analysis of the growth rate on the synthetical dataset. The plot shows the mean magnitude of selected vectors across time (in the reference space of the first image, t = 1). By studying the fit of a linear function, the growth rate can be estimated with the slope of the function. In this case, the estimated growth rate (0.812) is close the theoretical one (0.825). Our work is also related to the recent study of growth 28 seeds performed by [START_REF] Lefèvre | Identification of growth seeds in the neonate brain through surfacic Helmholtz decomposition[END_REF], where an analysis

Influence of the bandwidth parameter h

8

  brain folding studies are of great importance to improve 9 our understanding of malformation of the cortex. 10 In this work, we focus on the study of the evolution of 11 brain folding during later stages of intrauterine life. In-12 stead of computing scalar features such as surface curva-13 ture to provide global markers of brain development, we 14 have chosen to investigate the use of a data driven ap-15 proach. The underlying question is: what is the smallest 16 most discriminative set of features reflecting fetal brain 17 folding? Following the work of Rajagopalan et al. (2012), 18 a deformation based morphometry approach is adopted to 19 capture directional growth information on the white mat-20 ter / cortex interface. However, the voxel-based statistical 21 testing method applied in (Rajagopalan et al., 2012) does 22 not lead to easily interpretable results of growth patterns. 23 In this study, we have considered the use of feature selec-24 tion techniques in order to extract the sparsest set of de-25 formation fields describing the brain folding process. By 26 studying brain development as a shape modeling problem, we are able to extract the most discriminative set of points 28 related to brain folding during in utero development.

  Subjects and fetal MRI acquisition 31 This study has been conducted on a population of 22 32 fetus aged from 26 to 34 weeks of gestational age at scan 33 time. The data set includes 23 MRI T2 weighted images. 34 Fetal MRI images have been obtained on a 1.5 T Siemens contrast variation induced by motion of the fetal brain 54 with respect to the imaging coil in the magnet. More-55 over, a super-resolution technique is applied in order to 56 remove the effects of the blurring convolution and to in-57 crease the voxel grid density. The resolution of the recon-58 structed images is: 0.74 × 0.74 × 0.74mm. A topological 59 based clustering technique is then applied on the motion-60 compensated high-resolution images to provide segmen-61 tation maps (ventricles, CSF, cortical plate, non-cortical 62 plate -including white matter, subplate, intermediate zone 63 and deep gray nuclei)

  approach for deformation based 66 morphometry to compare anatomies at different stages of 67 brain development, a mapping is estimated to bring every 68 subject's anatomy into correspondence within a common 69 coordinate system. This common space corresponds to 70 the average space of the subjects's anatomies. For this 71 purpose, an iterative registration approach proposed by 72 Guimond et al.[START_REF] Guimond | Average brain models: A convergence study[END_REF] has been used in 73 this work. The procedure to estimate a common space is 74 as following. First, a reference image is chosen among the 75 dataset and the non linear transformations from this refer-76 ence to the remaining images of the dataset are estimated 77 using ANTS diffeomorphic registration technique (Avants 78 and Gee, 2004). Then, the choice of a reference induces 79 a bias that is corrected by using the average of all trans-80 formations. Finally, this average transformation maps the 81 chosen reference to a new space which corresponds to the 82 average space of the subject's anatomies. The reference for 83 the next iteration will be this new space. These steps are 84 repeated iteratively, until convergence of the algorithm, 85 which usually occurs within a few iterations (Guimond 86 et al., 2000). 87 Parameters for ANTS registration are: cross correla-88 tion as similarity measure, gaussian regularization with a 89 symmetric diffeomorphic transformation model. To avoid 90 any possibility of mis-registration due to brain tissue con-91 trast changes, the similarity criterion used during the reg-92 istration process is a weighted combination of image inten-93 sities and tissue label maps (a similar strategy has been 94 proposed by Habas et al. (2010a)). Weights have been set 95 to 0.5 for each feature. 96 To assess the quality of the atlas building step based on 97 image registration, we have computed the average DICE 98 coefficient of the white matter maps between the estimated 99 template and the population. This average DICE coeffi-100 cient is equal to 0.95 ± 0.06, showing the validity of the 101 approach.

103

  Our objective concerns the selection of the most dis-104 criminative deformation-based features reflecting fetal brain 105 development. Please note that although we focus in this 106 work on brain folding, the proposed approach is versatile 107 and it could be used with any other kind of relevant features dedicated to brain anatomy study.

  Such a dimension reduction issue can be tackled with either feature extraction or feature selection techniques. Feature extraction methods transform the high-dimensional data into a space of fewer dimensions. The main linear technique is the principal component analysis (PCA) which performs a linear mapping of the data to a lowerdimensional space in such a way that the variance in the low-dimensional representation is maximized and the covariance is minimized. Although PCA may be applied to study affine deformation, the linear assumption makes it less adapted for non linear deformations. Non-linear feature extraction techniques have been proposed by preserving local data structures. An important example of such non linear techniques is Isomap

3Algorithm 1 :

 1 displacement vector of the j th subject. 4 Feature selection algorithm Input: Y Set of deformation fields P Maximum number of parameters to select E Error threshold Output: X Selected features

6Figure 2 :

 2 Figure 2: Selected points and cost function evolution of the proposed algorithm on a synthetic dataset. The colors reflect the energy of the parameters as a log scale. The first selected parameter (in red) is located on the border of the ellipse and on the varying axe of the dataset. The cost function decreases dramatically after one selected parameter and no parameter can be added without increasing the cost function after 2 iterations.

  11 the reconstruction function (see equation (3)) and the max-12 imum number of parameter is set to M ′ = 100. The mean 13 parameter error threshold is set to E = 0.1 mm for ex-14 periments on brain MR images. For synthetic experiment 15 only, the threshold is set to E = 0 mm in order to study 16 the convergence of the algorithm.

  17

Figure 1 :Figure 3 :

 13 figure 1 can be seen in two different ways: 1) as a set 47

30Figure 5 :

 5 Figure 5: Five instances of the dataset used. The selection algorithm is applied on the cortical deformation fields mapping the reference space to population.

Figure 7 :

 7 Figure 7: The regions of interest used with the feature selection algorithm. These regions corresponds approximately to the frontal (red), temporal (blue), parietal (green) and occipital (yellow) lobes.

  33 few features (up to 6) are required to discriminate the sets 34 of images for each ROI. The coronal view confirms pre-35 vious experiments (see Section 3.2.1): the selected points 36 are located into sulci valleys.

  37 3.2.3. Evolution of deformation vector magnitude 38 We examine in this section the temporal evolution of 39 magnitudes of selected deformation vectors. Once the 40 most discriminant points are selected for each region of 41 interest, we compute the set V of deformation fields be-42 tween each subject and the mean image estimated at 26 43 weeks: V = { v x (t i )} where x is the location of each se-44 lected point and t i is the time difference between the age 45 of the subject i and 26 weeks. The evolution of the mag-

Figure 6 :Figure 8 :

 68 Figure 6: Study of the influence of the choice of the bandwidth parameter of the proposed algorithm. As the bandwidth increases, the number of selected parameters increases as well. An important property is the location stability of the solution with various bandwidth values.

Figure 9 :Figure 11 :

 911 Figure9: Selected cortical points of four lobe regions of the right hemisphere using the proposed algorithm. The left column shows the solution for a fixed bandwidth (h = 2). The other columns show the results for the optimal bandwidth. Top row: frontal lobe; second row: temporal lobe; third row: parietal lobe; last row: occipital lobe. The selected cortical points are depicted as red spheres. The color code in 3D mesh views maps the curvature of the cortical plate from blue color (positive curvature) to red color (negative curvature).

  27

29

  of depth maps of cortical surface based on the Helmholtz 30 decomposition of the deformation fields. In the work of 31[START_REF] Lefèvre | Identification of growth seeds in the neonate brain through surfacic Helmholtz decomposition[END_REF], which is also related to the works of 32 Cuzol et al. (2005) and Grenander et al. (2007), the main 33 idea relies in that a decomposition is used to make the 34 physical interpretation of the deformation fields (captur-35 ing the brain maturation process) easier. Our approach 36 is complementary to these techniques in the sense that 37 the relevant (i.e. most discriminative) displacement vec-38 tors are selected using a non-parametric approach. Further 39 work would consist in analyzing the impact of the defor-40 mation model on the selected features. 41 As shown in our experiments, the most discriminative 42 points lie in deep sulcal regions. These specific locations 43 could be used as consistent shape features in further stud-44 ies. Future work will consist in analyzing these patterns

nitudes of computed deformation fields at selected point 47 locations is estimated using a temporal polynomial of de-48 gree 3 regression approach (see figure 10 for an instance 49 of regression over time for one vector). Figure 11 shows [START_REF] Im | Spatial Distribution of Deep Sulcal Landmarks and Hemispherical Asymmetry on the Cortical Surface[END_REF], especially over time (Meng 55 et al., 2014), may help to set up new early markers of later 56 functional development [START_REF] Dubois | Primary cortical folding in the 104 human newborn: an early marker of later functional development[END_REF] or particular 57 cortical malformation (such as polymicrogyria) that might 58 be difficult to detect and quantify visually.


Aymen Jaziri
email: aymen.jaziri@orange.com

Ridha Nasri
email: ridha.nasri@orange.com

Tijani Chahed
email: tijani.chahed@telecom-sudparis.eu

Traffic Hotspot localization in 3G and 4G wireless networks using OMC metrics

In recent years, there has been an increasing awareness to traffic localization techniques driven by the emergence of heterogeneous networks (HetNet) with small cells deployment and the green networks. The localization of hotspot data traffic with a very high accuracy is indeed of great interest to know where the small cells should be deployed and how can be managed for sleep mode concept. In this paper, we propose a new traffic localization technique based on the combination of different key performance indicators (KPI) extracted from the operation and maintenance center (OMC). The proposed localization algorithm is composed with five main steps; each one corresponds to the determination of traffic weight per area using only one KPI. These KPIs are Timing Advance (TA), Angle of Arrival (AoA), Neighbor cell level, the load of each cell and the Harmonic mean throughput (HMT) versus the Arithmetic mean throughput (AMT). The five KPIs are finally combined by a function taking as variables the values computed from the five steps. By mixing such KPIs, we show that it is possible to lessen significantly the errors of localization in a high precision attaining small cell dimensions.

I. Introduction

One of the most important tasks that a telecommunication professional has is the network planning and optimization especially with the new technologies proposed in 3GPP standards like HetNet deployment [START_REF] Damnjanovic | A survey on 3GPP Heterogeneous Networks[END_REF] and the features of green networks [START_REF] Chen | Fundamental Trade-offs on Green Wireless Networks[END_REF]. RF planning and optimization mainly consist in the configuration and the densification of the network by modifying or adding new small cells and/or changing the network parameters and also saving energy by switching on/off [START_REF] Elayoubi | Optimal control for base station sleep mode in energy efficient radio access networks[END_REF] parts of the network according to the load in each zone. These tasks derive the need of finding hotspot zones [START_REF] Kifle | On The Potential of Traffic Driven Tilt Optimization in LTE-A Networks[END_REF] where the data traffic is very high leading thus to the congestion of the network and the degradation of its performances. Besides, if we don't localize these zones properly -like choosing wrong areasthe worst consequences range from an overall bad system performance (compared to what it could be) to cases where we need more equipment to meet the requirements in the same region. In other words, it implies in loss of CAPEX, OPEX and poor network quality. Heavy traffic localization often consists in using probes. Measurements are performed to optimize the network configuration and deploy new small cells or switch on/off some cells. However, it is possible to identify the heavy traffic areas in the network using KPIs taken from the OMC. The accuracy of localization can be very high and improved comparing to the recent used techniques. As a new alternative, key performance indicators can give an idea about the hotspot zones but with very low accuracy (with a granularity of the cell radius). Nevertheless, the combination of many KPIs would significantly increase the precision of traffic localization. The main contribution of this paper is to combine several KPIs according to a given new algorithm. This new algorithm is not clearly disclosed because of the intellectual property. Using this new method, the precision of localization is increased and the costs are significantly reduced compared to probes-based methods. Several traffic location techniques have been proposed to detect the hotspot zones of the data traffic. Currently, the operational methods used for optimization are based on probing, trace analysis and protocol decoding. Probing is an expensive solution because it needs tools for processing probes, storage servers... In the other hand, some solutions are under studies in the literature in order to give the highest accuracy with the lowest costs. To be based on the real data traffic evolution, authors in [START_REF] Ahmed | Traffic location in mobile cellular telecommunications systems[END_REF] provide a test transmitter which plays the role of a neighbor cell at first and then the role of a serving cell. Tests are realized within the existing cells in the network in order to assess the traffic density within the vicinity of the transmitter means. This solution needs to realize measurements in each area separately rendering it cumbersome and also comparable to probing in terms of cost. In the same context, the patent [START_REF] Ho | Process for analysing traffic location in a cellular radio communication network[END_REF] discloses a method where the user equipments (UEs) send periodically a radio measurement report taken from both the serving cell and the neighbor cells. A recording unit is installed on the interface between the base station (BS) and studied base station controller (BSC) interface (A-BIS), to examine the messages exchanged on this interface, and based on these measurements, the traffic distribution is calculated. This method is not very accurate and the dimension of precision doesn't go to the small cell dimensions. In [START_REF] Randriamasy | Geographic representation of traffic load in a mobile radio communication network[END_REF], authors propose to divide each cell of the network in a plurality of areas, and calculate a traffic value for each of these areas from the measurements in the network. The useful indicators to compute this traffic values are the channel occupation time, the call attempts, the number of handovers. This method is very simple since it is easy to get the inputs needed and it doesn't need additional tests or equipment. However, the localization of hotspots with this method is not very accurate and also the proposed optimization problem could be hard to solve. Another alternative is presented in the patent [START_REF] Vaara | Traffic hot spot locating method[END_REF] which treats the hotspot traffic localization via statistical analysis of the timing advance measurements. Then, the accuracy is increased with providing the related neighbor cell measurements. These measurements are collected and processed in a base station controller BSC and the related accuracy is better than the previous solutions but still not enough to meet with dimensions of small cells. With probing like in [START_REF] Ahmed | Traffic location in mobile cellular telecommunications systems[END_REF] [6] [START_REF] Vaara | Traffic hot spot locating method[END_REF], the costs are important due to the equipment used to collect information from calls and the need for servers to store the traces. Also, the probes could not collect all the traffic transmissions, some of the transmissions are not captured and lost, so it is still not too efficient for traffic localization. However, a solution based on KPIs like in [START_REF] Randriamasy | Geographic representation of traffic load in a mobile radio communication network[END_REF] will not cost too much; we have only to provide an algorithm to localize the heavy traffic with high accuracy. Note that the localization of each UE is possible and with a good accuracy thanks to many techniques cited in [START_REF] Roxin | Survey of Wireless Geolocation Techniques[END_REF]. For example, triangulation with AoA measurement (LTE Rel 9) is used and means that the AoA of the UE signal is calculated from two reference points which are generally the position of two BSs then the location is determined by the intersection of the two lines bearing from the known reference points. However, these individual locating methods are cumbersome and consume a lot of time to extract from them the distribution of traffic data transmission because it deals with a large number of users in the network. This paper is organized as follows. First, the KPIs, from which we extract the information of traffic distribution in the network, are defined. Next, a modular algorithm is proposed and in each step, only one KPI is analyzed before combining them in the last step. Then, the performance of the estimation is analyzed by simulating an LTE network with the configuration recommended in 3GPP-LTE A standard.

II. Description and design of the traffic localization algorithm

A. Network used metrics

In order to better understand the utility of all the aforementioned KPIs, each KPI is calculated and can be obtained a follows:

1) Timing Advance: For the timing advance (in GSM, LTE and LTE-A) or the propagation delay (in WCDMA), a time offset measurement is realized by the BS between its own transmission and the transmission received from the UE. According to the calculated offset, the serving BS determines the suitable timing advance for the UE [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA)[END_REF]. Then, from this timing advance, the BS calculates the distance traveled by the radio signal. In practice, depending on the resolution (or granularity) of timing advance, a specific distance range where the UE is located will be calculated. Actually, timing advance is used as a KPI for network supervision and analysis and with a granularity of 78.125 meters in LTE networks [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA)[END_REF]. For the timing advance, each cell is divided into intervals according to the interval of distances like in figure 1. Then the distribution of the UEs in the cell is calculated according to these intervals.For example, 30% of the UEs are in the range of TA1, 20% in the range of TA2, 40% in the range of TA3 and 10% in the range of TA4. 2) Angle of Arrival: Angle of Arrival (AoA) is defined as the estimated angle of a UE with respect to a reference direction which is defined as geographical North; the value of AoA is positive in an anticlockwise direction [11]. In general, any uplink signal from the MS can be used to estimate the AoA, but typically a known signal such as the Sounding Reference Signals (SRSs) or DeModulation Reference Signals (DM-RS) would be used [START_REF] Sesia | LTE -The UMTS Long Term Evolution: From Theory to Practice[END_REF]. In order to construct the shape of the downlink beam, direction of arrival estimation already exists and is supported in WCDMA but it is not standardized [START_REF] Singla | Beam forming Algorithm for Smart Antenna in WCDMA Network[END_REF]. According to the angle of arrival, the cell is divided into ranges like in figure 2 relative to the angle between the UE, the BS and the geographical north. With this classification, the corresponding spatial distribution is calculated. An example of distribution that we can have is 30% of UEs with range corresponding to AoA1, 40% for AoA2 and 30% for AoA3.

3) Neighboring Cell Level: In each access or handover process, every UE measures the signals of the detected cells and sends a report of these measurements to its serving cell. In GSM, the detected cells list is reported periodically by the UE to the BS. However, in 3G and 4G networks, it is reported to the BS only when a special event is activated like handovers or cell selection. For each neighboring cell, a percentage of UEs, connected to a serving cell and declare this neighboring cell as the first best candidate for handover, is calculated. The corresponding KPI is the percentage of UEs declaring the same neighbor cell as their second best serving cell (best candidate cell for handover). The information extracted from this KPI is exploited to recover the errors of AoA estimation which is due to reflexions in real networks. The exploitation of this KPI is equivalent to the triangulation technique used in the probes-based methods.

4) Load Time: The load time represents the percentage of time when the serving cell's resources are fully occupied. This KPI is a gauge counter (means it is a percentage in time) already defined in the actual tools and is calculated during busy hours. In hotspot zones, the traffic evolution almost has the same behavior in the neighbor cells due to the load balancing and forced handovers [START_REF] Nasri | Handover Adaptation for Dynamic Load Balancing in 3GPP Long Term Evolution Systems[END_REF]. Then, if a cell is congested and one of its neighbors is congested, it is most probable that most of the traffic is located in the edge between the two cells and if there is no correlation with any of the neighbor cells, the traffic is most probably generated from users close to the serving cell.

5) Harmonic and Arithmetic Mean Throughput :

The harmonic mean throughput gives more importance to the UEs with bad radio conditions and which are generally in the cell edge. From the formula of the cell throughput defined into the OMC by some constructors, it is important to notice that the cell throughput has approximately the same value as the harmonic mean throughput. On the other side, the arithmetic mean throughput is the mean of the throughputs of all the connected users. The utility of these KPIs comes from the following possible scenarios; if most of the UEs are close to the BS (so, in very good radio conditions), the arithmetic mean throughput is high and the harmonic mean throughput is low. In this case, the difference between the harmonic and the arithmetic mean throughput will be high. However, when most of the UEs are in the cell edge, the arithmetic mean throughput decreases and the harmonic mean throughput increases comparing to the first scenario. So, the difference between the two mean throughput is low. In sum, the difference between the arithmetic and the harmonic mean throughput gives an idea about the location of the hotspot of traffic. This KPI is expected to give the same information like the correlation of load time but with different parameters and this can ameliorate more the estimation and correct the errors due to time dispersion or reflexions problems.

B. Algorithm formulation

Estimating the hotspot zones in the network is realized as follows: First, the OMC reports the statistics of each KPI corresponding to the cells in its controlled area. Then, each KPI is analyzed into an algorithm in order to extract weights indicating the spatial distribution of the traffic.

Step 1: calculate the spatial distribution of weights according to the timing advance: Actually, the cell is divided into many areas depending on the distance (or the timing advance) between the pixel's position and the serving BS's position. Then, each pixel takes as a weight the percentage of UEs in its range of TA.

Step 2: calculate the spatial distribution of weights according to the Angle of arrival: Each cell divides its covered area into many sub-areas depending on the angle between each pixel's position, the serving BS's position and the geographical north. So, each pixel takes as a weight the percentage of UEs in the same range of angle of arrival defined in the corresponding KPI.

Step 3: calculate the distribution of weights according to the Detected neighbor cell: According to the neighboring cell level, the weight in each pixel will be the percentage of UEs reporting as a second best serving cell the same one as the cell having the best RSRP measurement in this pixel. It is indeed possible to define an RSRP map, and define the second best serving cell for each pixel with the RF fingerprinting technique [START_REF] Pereira | Evaluating location fingerprinting methods for underground GSM networks deployed over Leaky Feeder[END_REF]. In this context, RF fingerprinting consists in matching each pixel to a range of statistics pertaining to RF topology. With the grid calculated by this technique, second best serving cell is defined in each pixel.

Step

4: calculate the distribution of weights according to the traffic load time and correlation between cells:

The correlation in load time between cells is tested and the distribution of weights is calculated. In the following step, for each pixel, we identify the set of cells which are candidate to be the serving cell and with a load time near to the load time of the best serving cell. Then each pixel in the network map takes as a weight the sum of all the load time of the set divided by the cardinality of this set.

Step 5: calculate the distribution of weights according to the Arithmetic mean throughput vs Harmonic mean throughput:

The weights are the difference between the arithmetic mean throughput of the cell and the harmonic mean throughput. This quantity is divided by a constant in order to have all the calculated weights in [0 , 1]. The weight in each pixel depends on its position relative to the position of its best serving cell and also on the value of the difference between the two mean throughput. So, the pixels in the cell edge take the value of the difference between the HMT and the AMT divided by a constant if this difference is less than a certain threshold which should be calculated properly and they take zero if this condition is not satisfied. In contrast, the weights in the pixels close to the BS take the value of the difference between the HMT and the AMT divided by a constant if this difference is higher than a certain threshold and zero else.

Step 6: include all KPIs in a metric and estimate the hotspot traffic zones: After analyzing all the KPIs and extracting values from them, the most important step is how to transpose all the outputs from the previous steps of the main algorithm into the same map. In fact, a metric including all the KPIs is defined with the different weights depending on the importance of the information given by each KPI. The metric used can be any function including all these projections and it is possible to give more importance to some KPIs and to reduce the impact of the other KPIs.

Finding the best function can be an optimization problem that gives the optimal coefficients for the function. In order to improve the estimation of the traffic distribution in the network, the values estimated in each pixel are combined with the values calculated in the neighboring pixels. The function of the smoother is distance decay because the correlation between two pixels is smaller when the distance between them is high.

III. Simulation and analysis A. Simulation

In order to test the proposed algorithm, we consider a realistic network with 23 tri-sector sites (see figure 3) with the available bandwidth of 20MHz. At each second, UEs are generated depending on the arrival rate. We suppose that each UE has a file of size 1000 Kbit to download. Moreover, a part of UEs are moving during the transmission with mobility of 8.33 km/h meaning that the simulation supports also events of handover. UEs are scheduled according to the round robin model. Finally, KPIs are updated in the network and are filtered approximately every quarter of an hour which corresponds to the time granularity of most OMC vendors.

B. Numerical results

In figure 4, we plot the real distribution of the traffic. In figure 5, the estimated traffic distribution is represented. Thus, the first conclusion is that the hotspot zones are found in the entire network map and the accuracy of the estimation is 25 meters. The figure 6 shows that using all the five KPIs in a real network topology gives better estimation of traffic distribution comparing to the case when using some of them. Moreover, from the figure 6, the TA, AoA and the neighbor cell level are useful KPIs but each one is with low accuracy, however the Mean throughput and the correlation of the load time are less important but they have also an impact on the estimation and they improve significantly the estimation when using the other KPIs and it doesn't improve the estimation when using some of these KPIs. The error of estimation using all the proposed KPIs is less than 10% while in the case where one of the KPIs is eliminated; the error becomes more than 12%. For a moderate traffic, the performance of our estimation can be observed in the figure 6. However, there stills some problems related to these KPIs especially when the data traffic is heavy. So, some of the KPIs can be biased by the UEs not accepted in the network and trying again and again to be served. This problem can be treated with analyzing only the performance in term of elapsed traffic. Then, the figure 7 represents the CDF of elapsed traffic during the simulation. As a result, simulations proved that exploiting all the proposed KPIs gives the smallest error of less than 8% in the estimation of elapsed traffic and the percentage of error is less than the case where one of these KPIs is eliminated.

We compare in the following table the percentage of detected hotspots which has the highest weights to the real hotspots in the network. From the previous table, we observe that the first hotspots which have the highest weights are well estimated. In other words, the weights attributed to 0.5% of most important hotspot zones represent 0.42% hotspot zones in the estimation. However, when the threshold of detection is decreased, the estimation presents less efficiency and this is due to the fact that some zones near to the hotspots could take a significant weight without carrying really heavy traffic. In the other side, the results of the proposed algorithm are better than in the case where only TA is exploited which is a solution proposed in several studies [START_REF] Vaara | Traffic hot spot locating method[END_REF].

In sum, the simulations show encouraging results and the proposed solution is very promising to meet with the exact results with optimizing the function including all the KPIs and 10% of error in the hotspot zones is near to the percentage of error generated by the actual used techniques such as probing which costs too much in term of OPEX and CAPEX.

IV. Conclusion

Traffic hotspot localization is an important issue in the network densification and also in the efficient policy of switching on/off small cells. Therefore, proposing a new technique based only on OMC metrics could be a very interesting solution that can save more CAPEX and OPEX and give a good precision. In this work, a new algorithm consisting in exploiting KPIs is built and the accuracy that it provides is high comparing to the actual and operational tools. It can be considered as a new tool for dimensioning and planning small cells. KPIs' based location solution shows a good precision with a good level of certainty. For future works, we aim to find the suitable function including the KPIs and this can be done by formulating an optimization problem which gives the best function mixing these KPIs with the optimal weights. Furthermore, potential hotspots and additional information layer about commercial or industrial zones would be treated in order to optimize much more the proposed localization technique.

Fig. 1 .

 1 Fig. 1. The division of the covered area according to TA

Fig. 2 .

 2 Fig. 2. The division of the covered area according to AoA

Fig. 3 .

 3 Fig. 3. Map of the network layout and the best serving cell in each pixel

Fig. 4 .

 4 Fig. 4. Real data distribution

Fig. 5 .

 5 Fig. 5. Estimated data distribution

Fig. 6 .

 6 Fig. 6. CDF of real and estimated weights in term of access demand

Table :

 : Percentage of detected hotspots

		1
		0.9
		0.8
		0.7
		0.6					Real weights
							TA
	F(x)	0.5					AoA Neighbor cell level
		0.4					TA and AoA
							TA and Neighbor cell level
		0.3					AoA and Neighbor cell level
							TA, AoA and Neighbor cell level
		0.2					TA, AoA, Neighbor cell level and Load
							TA, AoA, Neighbor cell level and Mean throughput
		0.1					All KPIs
		0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								x				x 10 -3
	Fig. 7. CDF of real and estimated weights in term of elapsed traffic
			Real % of first Estimated % with Estimated % with
			hotspots			all KPIs			only TA
			0.5%			0.42%			0.12%
			1%				0.57%			0.26%
			2%				1.1%				0.47%
			5%				2.51%			1.12%
			10%			4.74%			2.29%
			20%			9.55%			4.87%
			50%			26.92%			14.22%
			70%			41.36%			23.23%


1) Inception of localization of ductile deformation within the Gialiskari detachment located between the Kakinagrion granite related to migmatite and the Upper Cycladic unit (UC).

2) Evolution of the deformation through the brittle-ductile transition. Cataclastic deformation starts to localize in the Gialiskari detachment.

3) Intensive localization of cataclastic deformation in the Gialiskari detachment which is overprinted on the ductile deformation.

4) Distribution and quantification of ductile deformation

The Raches and Karkinagrion granites show intense ductile deformation. In the field, a gradient of strain has been observed from south to north. In order to characterize this gradient, we used a scale of ductile deformation intensity, divided into seven grades. Each grade of deformation was based on macro-structural criteria adapted to fieldwork.

In addition, we have reconstructed by image analysis the strain ellipsoid of the Raches granite based on the shape fabric of biotites. A good correlation was observed between the degree of anisotropy (P') and the quantitative map of deformation elaborated in the field. ii) The Gialiskari detachment.

iii) The contact between the Raches and Karkinagrion granites. Our new mapping, reveal a much larger volume of the Karkinagrion granite than previously described by Ring (2007) andBolhar et al. (2010a). We also highlight the presence of migmatites related to the Karkinagrion granite. These observations led to redefine the nature of the contact between these two granites described by Bolhar et al. (2010a). It rather seemed that the Raches granite intruded the complex migmatites-Karkinagrion granite and not the opposite.

The objective of this work was to study interactions between syn-tectonic granites and MCCs development. Based mainly on field observations, we conclude that the Karkinagrion and Raches granites emplacement postdated the initiation of the Ikaria MCC. While the exhumation of the metamorphic dome involved the Agios Kirykos detachment during the first stages, the plutons intruded this detachment and deformation migrated upward when granites were emplaced.

This complex interaction between detachments and intrusions were found also on Tinos and Mykonos within the NCDS. It suggested that the plutons did not initially localize detachments, instead that the detachments localized and controlled the ascent of the plutons. The intrusions then interacted with detachments, making them migrating upward.

2) Geological setting

The Aegean domain is an ideal place to investigate the development of MCCs and to study the role of syn-tectonic granites on their development.

Cycladic islands are located in the center of the Aegean domain and correspond to the deepest exhumed part of the Hellenides-Taurides chain. Three tectonic units are recognized in the Cyclades:

1) the Upper Cycladic unit (UC)

2) The Cycladic BlueSchist unit (CBS)

3) The Cycladic Continental Basal unit (CCB). Shear zones control the ascent and emplacement of plutons.

Model 1: Brown and Solar (1998a)

Emplacement of plutons allow the initiation of shear zones.

1) Introduction

The goal of the present study was to investigate the interactions between syn-tectonic granites and Metamorphic Core Complexes (MCCs) development. Two different models were described:

-Model 1: Shear zones control the ascent and emplacement of plutons.

-Model 2: Emplacement of plutons allow the initiation of shear zones.

Interactions between extensional shear zones and syn-tectonic granite intrusions: the example of Ikaria Island (Cyclades, Greece)

ii) Gialiskari detachment iii) Contact between Raches and Karkinagrion granites

Valentin LAURENT 1,2,3 , Alexandre BEAUDOIN 1,2,3 , Laurent JOLIVET 1,2,3 , Laurent ARBARET 1,2,3 , Romain AUGIER 1,2,3 , Aurélien RABILLARD 1,2,3 (1) Université d'Orleans, ISTO, UMR 7327, 45071, Orléans, France, (2) CNRS/INSU, ISTO, UMR 7327, 45071 Orléans, France, (3) BRGM, ISTO, UMR 7327, BP 36009, 45060 Orléans, France B771

  deformation reaching the upper part of the crust.

  -scale mapping and major structures A new field study of the granitic plutons of Ikaria Island (Cyclades, Greece) was conducted and resulted in the proposition of a new geological map. Our survey shows three major structures : i) The intrusive contact between the Raches granite and the metamorphic dome, in opposition to a tectonic contact like described by Photiades (2002a), Kokkalas and Aydin (2013) or Kumerics et al. (2005).

  controlling the ascent and emplacement of plutons.