Nabil Amari

David Folio

Antoine Ferreira

Robust Nanomanipulation Control based on Laser Beam Feedback

This paper presents a study of the control problem of a laser beam illuminating and focusing a microobject subjected to dynamic disturbances using light intensity for feedback only. The main idea is to guide and track the beam with a hybrid micro/nanomanipulator which is driven by a control signal generated by processing the beam intensity sensed by a four-quadrant photodiode. Since the pointing location of the beam depends on real-time control issues related to temperature variation, vibrations, output intensity control, and collimation of the light output, the 2-D beam location to the photodiode sensor measurement output is estimated in real-time. To overcome this problem, the Particle Filter (PF) algorithm is used to estimate the position of laser beam. The dual manipulators are controlled by combining different performance dynamics (micro and nano manipulators) in order to track a laser beam with very-high precision in a dynamics operating mode.

I. INTRODUCTION

Beam analysis (laser beam, or X-ray radiation) represents a great tool to improve the comprehension of physical phenomena of different types of objects at the micro or nanoscale [START_REF] Stuhrmann | Automated tracking and laser micromanipulation of motile cells[END_REF] [START_REF] Nahmias | Analysis of Radiation Forces in Laser Trapping and Laser-Guided Direct Writing Applications[END_REF]. A typical application is in biology where it is needed to diagnose different pathologies of cellular samples using an infrared beam. However, the main difficulty is to maintain the microsample in the focus of the laser beam with a tradeoff between robustness of the tracking and precision of the focal point (e.g., in few micrometer size spots). Usually, an atomic force microscope (AFM) nanogripper constituted by two AFM tips is used to handle and to maintain the cell under the field of focus of the laser beam [START_REF] Sitti | Atomic Force Microscope Probe Based Controlled Pushing for Nanotribological Characterization[END_REF][4] [START_REF] Clévy | A micro-manipulation cell including a tool changer[END_REF]. For characterization purposes, the sample should be maintained in a stable way during for long periods of time with the capacity to reject the external disturbances. These perturbations are produced by the nanomanipulator platform that is subjected to mechanical microvibrations, actuator thermal drifts, photodetector noises, and brownian motion of the laser beam. Furthermore, the atmospheric turbulences are critical for the success of such AFM-based nanomanipulation tasks [START_REF] Shen | Tracking Nanometer-Scale Fluorescent Particles in Two Dimensions With a Confocal Microscope[END_REF][7] [START_REF] Nestor | Observer-Based Intensity-Feedbck Control for Laser Beam Pointing and Tracking[END_REF]. Classically, the task sequencing is defined by three main tasks (see Fig. 1): i) detection and selection of the microobject, ii) stable nanohandling, and iii) transportation of the microobject in the focus of the laser beam. Such procedures usually involve a two-fingered AFM-based nanomanipulation system with different types of coarse and fine stages with different bandwidth dynamics, position precision, travel range (from few nanometers to hundreds micrometers) and actuation principles (magnetic and piezoelectric-stack N. Amari, D. Folio and A. Ferreira are with the Laboratoire PRISME EA 4229; Ecole Nationale Supérieure d'Ingénieurs de Bourges, 88 boulevard Lahitolle, 18020 Bourges, France. Corresponding author: Antoine Ferreira (Email:antoine.ferreira@ensi-bourges.fr, Tel: +33 2 4848 4079) Fig. 1. High-precision position measurement based on laser beam transmission are commonly used in nanomanipulation applications. actuated nanomanipulators) [START_REF] Iwata | Nanomanipulation of biological samples using a compact atomic force microscope under scanning electron microscope observation[END_REF][10] [START_REF] Seo | Dual Servo Control of a High-Tilt 3-DOF Microparallel Positioning Platform[END_REF]. Different control architectures have been developed for dual-stage servo control schemes to deal with uncertainties (modeling errors, sensor limitations) and environment noises [12][13]. They can be classified into two categories: those based on decoupled or sequential single-input-single-output (SISO) designs, and those based on modern optimal design methodologies, such as LQG, LQG/LTR, H ∞ , and µ-synthesis, in which the dual-stage controllers are obtained simultaneously. Another problem to deal with is to ensure that the microsample is maintained in the focus of the beam that possess a section area of few squared micrometers. Due to brownian motion and micromechanical vibrations perturbations, the beam tracking is a difficult problem since the microsample should be illuminated by the maximum beam intensity at the focal point (as shown in Fig. 1)) during the characterization task.

We studied in this work the robust control issues of the dual AFM-based micro/nano manipulator motions to ensure the tracking of the laser beam with micro/nanometric resolution. First, robust control strategies for magnetic and piezoelectric-stack actuated nanomanipulators are implemented to deal with uncertainties (modeling errors, sensor limitations) and environment noise. Then, efficient robust algorithms are proposed to track the handling position variations due to beam exposition (electrostatic forces, Brownian motion, scattering) based on the Particle Filter (PF) algorithm. Secondly, a H ∞ controller is designed to deal with uncertainties (modeling errors, sensor limitation, non linear effects) and environment noise. Moreover, the cooperative control scheme has to coordinate two nanogrippers mounted in real-time on a 6 degree of freedom (dof) dual-stage micro/nanomanipulator system.

The paper is divided into five sections. Section 2 describes the experimental setup composed of a two-fingered AFMbased nanomanipulation system operating under the field of focus of a laser beam emitting at infrared. Then, section 3 describes the dynamics modeling of the system. Section 4 describes the master-slave controller scheme with a decoupling structure for the laser beam tracking and nanogripper handling. Finally, section 5 presents experimental results that demonstrates the robust performances of the laser beam tracking system.

II. EXPERIMENTAL SETUP

The two-fingered AFM-based nanomanipulation system is illustrated in Fig. 2. First, the handling system comprises two atomic force microscope (AFM) probes (Kleindiek FMT-400) each mounted on two computer-controllable micro/nano manipulators facing each other. Each micro/nano manipulator is composed of 6 DoF high-precision dual-stages: a x-yz closed-loop nanostage (P-611.3S NanoCube from Physics Instruments) with a fine motion in the range of 120 × 120 × 120 µm which is mounted on a x-y-z closed-loop microstage (F-130 DC from Physics Instruments) with a coarse motion in the range of 15 × 15 × 15 cm. On the side view, the laser beam is generated from a 635 nm laser diode (Red), and a four-quadrant position sensing device (PSD) measures the position of the image that the laser beam forms on a fixed plane. The imaging system is composed of a top-view (optical microscope -Mituyo ×50) and side-view (digital microscope -TIMM ×400) used for localization and guidance. Fig. 3 shows the control scheme of the two AFM probes to detect, grab, transport a microball in the focus of the laser beam, and the laser beam motion tracking (Brownian or stochastic trajectory) are processed in real-time using MATLAB © xPC software. For laser beam intensity maximum detection and beam tracking a highspeed data acquisition (DAQ) (NI 6289) card is used to register photodiode voltage output, generated when the PSD is illuminated. A multi-thread structure is developed to independently control and coordinate the two AFM tips for parallel manipulation.

III. DYNAMICS MODELING

This section reviews the different model dynamics of the different system components.

A. Dynamics of Piezoelectric and Magnetic Actuators

The micromanipulation process needs the cooperative control of both AFM probes in order to handle and to track, in real-time, the microsample under the field of focus of the laser beam detected by the PSD. The cooperative controller sends commands to the micro/nano stages along the x-yz directions. As no parametric information on drivers are available, an identification phase is needed to set up dynamic modeling of dual micro/nanostages. We identified the dynamics of our dual micro/nano manipulators in x-y-z-directions where the dynamic models of the micro and nanostages are chosen as a third-order approximation reduced order for each x-y-z axis, respectively:

G micro(x,y,z) (z) = b 0 + b 1 z -1 + b 2 z -2 + b 3 z -3 1 + a 1 z -1 + a 2 z -2 + a 3 z -3 G nano(x,y,z) (z) = b 1 z -1 + b 2 z -2 + b 3 z -3 1 + a 1 z -1 + a 2 z -2 + a 3 z -3 (1)

B. Dynamics of Four Quadrant Detector and Laser Beam

A four quadrant photo sensitive detector (PSD) has four photosensing parts arranged in four quadrants, respectively. Usually, the laser beam is pointed towards the dead center between the 4 quadrants and the beam diameter is chosen to fit within the PSD area. When light falls on all PSD quadrants, they generate currents for each quadrant according to the light intensity and then amplified into voltage signals V 1 to V 4 . The difference between the left and right quadrants (V x ) and top and bottom quadrants (V y ) can be used to indicate the offsets of the spot and be adjusted to zero by centering the beam, whereas the sum quadrants voltages corresponding at intensity laser beam (V s ) is at a maximum, that is:

V x = (V 1 + V 4 ) -(V 2 + V 3 ) V y = (V 1 + V 2 ) -(V 3 + V 4 ) V s = V 1 + V 2 + V 3 + V 4 . (2)
Fig. 4 illustrates the voltage signal V x (V y is similar as V x ) and V s . Thus V x and V y channel outputs are directly related to the energy of the laser beam that falls in each quadrant. It is assume that the light intensity on the laser's beam cross section obeys Gaussian distribution. The current generated CONFIDENTIAL. Limited circulation. For review only. by each sensing element can be described as given in:

I = k 1 2E l π 2 r e 2(x 2 1 +y 2 1 ) r 2 dx 1 dy 1 ( 3
)
where I is the current, r the radius of the laser light spot, E t is the energy of the laser beam, (x 1 ,y 1 ) is the coordinate of a point on the light spot in a coordinates system located at the center of the light spot, and k 1 is a coefficient. Furthermore, the AFM probe can be used to measure changes in the beam intensity, and used to correct the V x and V y output values for voltage changes due to intensity fluctuations rather than actual beam deviations by filtering. Moreover, the laser beam motion is assumed similar to the Brownian motion of a particle subjected to excitation and frictional forces. The Brownian motion is given by the generalized differential equation:

d 2 x(t) dt 2 + β x dx(t) dt = W x (4)
were β x coefficient of friction and W x is zero-mean Gaussian random variable with variance δ 2

x . The y-axis can be modeled in the same manner as the x-axis, though with different dynamics. The discrete state space of (4) laser beam is given by:

X k = AX k-1 + BW k (5)
Y k = CX k-1 (6)
with

X k = [x k y k ẋk ẏk ] T , A =      1 0 ∆T 0 0 1 0 ∆T 0 0 a x 0 0 0 0 a y      , B =

IV. CONTROL SCHEME OF A TWO-FINGERED NANOMANIPULATION SYSTEM TRACKING A LASER SPOT

The analysis of the robotic tasks demonstrates that the handling task, realized by the two-fingered nanogripper, necessitates that the maximum beam intensity should be focused on the microsample in a stable and robust way during several minutes without releasing it. It requires a robust control motion controller of the dual micro/nano manipulators coupled to a laser beam tracking controller to compensate localization uncertainties of laser beam motion. Fig. 5 presents the block diagram of the dual-stage controller using the proposed decoupled control structure. A stochastic filter is used to detect and estimate the maximum laser beam intensity in order to anticipate the intensity variations.

A. Particle filter approach

The first step is to guaranty a real time estimation of the maximum beam intensity during motion. For that, we use a Particle Filer(PF). It approximates the real distribution by generating a set N samples x i k distributed according to the posterior distribution, p(x t |y t ), and they associate a normalized importance weight. The adaptive multinominal resampling [START_REF] Van Der Merwe | The Unscented Particle Filter[END_REF]., based on bootstrap method [START_REF] Efron | An Introduction to the Bootstrap[END_REF], is used to estimate and predicate the laser beam positions. The corresponding PF approach is presented in Algorithm 1.

Algorithm 1 SIR particle filter approach 1: Generates N samples {x i 0 } N i=1 form the initial distribution 2: Calculate the weights

w i k αp(y k |x i k ), i = 1, ...., N and normalize w(i) k = w (i) k / N j=1 w (j) k 3: Generate a new set {x i * k } N i=1
by performing resampling 4: Predict new particles x i k+1 = f (x i * k , v k ), using different noise realization for the particles.

5: Compute the output of the SIR filter by:x

N k = N i=1 w i k x i k B. H ∞ Controller
The MIMO H ∞ standard controller is adopted here in order to provide robustness performances against model uncertainties(modeling errors, sensor limitations, nonlinear effects) and environment noises [START_REF] Ahmad | H∞ Controller Design for High-Preformance Scanning Tunneling Microscope[END_REF] of the piezoelectric actuated nanomanipulators. The magnetic actuated micromanipulators are controlled using a linear PID controller. The Fig. 6 shows the H ∞ block diagram where W i represents the weighting function. The H ∞ standard problem consists to find a controller K which stabilizes the system and determine a positive number γ 0 to satisfy the following condition:

T (s) ∞ < γ(γ ≥ γ 0 ) (7)
where T (s) represents the closed-loop transfer function of the system represented in Fig. 6, given by:

T zw = W 1 S W 1 SGW 3 W 2 KS W 2 KSGW 3 (8) S = (I + G(s)K(s)) -1
is the sensitivity function. To solve the H ∞ standard problem, we need the appropriate choice of the structure weighting functions with optimal parameters which reflect the robust stability and performance requirements. Considering that the axis coupling is negligible, lets us to choose:

W i = diag(w i1 , w i2 , w i3 ) , with i = 1 : 3 (9)
and

w 1i = a1ip+a2i p+b1i w 3i = p+a3i b2ip+b3i w 2i = const i (10)
The inverse of w 1i and w 3i are an upper bound on the desired sensitivity loop shape and complementary sensitivity function T (s), the inverse w 2i will effectively limits the controller u [START_REF] Skogestad | Multivariable feedback control :analysis and design[END_REF].

C. Parameters optimization

The GA approach [START_REF] Haupt | Practical genetic algorithms[END_REF] is used to compute the parameter weighting functions namely: a 1i , a 2i , a 3i , b 1i , b 2i , b 3i and const i , and to obtain an optimal H ∞ controller. The concept is to generate a population of chromosomes representing the parameters to optimize, and subjected to check whether it satisfies the performance index or not. If any the chromosomes does not satisfy the performance index a new chromosome is generated repeatedly until satisfied. The controller synthesis is realized with the following computing procedure form (see Algorithm 2):

For displacement on x-direction, the computing procedure allows to obtain the following weighting function:

W 11 = s+1 s+0.9258 W 31 = s+0.005935 0.9613s+0.9032 W 21 = 0.7097
The optimized discrete transfer function K(z) is given by: Fig. 6. H∞ structure Algorithm 2 GA approach 1: Parameters initialization: Convert binary population to a decimal population:

a 1i ∈ [0 1], a 2i ∈ [10 100], a 3i ∈ [0 1],b 1i ∈ [0.001 0.01], b 2i ∈ [0 1, b 3i ∈ [0 1],const i ∈ [0.
q(r) = [max(r) -min(r)] /2 N -1 , 1 ≤ r ≤ N 9:
Evaluate the individual fitness function of all individuals: fori ← 1, N p do f i = T (s) ∞ end for (where f i =array of fitness function of chromosomes) the chromosome i is reproduced integer part of e j 15:

Breed new individuals through crossover and aleatory mutation operations of genes to obtain new chromosomes After obtain optimal weights, use Riccati approach to synthesize H ∞ controller and use a reduction method based on singular values to reduce the controller order. 18: until last generation K x (z) = 1.186z 4 -4.01z 3 +5.623z 2 -3.953z+1.154 z 4 -3.25z 3 +4.405z 2 -3.06z+0.9044 then, we obtain u x (k).

Fig. 7 represents the fitness evolution for each chromosome of the population by generation. As one can see, Fig. 7(a) demonstrates a rapid convergence of the fitness. For each generation is selected the chromosome population having a maximum fitness as shown in Fig. 7(b). Once again, a rapid convergence of the fitness parameter to an optimal value γ opt = 1.01 is observed, which assumes that many solutions of parameters are envisaged. As GA exhibits an excellent characteristic of global search and selection technique, the chromosomes converge to the same optimal value. Fig. 8 represents the W 31 parameters optimization evolution and illustrates the parameters optimization convergence.

V. EXPERIMENTAL RESULTS

A. Particle filter validation

To evaluate the performances of the proposed PF, used to estimate the laser beam position, we controlled the laser beam trajectory following a synthetic trajectory generated randomly without working zone. For particle filtering, the particles number is determined such as the computational remained suitable for real-time applications. After several CONFIDENTIAL. Limited circulation. For review only. The results of the laser beam motion prediction using the PF are presented in Fig. 9. At first glance, the filter succeed to follow the true trajectories very closely. Performance tests were realized by comparing the performances achieved by the PF with a classical Kalman filter. The evolution of the mean quadratic error ε(n) is represented in Fig. 10. The results exhibit that the performances of the PF are better than the classical Kalman filter in terms of position precision. Once again, Fig. 11 illustrates that the PF commits less error on the estimation of the beam velocity.

B. Nanomanipulation Validation

The control strategy presented in section IV has been implemented in the experimental platform described in Section II. Fig. 12 shows the microhandling task of a microsphere with a diameter of 20 µm in the focus on the laser beam. First, we positioned the laser beam in the sensed region of the four-quadrant PSD sensor. Once the localization and the position estimation of the laser beam intensity sensed by the PSD, the particle filter allowed a maximum local search of the beam intensity reflected on the microshpere and projected on the PSD surface. Fig. 12(a) shows the detection of one AFM tip by the laser beam, and illustrates the robustness of the PF to determine the localization of the maximum intensity. Then, the gripping and the transport of the sample

VI. CONCLUSIONS

The work presented in this paper deals with the focalization problem of a microsample in the laser beam for radiation applications. The proposed control methodology is based on a two-fingered AFM-based nanogripper handling a microsample illuminated by a laser beam. By processing the maximum beam intensity sensed by a four-quadrant photodiode, the position of the laser beam is detected and located. The main idea of the controller is to guide and to track the beam reflected on the PSD sensor with the cooperative control commands sent to the hybrid micro/nanomanipulators. The experimental results demonstrate the tracking performances when the incident laser beam spot is subjected to external disturbances such as random noises. We used stochastic filters (particle filter and Kalman filter) to localize the laser beam in a three-dimensional space. The preliminary results show that the stability of the microhandling task is preserved during the laser tracking operation.

Fig. 2 .

 2 Fig. 2. Experimental setup.

Fig. 3 .

 3 Fig. 3. Schematic diagram of the architecture of the laser beam tracking control system.

Fig. 4 .

 4 Fig. 4. Output voltage curve Vx and Vs of the four-quadrant PSD.

  k , y k ) and ( ẋk , ẏk ) are the source potion in the plane x-y and velocity respectively.

Fig. 5 .

 5 Fig. 5. Master-slaver controller with decoupling structure for maximum light tracking.

  CONFIDENTIAL. Limited circulation. For review only. Preprint submitted to 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems. Received February 5, 2014.

01 1] 2 :

 2 Parameters and generations number N = 21 G n = 100. 3: Bits and chromosomes number N b = 5 N p = 7. 4: Population size P z = (N p , N • N b ). 5: Crossover rates=0.75 and mutation rates=0.3. 6: Generate randomly the initial binary population of weights parameters. 7: repeat ⊲ for each generation : 8:

10 :

 10 Apply the Elitism selection technique to select the best-fit individuals for reproduction11: compute the means fitness:F = Np i=1 f i , F m = fi F 12:for each chromosome, compute: 13: forj ← 1, N p do e j = fj Fm end for 14:

16 :

 16 Replace least-fit population with new individuals17:

  Preprint submitted to 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems. Received February 5, 2014.

Fig. 7 .

 7 Fig. 7. Fitness evolution: (a): Population fitness computed for each generation, (b): Optimal fitness

Fig. 8 .Fig. 9 .

 89 Fig. 8. w 31 parameters evolution for x axis displacement

Fig. 10 .

 10 Fig. 10. Mean quadratic position error

Fig. 12 .

 12 Fig. 12. Microhandling task of a 20 µm spherical bead using a dualnanomanipulation system: (a) detection of the AFM tip#1 by the laser beam, (b) approach of the AFM tip#1 and AFM tip#2 close to the microobject and (c) stable and reliable handling of the microobject.

Preprint submitted to 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems. Received February 5, 2014.


Zoubeyr Boudene

Adrien Boudier

Pierre-Alain R Breuil

Hélène Olivier-Bourbigou

Pascal Raybaud

Hervé Toulhoat

Theodorus De Bruin
email: theodorus.de-bruin@ifpen.fr

Understanding the Role of Aluminum-Based Activators in Single Site Iron Catalysts for Ethylene Oligomerization

Keywords: Homogeneous catalysis, oligomerization, ethylene, catalyst activation, cocatalyst, methylaluminoxane, mechanism, DFT, molecular modeling

In a combined experimental and theoretical study the activation process of a single site ethylene oligomerization catalyst with aluminum based activators has been studied.

The results put forward a plausible deactivation reaction path of the catalyst for trimethylaluminum, while for methylaluminoxane and a novel phenoxyaluminum-based activator the experimental catalyst's activity correlates with the energy barrier for the ethylene insertion.

Introduction

Methylaluminoxanes (MAO) [1][START_REF] Kaminsky | Process for Producing Polymers and Copolymers of Ethylene[END_REF][3][4] have been widely used as activators of ethylene oligomerization and polymerization iron-based catalysts, especially those belonging to the bis(imino)pyridine family, affording highly productive catalytic systems [5][6][7][8][9][10][11][12]. As high MAO/metal ratios are required, leading to high costs, considerable efforts have been undertaken to replace them in oligomerization reactions [13][14][15]. This requires a good understanding of the activation process of the catalyst. However, the identification of the precise role and active center of MAO during the activation process is hampered by the multitude of aluminum species present in solution [16,17], the control of the reaction exothermicity, along with the catalyst degradation and polymer byproduct formation.

In this communication we report a well-defined dimeric aluminum complex [PhOAlMe 2 ] 2 and its activity for ethylene oligomerization using ({2,6-(2-(CH 3 )C 6 H 4 N=C(CH 3 )) 2 -C 5 H 3 N}Fe(II)Cl 2 ) (abbreviated as LFeCl 2 ) precursor in comparison with the MAO and trimethylaluminium (TMA) cocatalysts. We furthermore present the DFT simulations that helping to explain the observed experimental activities. [PhOAlMe 2 ] 2 and TMA are structurally well-defined aluminum species, whereas MAO is not. For the latter, a previously described model is therefore used (Scheme 1, structure (iii) in the lower rectangular box) [18][19][20].

While the iron-catalyzed oligomerization mechanism is largely accepted to follow a Cossee-Arlman type mechanism, leading to a Schulz-Flory distribution of linear -olefins, the debate concerning the nature of the metallic active species and the corresponding activation path remains open and is fuelled with numerous studies regarding the oxidation state of the active species or the non-innocence of the ligand towards redox reactions [21][22][23][24][25][26]. For example, [LFe(II)Me] + was observed by mass spectrometry along with other cationic species after activation of [LFe(II)Cl 2 ] by MAO. We here assume a high-spin ferrous iron cationic species after activation with a redox-innocent neutral bis(imino)pyridine chelate and suppose no electron transfers between the metal and ligand.

Al

R X Cl N N N Fe Al R X Cl N N N Fe N N N Al R X Cl Fe N N N Al R X Cl N N N Fe Al R X Cl N N N Fe N N N Fe N N N Fe Al Z Al R Cl Z Al R X Cl Fe 1 2 1 2 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 1 2 C 2 H 4 3 4 + Al R X Cl + Z= -Me or -OPh (I) (II) (II-a) (III-a) (III-b) (III-c) (II-b) (III-d)

Materials and Methods

Methylaluminoxane (MAO 10 wt % in toluene) and TMA (trimethylaluminum) were purchased from Chemtura, while [PhOAlMe 2 ] 2 was obtained by addition of an equimolar amount of phenol to TMA in n-heptane in good yield (97 %) and was characterized by 1 H and 13 C NMR and XRD [27]. [LFe(II)Cl 2 ] was prepared according to the procedure described in ref [28].

A toluene solution containing [LFe(II)Cl 2 ] (10 μmol) was injected in the autoclave under an inert nitrogen atmosphere, followed by the cocatalyst solution. The ethylene pressure was immediately increased to 30 bar, the temperature to 50 °C and the mixture was mechanically stirred. Additional details are provided in the Supplementary Materials.

DFT calculations were carried out with the unrestricted M06 functional [29], proven to provide an accurate description of aluminoxanes [30] as well as the spin state of the iron catalyst (Table 5S), in combination with Los Alamos pseudopotential and associated valence basis set on iron atoms [31] and the basis set 6-31G(d,p) for the other atoms. We have verified that inclusion of solvent effects or expansion of the wave-functions over larger basis sets do not change the conclusions presented below from results of gas phase DFT calculations. This is in line with the theoretical findings of Zurek and Ziegler who also concluded that relative energies were not impacted by the inclusion of a solvent model [19]. Further details are presented in the supplementary materials.

Results and discussion

Zero oligomerization activity was observed with the [LFe(II)Cl 2 ]/TMA system, even with an Al/Fe molar ratio of 500, while at the same ratio an activity up to 1.4×10 5 g•(mol(Fe)•h) -1 was measured with [PhOAlMe 2 ] 2 . Ideally, lower Al to Fe ratios are sought after, as put forward in the introduction, but this leads to lower oligomerization activities.

A Schulz-Flory distribution (K = 0.68) of linear -olefins was afforded comparable to the one observed when the activation is performed with MAO (Al/Fe = 250, K = 0.69). This latter catalytic system showed a higher, up to 7.4×10 7 g•(mol(Fe)•h) -1 , activity. Nevertheless, the iron bis(imino)pyridine/MAO system is temperature-sensitive and coupled to the important exothermicity of the reaction, processing issues are encountered even at low catalyst loading. In our case, more than 40 wt% of the formed products are waxes and polyethylene (PE) with this system, while none of these heavy products were obtained with [PhOAlMe 2 ] 2 .

As [PhOAlMe 2 ] 2 is a structurally well-defined cocatalyst, it serves as a very good candidate in this theoretical study to help to elucidate the activation mechanism of the precatalyst using MAO.

[LFe(II)Me] + may be obtained from [LFe(II)Cl 2 ] after a transfer of the two chloride anions to the cocatalyst, which becomes either mono-or dichlorinated, while Fe(II) gets mono-methylated.

These initial steps have been studied with DFT calculations for the TMA and [PhOAlMe 2 ] 2 cocatalyts. Only minor energy differences between the reaction intermediates were found for both systems, thus not explaining the experimental activity difference. Since moreover, the complete dissociation of the created cationic iron species and the anionic cocatalyst is prohibitively costly, we will focus on the reactions of the ion pair, i.e., cationic monoalkyl [LFe(II)Me] + and investigate its interactions with the three cocatalysts, each being mono-or dichlorinated (Figure 5S) as well as with ethylene.

Table 1 shows that only the formation of the monochlorinated TMA anion is an exergonic reaction. For [PhOAlMe 2 ] 2 and the MAO model chloride anion formation is endergonic and moreover the dichloride anion is thermodynamically slightly more stable. Since the energy difference is small between the mono-and dichlorinated anion both intermediates will be considered in this study. This endergonicity puts forward an additional reason 1 The ethylene uptake by (I) yields intermediate (II), where the coordinated ethylene can be inserted in the Fe-C 1 bond to give (III-a). The intermediates and transition states (II-a ‡ ) for the insertion of ethylene into the Fe-C bond were optimized in the presence of the different chlorinated anionic activators (Figure 6S). For each ion-pair, an equatorial and axial approach (Figure 4S) of the ethylene molecule has been evaluated and only the most stable -complex was 1 Other evoked reasons are the partial "loss" of MAO due to reactions with impurities (scavenger effect, ref. [4]) and minor presence of reactive aluminum sites in MAO oligomers, as shown in ref [30].

retained. The calculated Gibbs energies for different optimized intermediates and transition states are reported in Figure 1. The ethylene uptake is endergonic for all ion pairs obtained from the activation with the three different mono-or dichlorinated anions, note that the reported energies are relative to the energies of the ion pair and ethylene molecule at infinite distance. Considering the monochlorinated counter ions, the Gibbs energy barriers for the ethylene insertion are 16.8, 24.9 and 29.4 kcal•mol -1 for respectively [AlMeO] 6 [AlMe 3 ], [PhOAlMe 2 ] 2 and TMA. Interestingly, the energy barrier for the insertion with the MAO model is even slightly lower than the value obtained for ethylene insertion inhibition path 1 the bare cationic -complex (DG ‡ = +19.1 kcal•mol -1 ). However, the relative energy is lower for the bare complex.

The energy barrier for ethylene insertion is influenced by different parameters, like the interaction of the -complex with the counterion. Hence, the smaller dissociation energy of the anion, i.e. 60.7, 67.8 and 81.0 kcal•mol -1 for respectively, MAO, [PhOAlMe 2 ] 2 and TMA, the lower the barrier. Also the change in the C-C ethylene distance in going from the -adduct to the transition state, convincingly correlates with the energy barrier: the smallest changes are observed for the smallest barrier energies, where the ethylene is weakly coordinated, e.g. in the bare complex and in presence of the MAO (see Table 3S and Figure 2S).

The effect of the presence of dichlorinated anions on the uptake of ethylene and the insertion path has also been calculated. Especially, the formation of the -adduct in presence of the dichlorinated phenoxy species, [Al 2 Ph 2 O 2 Me 3 Cl 2 ]¯, becomes disfavored: DG r = +13.7 kcal•mol -1 (Figure 1S) with respect to the monochlorinated form (DG r = +13.7 kcal•mol -1 ). With a barrier of 20.7 kcal•mol -1 for the ethylene insertion, the TS consequently gets as well significantly destabilized, making the overall insertion step slightly endergonic (DG r = +2.0 kcal•mol -1 ). Since [PhOAlMe 2 ] 2 experimentally does activate [LFe(II)Cl 2 ], we conclude that the monochlorinated anion is involved in the activation process, as it exhibits the more favorable (kinetic and thermodynamic) pathway. The dichlorinated form of TMA and MAO influences considerably less the stability of the -complex and the insertion barrier, as can be seen from Figure 1S.

The free energy barrier for ethylene insertion is lower in case of [AlMeO] 6 [AlMe 3 ] than for [PhOAlMe 2 ] 2 or TMA, thereby providing an explanation for the higher experimental oligomerization activities. However, the calculated activation energy barriers do not explain the experimental zero activity with TMA. Therefore, we need to search for another origin and explore routes for the catalyst's inhibition. respectively. Note that energies are relative to the energy of the ion-pair and the ethylene molecule at infinite distance. Hence, this inhibition reaction is unlikely to happen in the presence of ethylene.

We then investigated an alternative inhibition path that could take place before the ethylene uptake, e.g., the iron methylation: reaction (I) to (III-c) in Scheme 1. However, the calculations reveal that product (III-c) is unstable for all three monochlorinated anions and spontaneously dissociates to (I).

A third inhibition path, (I) → (III-d) in Scheme 1 has been considered. The chlorinated dimeric TMA anion is already partially dissociated in the presence of the cationic iron species. The complete dissociation of TMA (Eq. 1) requires DG diss = +4.5 kcal•mol -1 for the monochlorinated TMA dimer (R=Me) and even is exothermic for the dichlorinated anion (R=Cl): DG diss = -7.5

kcal.mol -1 .

[

LFeMe] + ••• [AlMe 3 ][AlMe 2 RCl]¯ → AlMe 3 + [LFeMe] + •••AlMe 2 RCl]¯ (1)
The subsequently generated AlMe ¯. This coordination of AlMe 3 could therefore be responsible for the poisoning of the catalyst. Moreover, ethylene is not sufficiently nucleophilic to replace the TMA monomer, as was concluded from a relaxed coordinate scan 2 , where the interchange of TMA with ethylene costs more than 18 kcal•mol -1 . These results seem to contradict the work of Bryliakov et al [25].

Nonetheless, they studied the iron 2,6-bis(2,6-diisopropylphenylimino)ethyl-pyridine catalyst, which experimentally can be activated by TMA, whereas we focused on the 2,6-bis(2methylphenylimino) analogue, which is experimentally inactive with TMA. Accordingly, our calculations explain these experimental findings by the fact support that the TMA association (eq.

2) is exergonic (-4.7 kcal•mol -1 ) for the bis(methylphenyl) complex, while it is endergonic (+5.6 kcal•mol -1 ) for the bis(diisopropyl) complex, likely due to the steric hindrance of the bulky isopropyl groups.

[

Al 2 Me 6 Cl)]¯•••[LFeMe] + → [AlMe 3 Cl)]¯•••[LFeMe.AlMe 3 ] + (2)
The dissociation reaction (eq. 1) with [PhOAlMe 2 ][PhOAlClRMe 2 ] is unlikely to occur since the complete dissociation costs +15.3 kcal•mol -1 for the monochlorinated (R=Me) and +19.9 kcal•mol -1 for the dichlorinated anion (R=Cl). Moreover, in the former case a relaxed coordinate analysis, where the Al-Al distance is stepwise increased, shows a barrier of about 25 kcal•mol -1 2 In this relaxed coordinate scan the r(Fe-Al) distance is constrained but systematically increased in steps of 0.5 Å, while all other atoms are relaxed (Figure 3S). Even at r(Fe-Al) = 10 Å the scan shows a tendency of increasing energy, taking the fully relaxed structure as a reference. We were not able to locate the transition state on the potential energy surface.

to generate the mononuclear Lewis acid Me 2 AlOPh that could deactivate the catalyst. 3 Since this inhibition reaction is in direct competition with the ethylene insertion (24.9 kcal•mol -1 energy barrier, but exergonic (DG r ins = -10.9 kcal•mol -1 ) the insertion reaction clearly is favored. We can therefore reasonably assume that this inhibition route is unlikely to take place. In the case of the MAO all aluminum atoms are tetravalent, thus lacking the possibility to coordinate to methyl on the iron center, making this inhibition reaction improbable.

Conclusions

In summary, we have identified [PhOAlMe 2 ] 2 as cocatalyst able to activate iron bis(imino)pyridine precatalysts. This well-defined cocatalyst allowed us to unravel the ethylene insertion and inhibition mechanisms of the active species by DFT calculations. These key mechanisms were further extended to a representative MAO model and TMA and furnish a quantified analysis of the experimental ranking of the three cocatalysts: MAO > [PhOAlMe 2 ] 2 >> TMA.

We conclude that efficient cocatalysts should be able to generate a cationic alkylated iron species, whereby the cocatalyst itself is turned into a chlorinated anion. It seems that this anion then should weakly interact, since lower barriers are observed for anions with smaller cationic complex/anion dissociation energies and having higher activities. Additionally, the cocatalyst should not contain (or easily generate) trivalent reactive aluminum sites (such as in TMA) that could form adducts with the already generated cationic methylated iron species, thereby 3 We were unable to locate the transition state on the potential energy surface.

  we first show how the nature of the anion influences the ethylene insertion mechanism: (II)  (III-a) in Scheme 1. In parallel, we explore possible routes leading to the catalyst deactivation: (II)  (III-b), (II)  (III-c) and (II)  (III-d).

Figure 1 .

 1 Figure 1. Gibbs energy profiles (kcal•mol -1 ) for the ethylene insertion step (right) and inhibition

  The inhibition reaction path (II) to (III-b) in Scheme 1 involving ethylene, was studied for the three monochlorinated anions (Figure2, left). The activation energies show that this path is kinetically disfavored with respect to the oligomerization: the difference between the two barriers (DG ‡ b -DG ‡ a ) is +6.9, +9.0 and +23.4 kcal•mol -1 . Moreover, all reactions are endergonic: 9.1, 18.0 and 21.6 kcal•mol -1 for TMA, [PhOAlMe 2 ] 2 and [AlMeO] 6 [AlMe 3 ]

  3 easily coordinates to the [LFe(II)Me] + in presence of the anionic [AlMe 2 RCl] ¯ species by forming two bridging methyl groups between the Al and iron center (the methyl group of [LFe(II)Me] + and one of the AlM 3 methyl groups, Figure 7S) forming III-d, thereby consuming the vacant site that is required for ethylene coordination. The overall reaction is exergonic by -4.7 in the presence [AlMe 3 Cl]¯ and -8.6 kcal.mol -1 in the presence of [AlMe 2 Cl 2 ]

(i) (ii) (iii) Scheme 1. Pathways for ethylene insertion and catalyst inhibition.

			insertion
							Inhibit.1					Inhibit.2
														Inhibit. 3
							R= Me or Cl
	Me Me	Al	O O Ph Ph	Al = X Me Me	Me Me	Al	Me Me = X Al	Me Me	Me Me Al O	O Al	O Al Me	Al O	Al O Al O Me Me Me	= X Me Al Me

Table 1 .

 1 why such large excesses of [PhOAlMe 2 ] 2 and MAO are needed to activate. Calculated Gibbs energies (kcal•mol -1 ) to form the mono-and dichlorinated anionicIn the following we assume that each cocatalyst is able to generate the [LFe(II)Me] + species;

	cocatalyst.
		[PhOAlMe 2 ] 2	TMA dimer [AlMeO] 6 [AlMe 3 ]
	Monochlorinated	13.2	-0.6	5.9
	Dichlorinated	11.6	6.5	2.4

poisoning the catalyst. We hope that these results help for a more rational understanding of the role of cocatalysts in oligomerization reactions and may open new routes for the molecular design of improved cocatalysts.

Supporting Information.

Computational details, cocatalyst synthesis, catalytic experiments, geometrical features of intermediates and transition states, and Cartesian coordinates of all fully-optimized geometries.


Nathalie Courtois
email: n.courtois@brgm.fr

Ariane Grisey
email: grisey@ctifl.fr

Dominique Grasselly

André Menjoz

Yves Noël

Vincent Petit

Dominique Thiéry

Application of Aquifer Thermal Energy Storage for heating and cooling of greenhouses in France: a pre-feasibility study

Keywords: Heating and cooling, Aquifer Thermal Energy Storage (ATES), greenhouses, numerical modeling

In the present context of predicted shortage of fossil energies and awareness of environmental preoccupations (greenhouse gases for example), saving energy and converting it to renewable one become priorities. As energy-intensive consumers, greenhouses are directly concerned. That's why Ctifl (Technical Institute for Fruit and Vegetables) would like to develop the concept of "sustainable greenhouses" in France, using renewable energy. In this case, the method which retains the attention is the concept of reversible heating and cooling through Aquifer Thermal Energy Storage (ATES). This technique is already in application in other European countries, especially in the Netherlands, but has never been tested in France for greenhouses. This paper presents the research project actually led by Ctifl and Brgm (French Geological Survey). The objective of this project is to determine the parameters to take into account (and their relative importance), to evaluate the pre-feasibility of this technique on an agricultural site. This was done through numerical modeling of theoretical aquifer. The model used is a finite volumes software developed by Brgm, which can treat both water flow and thermal transfers. A sensitivity analysis is led on several sets of parameters sets, which depend either on exploitation conditions (pumping and injection discharge, distance between wells, etc.), or on aquifer conditions (geometry, thermal and hydraulic characteristics, etc.).

INTRODUCTION

In the now urging context of predicted shortage of fossil energies, and with the increasing awareness of environmental preoccupations (greenhouses gas for example), saving energies and converting to renewable ones become priorities. As energy-intensive consumers, greenhouses are directly concerned. In the past few years, the price of fossil energies as gas and fuel has increased in such a way that the annual costs of heating become a larger and larger part of the total exploitation charges of an agricultural site (about 20 to 35%).

That's why the Ctifl (Technical Institute for Fruit and Vegetables) would like to develop the concept of "sustainable greenhouses" in France, promoting the use of renewable energy. In this case, the method which retains the attention is the concept of reversible heating and cooling through Aquifer Thermal Energy Storage (ATES). This technique is already in application in other European countries, especially in the Netherlands. A status beginning of 2005 reports that over 400 projects were operational in the Netherlands, concerning office and commercial buildings, hospitals, housing, industry and agriculture [START_REF] Snijders | Aquifer Thermal Energy Storage in the Netherlands, Status beginning of 2005[END_REF]. Other examples of applications can also be mentioned, as the new parliament building in Berlin (Germany), the Sussex hospital in Canada, large scale experiences in Sweden (heating and cooling of commercial buildings) and an experimental greenhouse in Turkey [START_REF] Turgut | Aquifer thermal energy storage application in greenhouse climatization[END_REF]).

This paper deals with a pre-feasibility study of Aquifer Thermal Energy Storage, especially in relatively few deep aquifers (10m to 100-150m depth), applied for heating and cooling of greenhouses in France. This study is actually led by Brgm and Ctifl. Its objective is to determine the parameters to take into account (and their relative importance), to evaluate the pre-feasibility of this technique on an agricultural site, for a given range of energetic needs, both for cooling and heating. This study is led through numerical modeling of "theoretical aquifers". Different sets of parameters, depending both on aquifer characteristics and on exploitation conditions, are tested and compared through their effects on ATES efficiency.

USING GREENHOUSE AS "SOLAR CAPTOR" -DIMENSIONING THE NECESSARY FLOW RATE FROM AQUIFER

With large glass surfaces, greenhouses act as real "solar captors". The concept of "solar greenhouse" consists in exploiting heat surplus for heating in winter by storing this surplus in aquifer. The objective is to calculate the net energy surplus which is the result of energy inputs (solar radiation) minus the heat losses. The energy balance allows then to deduce the flow rates from the aquifer. Table 1 shows an example of energy balance for climate conditions in the South of France, with a temperature set point of 17 °C in the greenhouse during the night. The maximum cooling capacity is around 500 W/m² to have a temperature in greenhouse below 28°C during the summer period. The first calculation shows that the maximum flow rate demand from the aquifer is around 260 m 3 /h/ha during the cooling period.

PRINCIPLE OF AQUIFER THERMAL ENERGY STORAGE (ATES) USING A REVERSIBLE GEOTHERMAL DOUBLET

The principle of Aquifer Thermal Energy Storage is to take advantage of the thermal capacity of both the geological formations and the water they contain: groundwater is used both as a reservoir and a vector of energy. Geological materials constitute favourable environment for energy storage as they present low thermal conductivities leading to a slow diffusion of energy and moderate thermal losses [START_REF] Gauthier | Evaluation du potentiel géothermique des nappes aquifères avec stockage thermique : revue de littérature et application en deux contextes du Québec[END_REF]). The environments of consideration are mainly quite shallow aquifers, lying in few ten meters depth, and where the groundwater temperature remains quite constant over the year (close to the annual mean temperature of the outside air at the site). This low temperature (<30°C) geothermal energy is largely exploited for the heating and cooling of houses, through the use of a geothermal pump (as in the Ile-de-France region for example).

The Aquifer Thermal Energy Storage installation envisaged in the study is based on the use of a geothermal doublet, that means a pair of water wells, one "hot well", and one "cold well". The system is said reversible, as each borehole is used alternatively in pumping or injection according to the season. The principle scheme of a reversible geothermal doublet is shown in Figure 1:

in cold season, groundwater is pumped from the "warm well", gets cold while heating the greenhouse, and is injected in the "cold well", -in the hot season the flow is inverted: groundwater is pumped from the "cold well", gets heated while cooling the greenhouse, and is injected in the "warm well".

This use of a reversible geothermal doublet presents some advantages:

storage of heated and cold water increases the temperature contrasts, and so the efficiency of the system, -pumped groundwater is re-injected in the same aquifer, limiting the risks of hydraulic overexploitation, -system is inverted seasonally, limiting the risks of longterm warming or cooling of the aquifer, that could cause both environmental degradations and a drop in efficiency [START_REF] Bridger | Designing Aquifer Thermal Energy Storage Systems[END_REF]).

This functioning involves also some counterparts:

the investments costs are higher than in a non reversible system, as boreholes have to be very carefully designed to function in pumping as well as in injection, -injection in a borehole is always more difficult than pumping (higher pressure necessary to inject than to pump the same water discharge, clogging risks, etc.). This implies an over-dimensioning of the boreholes (diameters, strainer, etc.) compared to wells acting only for injection.

It should be noted that a particular attention would have to be taken on all the surface installations to avoid any contaminant risks of the re-injected groundwater.

PROCESSES INVOLVED IN ATES SYSTEM

Thermal energy transfer in the aquifers

The thermal energy transfer in the aquifers is governed by thermal diffusion, advection, and dispersion. These processes are summarized in Figure 2. Diffusion is the transfer due to thermal conductivity of the aquifer and its upper and lower limits. Diffusion in one direction is proportional to the thermal conductivity coefficient λ a (expressed in W/m/°C) and to the temperature gradient in that direction. The thermal conductivity of an aquifer depends both on the type of geological material (sand, gravel, etc.), and on its saturation degree by water (or porosity ε. for confined aquifers). Current values of thermal conductivity of aquifers are given in the range 2.0-2.5 W/m/°C, with porosity between 5% and 20%.

Aquifer acts also with a thermal inertia due to its calorific capacity γ a which represents the heat quantity necessary to increase the temperature of a 1 m 3 volume by 1 °C. Current values of aquifers calorific capacity are given in the litterature in the range: 2.0-2.9.10 6 J/m 3 /°C.

Most of the aquifers are subject to a regional natural flow V which depends both:

on the aquifer hydraulic conductivity K (in m/s), expressing its aptitude to let the water goes through under the effect of a hydraulic gradient i, -on the value of the hydraulic gradient i, which represents the slope of the piezometric surface of the aquifer, and acts as a potential gradient.

i K V . = (1)
where V, K, i are aquifer Darcy's velocity, hydraulic conductivity, and hydraulic gradient, respectively.

The regional flow V can be very low (few meters per year), to relatively large (few meters per day). Current values of hydraulic gradient are usually in the range 0.1 -10 ‰.

Advection represents the moving of the thermal stock due to the natural flow of the aquifer. To this ensemble movement adds a spreading out of the thermal stock due to spatial heterogeneity of the velocity field. This phenomenon is called dispersion, and leads to an increase of the "global" aquifer thermal conductivity as follows:

V V f T a global T f L a global L ⋅ ⋅ + = ⋅ ⋅ + = γ α λ λ γ α λ λ (2)
where λ a , γ a , γ f , V, α L , α T are aquifer thermal conductivity, aquifer calorific capacity, fluid calorific capacity, regional flow of the aquifer (Darcy's velocity), aquifer longitudinal (in the mean flow direction) and transverse (transverse to the mean flow direction) dispersivities, respectively.

Factors influencing the efficiency of ATES

As shown in Figure 2, different processes control the thermal energy transfer processes in aquifers, and will influence the efficiency of the thermal energy storage:

ensemble movement and spreading out of the "cold" and "warm" stock around the wells depending of the aquifer velocity field, -thermal inertia depending on the aquifer calorific capacity, -thermal losses trough the upper limit depending on the temperature gap between outside air and groundwater, and on the thickness of the cover.

Other factors concerning the exploitation conditions that have to be dimensioned according to the aquifer characteristics (and to other considerations as energetic needs, costs, available space, etc.) to optimise the ATES efficiency:

pumping and injection discharge Q, -pumping and injection cycles (duration, rest periods, etc.), -injection temperatures in the cold and in the hot wells, -distance L between the hot and the cold wells (to avoid interference between the heated and cold stocks), -disposition of the axis of the geothermal doublet towards aquifer flow direction.

SENSITIVITY ANALYSIS THROUGH NUMERICAL MODELLING

The objective of the pre-feasibility study was to evaluate and prioritize the effects of different factors that can influence the feasibility and the efficiency of ATES. Some theoretical studies have already been done to determine the influence of physical parameters, and have led to the edition of graphs for dimensionless variables [START_REF] Sauty | Du comportement thermique des réservoirs aquifères exploités pour le stockage d'eau chaude ou la géothermie basse enthalpie[END_REF]). If these curves can be used to obtain a quick evaluation of feasibility of ATES for a given configuration, they are nevertheless of limited help as they are generally based on oversimplified hypothesis (simple geometry of aquifer, etc.). Numerical modeling is then essential to take into account more complex aquifer geometries, density effects, etc., and to evaluate the evolution of heated and cold groundwater stocks in space and time.

Model definition

Grid

A numerical study was led with the MARTHE model, a finite-volumes software developed by Brgm that can model both hydrodynamics and thermal transfers.

The model is composed of 10000 cells with varying size, as shown in Figure 3. A finer nested grid is included close to the well to obtain a better definition of the thermal storage.

Figure 3: Above view of the model grid

In the simulations presented below, the axis of the geothermal doublet was supposed to be transversal to the regional aquifer flow.

In the vertical direction, the terrain is subdivided into 17 horizontal layers with varying thickness. As shown in Figure 4, outside air temperature is prescribed on the first layer of the cover. The aquifer part is represented by three layers. The vertical water movements due to density phenomenon have been considered as negligible for the considered range of temperature (from 10 to 30°C). A pumping/injection schedule has also been defined over a one year cycle:

cooling: pumping of cold water and injection of heated water from June to September (summer), -rest (no pumping) in October / November, -heating: pumping of heated water and injection of cold water from December to March (winter), -rest (no pumping) in April / May. These exploitation conditions remained the same for all numerical simulations.

Modeling parameters

Simulations were realized with a weekly time step, for 15 yearly cycles (180 months) following the exploitation schedule presented above. The computer CPU time for 15 years of exploitation was about 3 hours on a standard desktop PC.

The initial temperature of the aquifer is 14.5°C (annual mean air temperature). The "basic" principle of a sensitivity analysis is to make the different parameters varying one by one and to compare the obtained simulations results. In this study, several methods have been retained to compare the simulations, in order to evaluate the relative influence of the different parameters on ATES efficiency.

The "global" thermal power that can be furnished by water is proportional to water discharge and temperature following the following equation:

T Q P f g ⋅ ⋅ =γ
(3) where P g , γ f , Q, T are "global" thermal power, fluid calorific capacity, discharge and temperature, respectively.

The model gives the time evolution of temperatures simulated in the cells containing respectively the hot and the cold wells. As these temperatures will condition the thermal power, the simulations are compared through:

the temperatures simulated at the end of the 4 months pumping cycle, -the deviation with the natural aquifer temperature, -the number of yearly cycles necessary to obtain the stabilization of the temperatures obtained in the hot and cold wells at the end of the 4 months pumping cycle.

The MARTHE model allows also the simulation of the spatial distribution of the temperatures in the aquifer. For the different simulations, these are compared for some key dates of the yearly exploitation cycle (end of pumping cycle, end of rest period, etc.).

In this particular case of thermal storage, the thermal efficiency can be measured through the deviation between the temperature of the stored (and pumped) groundwater and the initial and natural aquifer temperature. We introduce here the notion of "useful" thermal power, that allows to evaluate the surplus of power obtained by storage compared to a "simple" exploitation of groundwater (at constant temperature) without storage.

T T

T P T T Q P g f u 0 0 - ⋅ = - ⋅ ⋅ =γ (4)
where P u , P g , Q, T, T 0 are global and useful thermal power, groundwater discharge, groundwater temperature at wells, and initial aquifer temperature, respectively.

To compare storage efficiency, a recovery factor is also calculated, for the cold and the hot wells respectively, as the ratio between the quantity of energy pumped during a 4 months cycle, and this injected during the previous season:

stored u pumped u E E r = (5)
where r, , are recovery factor, pumped and stored "useful" thermal quantity of energy, respectively.

pumped u E stored u

E

The "useful" quantity of energy is defined as the integral of the "useful" thermal power on a 4 months cycle:

∫ = months u u dt t P E 4 0 ). ( (6)
where E u , and P u are "useful" thermal energy and thermal power, respectively. For all the simulations, this recovery factor is calculated for the 15 th year of exploitation.

Simulation results

Influence of aquifer regional flow

Aquifer regional flow is one of the most important parameters that will condition the efficiency of ATES. It will depend on the permeability K, and on the hydraulic gradient, i. The graphs below show the comparison of two simulations with and without aquifer regional flow.

Figure 6 shows comparison of cold and heated groundwater stocks for the first and the 15 th years of exploitation, and for different key dates of the yearly cycle. It appears clearly that the thermal storage has increased over the years, leading to an amelioration of the efficiency. The comparison of the first case without regional flow (which, in fact, never happens in reality) and the second one with a flow of 0.86 m/d is eloquent, as it shows the moving and spreading out of the heated and cold water. Figure 5 shows the temperatures simulated in the cells containing the cold and the hot wells. The strong effect of aquifer flow can also be clearly seen, through quicker increase (in the cold well) and decrease (in the hot well), of pumped water temperatures after the end of the 4 months injection period. For the 15 th simulation year, the gain in temperature (compared to the natural aquifer temperature) at the end of the 4 months pumping period is:

in the hot well, +6.3°C and +1.7°C, without and with aquifer flow, respectively, -in the cold well, +0.8°C and +0.4°C, without and with aquifer flow, respectively.

The influence of natural aquifer flow on the ATES efficiency can be evaluated on the curves of "useful" thermal power, shown below (Figure 7). It appears clearly that the recovered thermal powers are strongly affected by the aquifer flow. Compared to the theoretical case without flow, the recovery factors calculated on the 15 th simulation year decrease from 67% to 26% in the hot well, and from 53% to 22% in the cold well. Figure 8 and Figure 9 show that the results are slightly better in terms of recovered temperatures and "useful" thermal powers for a cover 20m thick compared to 5m. 5.2.3 Influence of the distance between the wells Distance between the two boreholes of the geothermal doublet is also an important parameter to take into account.

Figure 10 shows the evolution of simulated temperatures in the hot and cold wells for two distances of 150 and 200m, in a case of a 50 m 3 /h exploitation discharge. The interference between the heated and cold groundwater is greater when the boreholes are closer, this is to be seen almost in the temperatures of the pumped waters at the end of the 4 months pumping cycles. The interference is especially marked on the cold well, as shown on the "useful" thermal power (Figure 11). With a 150m distance, the recovered thermal power in the cold well is hardly decreased by the influence of the hot well, where a greater quantity of energy is stored. In this simulated case with a 50 m 3 /h discharge, we can conclude that it is better to install the boreholes on a distance greater than 150m.

The "minimal" distance between boreholes depends on the thermal radius of the stored water, which depends itself both on the aquifer thickness, and on the water discharge. A given distance could be sufficient for a discharge, and not for a greater one.

In real case, a compromise may have to be found between large distances to avoid reciprocal "thermal pollution" of the two stocks, space availability and disposition constraints on the agricultural site, and additional costs of equipments with the distance (piping lengths, greater hydraulic losses leading to larger pipes diameters, etc.).

Conclusions

The theoretical study led on numerical simulations has shown the variety of factors that influence the efficiency of ATES, dealing with the aquifer characteristics (regional flow, permeability, porosity, thermal parameters), and also with the exploitation conditions (discharge, injection temperatures in the cold and hot wells, pumping and injection cycles, distance between the geothermal doublet).

For this sensitivity analysis, the efficiency of the ATES has been evaluated through:

the temperature deviation between the pumped groundwater and the aquifer, at the end of the 4 months pumping period and for the 15 th year, -the recovery ratio between the recovered and stored quantity of energy on a 4 months cycle.

For the whole simulations at the hot well:

the temperature deviation varies from +7.9°C to -0.4°C, with an average of +2.5°C, -the recovery ratio varies from +74% to -1%, with an average of 36%.

For the whole simulations at the cold well:

the temperature deviation varies from +2°C to -3.2°C, with an average of +0.08°C, -the recovery ratio varies from +66% to -7%, with an average of 26%.

Thermal efficiency is on the whole better in the hot well than in the cold one, due to the fact that the stored quantity of energy is greater, and that, in certain cases, the stock of heated water has a strong influence on the cold one.

This sensitivity analysis allowed to evaluate the relative influence of the different parameters on thermal efficiency. The natural flow of the aquifer is the factor that will mostly condition the efficiency of the ATES. The best recovery ratios are those obtained for the "theoretical" (and never observed in the nature …) cases without flow. The distance between the boreholes is also an important parameter.

This sensitivity analysis also allowed to illustrate the complexity of the phenomena that govern the aquifer thermal energy storage. A same parameter can have positive and negative effects on the storage, or can have a positive effect in a given range of values, and a negative one in an other one, depending on the values of the other parameters (as an example, the increase of the exploitation discharge can produce a positive or a negative effect depending on the aquifer thickness and on the distance between the boreholes). This interdependence between the parameters is clearly shown on this sensitivity analysis, where the recovery ratio varies from negative values to very significant ones (+74%). The conclusion to retain is that each site will be a particular case, requiring a detailed dimensioning as a function of its hydrogeological context.

ATES SYSTEM DESIGN: CASE TO CASE STUDY

The main component of an ATES system is of course the presence, under the agricultural site, of a suitable aquifer. "Suitable" means here able to produce the required water demand for the project. First question is: how many wells might be needed to meet the demand, bearing in mind peak demand, average demand, and also the fact that a well is often less efficient in injection than in pumping [START_REF] Bridger | Designing Aquifer Thermal Energy Storage Systems[END_REF]).

Aquifer characterization for an ATES project requires a case to case study, as hydrogeological conditions can be variable at the local scale. Dimensioning a hydrogeological study for ATES project will hardly depend on the previous knowledge (previous studies, existing nearby water wells, geological maps, etc.), the complexity of the site (type of aquifer, spatial variations of regional flow, etc.), and the size of the ATES system itself. Two phases can be drawn, (i) the prefeasibility study, which is an indispensable preliminary step, and (ii) the detailed design of the whole ATES system. In these two phases, numerical modeling is a precious tool, to evaluate the capacity of thermal storage of the aquifer, and its evolution in space and time.

In the pre-feasibility step, aquifer characterization is almost done through exploitation of existing data, expertise, and eventually low costs investigations on nearby water wells (pumping tests, piezometric measure, etc.). A first set of numerical simulations can be led with parameters issued from previous studies, and literature review, in order to give some ranges of storage capacity and efficiency.

If the pre-feasibility study leads to a priori favorable conclusions, detailed ATES design can be engaged, with complementary aquifer investigations (drilling of test well, geological logging of boreholes cuttings, pumping test, tracing experiments, geochemical analysis, etc.). New sets of numerical simulations can be led with more accurate parameters issued from the local investigations on the site.

CONCLUSIONS AND PERSPECTIVES

Many ATES systems exist all around the world, and have proven to be viable, and energy-efficient technology. Hydrogeological aspects must be considered carefully during the system design to ensure a properly operating system on the long-term (avoid clogging, etc.). To be efficient and viable on the long-term, ATES systems require a really accurate design.

Aquifers at relatively shallow depths (10-100m) are present in large areas in France. ATES could be seen as an interesting perspective for greenhouses farmers willing to reduce their energetic bill and the environmental impacts linked to the consumption of fossil energies. In any case, a case to case study is necessary to evaluate the pre-feasibility on a given site. Numerical modeling appears as a precious tool, from pre-feasibility stage to detailed design.

NOMENCLATURE

Parameter

Figure 1 :

 1 Figure 1: Principle scheme of Aquifer Thermal Energy Storage applied to reversible heating and cooling, through geothermal doublet. Reversible functioning in cold and hot seasons.

Figure 2 :

 2 Figure 2: Principle thermal energy transfers

Figure 4 :

 4 Figure 4: Vertical view of the model grid 5.1.2 Pumping and injection cycles A first dimensioning of the energy needs for a theoretical greenhouse subjected to the climatic conditions of Nîmes in the South of France has led to the definition of exploitation conditions for the geothermal doublet:-

Figure 5 :Figure 6 :

 56 Figure 5: Comparison of the simulated temperatures in the cold and hot wells for i = 0 (V = 0), and i = 0.4‰ (V = 0.17 m/d) ; Q = 50 m 3 /h, Aquifer thickness = 25m, Cover thickness = 20m, K = 5.10 -3 m/s, ε = 15%.

Figure 7 :

 7 Figure 7: Comparison of simulated "useful" thermal power in the cold and hot wells for i = 0 (V = 0), and i = 0.4‰ (V = 0.17 m/d) ; other parameters same as in Figure 5.

Figure 8 :

 8 Figure 8: Comparison of the simulated temperatures in the cold and hot wells for cover thickness 5 and 20m ; Q = 50 m 3 /h, L = 200m, Aquifer thickness = 25m, Cover thickness = 20m, K = 5.10 -3 m/s, ε = 15%, i = 0.4‰.

Figure 11 :

 11 Figure 10: Comparison of the simulated temperatures in the cold and hot wells for L = 150m and 200m, i = 0.4‰; Q = 50 m 3 /h, Aquifer thickness = 10m, Cover thickness = 20m, K = 5.10 -3 m/s, ε = 15%.

Table 1 : Energy balance for climate conditions of South of France.

 1

			Winter	Spring Summer Autumn Year
	Energy		147	87	3	63	300
	consumption	for
	heating (kWh/m²)
	Energy input Solar	187	498	659	498	1607
	radiation (kWh/m²)
	Maximum	flow	118	242	260	80	21 (*)
	rate (m 3 /h/ha)
	1 ha = 10000 m²

*: the average number of m 3 per m² per hour during a year

Table 2

 2 presents the range of parameters that have been tested by modeling.

Table 2 : Range of parameters tested by modeling

 2

	Parameter	Unit	Range
	Cover thickness	m	5-30
	Aquifer thickness	m	10-30
	Substratum thickness	m	30
	Hydraulic conductivity	m/s	5.10 -4 -2.5.10 -2
	Porosity	%	5-15
	Hydraulic gradient	‰	0-2
	Regional flow	m/day 0-8.7
	Distance between the wells m	150-200
	Pumping/injection discharge m 3 /h	25-100
	Aquifer thermal conductivity W/m/°C 2.09-2.26
	Aquifer calorific capacity	J/m 3 /°C 1.84-2.09.10 6
	Longitudinal dispersivity	m	5-7.5
	Transverse dispersivity	m	1.7-2.5
	5.1.4 Methods to compare simulations

Table 3 : Summary of simulations results for i = 0 and i = 0.4‰ ; other parameters same as in Figure 5.

 3 Influence of the cover thickness Cover thickness will play a role of thermal insulator and reduce the exchange by conduction towards the surface.

	Parameters	Hot well		Cold well
	i (‰)	V (m/d) ΔT (°C)	r (%)	ΔT (°C)	r (%)
	0	0	+6.3	+67	+0.8	+53
	0.4	0.17	+1.7	+26	+0.4	+22

Table 4 : Summary of simulations for cover thickness 5 and 20m ; other parameters same as in Figure 8.

 4

	Parameters	Hot well		Cold well
	Cover thickness	ΔT (°C)	r (%)	ΔT (°C)	r (%)
	5 m	+1.6	+27	+0.0	+19
	20 m	+1.7	+26	+0.4	+23

Table 5 : Summary of simulations for L = 150 and 200 m ; other parameters same as in Figure 10.

 5 : by choice, a negative value traduces the fact that the temperature of the pumped waters in the cold well is greater than the initial temperature of the aquifer.

	Parameters	Hot well		Cold well
	Distance L (m)	ΔT (°C)	r (%)	ΔT (°C)	r (%)
	150	+2.7	40	-0.6(*)	+22
	200	+3.3	42	+0.6	+35

*

  ACKNOWLEDGMENTThis work was financially supported by Brgm, Ctifl, ADEME (French Agency for Environment and Energy Control), and VINIFLHOR (National Office for Fruit, Vegetables, Wine and Horticulture).

		Symbol Unit
	Hydraulic conductivity	K	m/s
	Porosity	ε	%
	Hydraulic gradient	i	‰
	Regional flow (Darcy velocity) V	m/j
	Longitudinal dispersivity Transverse dispersivity	α L α T	m m
	Groundwater temperature	T(t)	°C
	Initial aquifer temperature	T 0	°C
	Fluid calorific capacity Aquifer calorific capacity Aquifer thermal conductivity	γ f γ a λ a	J/m 3 /°C J/m 3 /°C W/m/°C
	Distance between the wells	L	m
	Pumping/injection discharge	Q	m 3 /h
	"Global" thermal power	P g	W
	"Useful" thermal power	P u	W
	"Useful" quantity of energy	E u	J
	Recovery factor	r	%