Georg Daniel Förster

Magali Benoit

Julien Lam
email: julien.lam@ulb.ac.be

Alloy, Janus and core-shell nanoparticles: Numerical modeling of their nucleation and growth in physical synthesis

While alloy, core-shell and Janus binary nanoclusters are found in more and more technological applications, their formation mechanisms are still poorly understood, especially during synthesis methods involving physical approaches. In this work, we employ a very simple model of such complex systems using Lennard-Jones interactions and inert gas quenching. After demonstrating the ability of the model to well reproduce the formation of alloy, core-shell or Janus nanoparticles, we studied their temporal evolution from the gas via droplets to nanocrystalline particles. In particular, we showed that the growth mechanisms exhibit qualitative differences between these three chemical orderings. Then, we determined how the quenching rate can be used to finely tune structural characteristics of the final nanoparticles, including size, shape and crystallinity.

I. INTRODUCTION

Envisioned applications of nanoparticles (NPs), including catalysis [1][2][3][4][5], optics [6][7][8][9][10][11], and medicine [12][13][14][15][16], depend on the ability to tailor their physical properties. A precise control of composition, shape, size distribution, and crystal structure should therefore be achieved during the synthesis [17][18][19]. Most recent research efforts went into raising the complexity of the material and multicomponent systems are increasingly researched in order to combine several properties within the same NP [20][21][22]. For binary systems, three different types of morphology are observed (1) alloy [23][24][25][26], (2) core-shell [27][28][29][30][31][32] and (3) Janus [11,[33][34][35].

Owing to their simplicity and their ability to generate a wealth of structures, physical synthesis methods represent an appealing alternative to more traditional soft chemistry routes and should be crucial to reach an ondemand synthesis of nanomaterials in general [27,32,[36][37][38][39][40][41][42][43][44]. Physical synthesis methods include laser ablation [45][46][47][48], flame pyrolysis [49][50][51], or magnetron sputtering [45,52,53] and always start with the production of an initially hot gas. The subsequent cooling renders the system unstable, which leads to the nucleation and growth of the NPs. These methods are highly versatile and easily controllable via experimental parameters. The out-of-equilibrium state during the growth process and the formation kinetics both play a crucial role in determining particle properties [54,55] and should therefore help obtaining new exotic structures and morphologies [27,32,38,39,43].

However, due to their complexity, details of the NP formation processes are still not very well understood. Experimentally, physical synthesis has been investigated mostly by optical techniques [51]. In particular, at the early stages, when gas and small clusters make up most of the system, the temperature can be probed by infrared [56,57] and Raman spectroscopy [58,59], laser-induced fluorescence [60], and laser-induced breakdown spectroscopy [59,61,62]. For instance, with the latter technique, the chemical composition of diatomic molecules and the electron density were recently measured for laser ablation in air and liquid [61,[63][START_REF] Lam | Pulsed Laser Ablation in Liquid : towards the comprehension of the growth processes[END_REF][START_REF] Amans | [END_REF]. The synthesized nanomaterials can also be characterized in terms of size and crystallinity using scattering methods based on optical or X-ray sources [66][67][68][69]. Understanding the mechanisms occurring in physical synthesis techniques has also benefited from numerical simulations [70]. For instance, models for laser-matter interaction and the subsequent plasma formation have been developed for the study of laser ablation [71].

Nevertheless, the growth of nanostructures in the cooling ablation plume appears to be much harder to study. On the one hand, nucleation and growth in such nonequilibrium processes are very fast and involve few atoms, so that it cannot be easily investigated using conventional experimental techniques. On the other hand, much numerical effort went into the investigation of this problem and gas condensation in general. The temperature quench is usually achieved by successive collisions between the reactive and inert gas atoms kept at room temperature. With such modeling, it was demonstrated that the condensation rate is almost independent of the interatomic potential and that it can be related to the partial pressure of the inert gas using an analytical model [72].

Various monoatomic materials were studied, including silicon [73], germanium [74,75], and metals [72] using semi-empirical potentials such as EAM, Tersoff, and Stillinger-Webber. Yet, for binary systems like alloy, Janus and core-shell nanomaterials, the atomistic mechanisms occurring during the synthesis are not very clear yet, especially in the context of physical methods [76][77][78].

In this work, we study three binary Lennard-Jones systems each tailored to favor one of these morphologies. On the one hand, we identify the differences in the synthesis mechanisms for each morphology by studying the temporal evolution. On the other hand, the final distributions of clusters are analyzed in terms of structural properties, including size, asphericity, crystallinity and, chemical ordering. Moreover, we show that these structural properties can be finely adjusted solely by changing the quenching rate, which can be controlled in most experimental set-ups. Using this simple model system, we aim at drawing very general conclusions regarding the underlying mechanisms which involve both thermodynamical and kinetic processes.

II. METHODS

Three different types of atoms make up the systems studied here. Truncated and shifted Lennard-Jones interactions are used throughout this work. The two first species are reactive, and their interactions are tailored to obtain (1) alloy, (2) Janus and (3) core-shell NPs [See Table I]. In the two latter cases, the coefficients are chosen by following the work of Mravlak et al. who computed the most stable structures of binary nanoclusters for a large number of combination of coefficients [79]. In the alloy case, the σ parameters are chosen to obtain the Kobb-Andersen system, which crystallizes in alloy BCC structures [80][81][82][83], while the energy parameters are set in a way to roughly match the melting and boiling temperatures obtained with the core-shell and Janus systems. The third species serves as an inert buffer gas and interacts with the reactive NPs and with one another also using truncated and shifted Lennard-Jones interactions, but with a cutoff equal to 2 (1/6) σ which leads to a purely repulsive interaction.

Initially, the reactive system consists of 10 5 atoms at a density of 0.02σ -3 . The same number of inert atoms is added in the simulation box. In our molecular dynamics simulations, the reactive atoms move freely (nonthermostated), a Berendsen thermostat controls the temperature of the inert gas. The quenching is obtained by linearly decreasing the inert gas temperature and thus inducing a subsequent cooling of the reactive gas. The temperature is decreased from k B T = 0.7 to k B T = 0.2 during quench durations ∆t, ranging from 10 6 t 0 to 10 7 t 0 . The simulation time step t 0 correpsonds to 0.005τ in Lennard-Jones units. To convert Lennard-Jones reduced units into physical units, one can use Lennard-Jones parameters obtained for aluminium where τ = 0.19 ps and /k B = 5886 K [84]. In this case, the corresponding quenching rates range from 3080 K/ns to 308 K/ns which is faster than experimental rates measured in the case of plasma quenching with laser ablation [START_REF] Amans | [END_REF]. This simulation protocol allows for a good control of the quenching rate without having to finely tune neither the inert gas pressure nor its interactions with the reactive atoms. Fig. 1, shows that initially, the reactive gas temperature follows closely the imposed constant quenching rate. Then, because of condensation and crystallization, potential energy is released that adds to the kinetic energy. This trend was also observed when instead of a linear decrease, the inert gas temperature was fixed at low temperatures throughout the entire simulations [72]. Final results are averaged over ten independent runs. All the simulations were performed using the LAMMPS software package [85]. For the analysis, atoms less than 1.5σ apart from each other are regrouped into clusters. In what follows, any clusters made of more than one atom is referred to as "NP". For each of the obtained NPs, several characterizations are carried out. The size is measured simply by counting the number of atoms and by means of the gyration radius of the NPs, i.e., the root mean square distances of the atoms from the center of mass of the NPs. Furthermore, we quantify the shape of the NPs with the help of an asphericity index,

χ = 3 (2I 1 -I 2 -I 3 ) /2 (I 1 + I 2 + I 3 ), with I 1 ≥ I 2 ≥ I 3
the principal moments of inertia of the NP. The asphericity index varies by construction between 0 and 1, with values close to zero corresponding to rather symmetric geometries such as spheres and larger values corresponding to rod-like and disk-like shapes. In order to determine whether or not an atom is part of a crystallite, we use the adaptive Common Neighbor Analysis method (a-CNA), as implemented in the ovito software package [86][87][88].

Alloy Janus Core-shell ij σij rcut(σ) ij σij rcut(σ) ij σij rcut(σ) A-A 0.

III. RESULTS

A. Chemical ordering and size distribution

In the following two sections, we focus on results obtained at an intermediate quench duration of ∆t = 5×10 6 t 0 , i.e., 4.8 ns in physical units. The final snapshots of Fig. 2 confirm that alloy, Janus, and core-shell were indeed obtained. The NPs can be assigned one of the three classes by defining two order parameters: (i) p 1 is the distance between the center of masses G of atoms of species A and B, and (ii) p 2 the ratio of the fraction of atoms of species A denoted x in the core (atoms within the gyration radius) and the shell (all other atoms). Fig. 4 shows that the core-shell NPs of this work consist of a core of atoms A and a shell of an alloy of A and B. Therefore, in our case, p 2 would equal two for perfect core-shell NPs, whereas it would equal one for ideal alloy or Janus NPs. The definitions of the order parameters are illustrated in Fig. 3(a). The NP is assigned alloy character if both order parameters are less than thresholds set to 1.1 σ and 1.1 respectively, if not the NP is of Janus (p 1 > p 2 ) or core-shell (p 1 ≤ p 2 ) type. For this analysis, we consider only NPs with more than 100 atoms. The majority of the NPs have the desired chemical ordering for all of the three Lennard-Jones models. This is shown in Fig. 3 for ∆t = 5 × 10 6 t 0 , but holds true also for all the other quenching rates (data not shown). However, we note that in some cases, a significant number of NPs solely made up of atoms B, form. As we will see in the following, these NPs are smaller and are formed at the end of the quenching process. Corresponding NPs containing only atoms of species A do not occur in the simulations reported here.

Fig. 4 shows the species-resolved atomic number densities within the largest NP obtained at the end of the simulations as a function of the distance from their respective center of masses. Alloy and Janus NPs ideally show uniform distributions, while in the case of the core-shell NPs, higher concentrations of species A can be found in the center as compared to the shell. However, neither, core nor shell are made up purely of either species A or B. In particular, the shell is composed of almost as many atoms of species B as of species A.

Regarding the size of the obtained NPs, probability distributions of gyration radii are given in Fig. 5 for the three system types. The narrowest distribution is achieved within the alloy system, where 59% of the NPs have a gyration radius in the range of 4.2-6.2 σ. The situation is qualitatively different for the Janus and core-shell systems: the distribution of the gyration radii is bimodal with relatively broad peaks below and above 5σ. The small size peak consists in NPs made mostly of atoms of species B that bind more weakly both with one another Species B Species A FIG. 4. Atomic number densities of atoms of species A (solid lines) and species B (dashed lines) as a function of distance from the center of mass of the NPs. Red, blue and green curves indicate alloy, Janus, and core-shell systems, respectively. The data correspond to the largest NP obtained at the end of the quenching simulations with ∆t = 5.0 × 10 6 t0. and atoms of species A.

B. Nucleation and growth mechanisms

In this section, we investigate the temporal evolution i.e., the emergence of the NPs. In Fig. 2, one can roughly observe three stages (1) binary gas, (2) coexistence of lower density gas and droplets, and (3) coexistence of very low-density gas and crystalline NPs with the desired chemical ordering. The evolution of the fraction of atoms of species A and B in the forming NPs are shown in Fig. 6. While in the alloy case species A and B behave similarly, in the Janus and core-shell cases, twostep mechanisms are observed in which the number of atoms A increases before that of atoms B. This shows that the nucleation and growth mechanisms are qualitatively very different for alloy, Janus and core-shell NPs. Inversely, the gas is dominated by monomers of species B during much of the simulation, and especially in the time frame of 2.0 -4.0×10 6 t 0 in the case of the Janus and core-shell systems. Droplets condensed from this remaining gas make up the small size peak observed in the size distribution [See Fig. 5].

In order to get a more detailed picture, two parameters are evaluated in time, the size of the largest NP and the fraction of atoms that are found within the crystalline parts of NPs [See Fig. 7]. Firstly, the size of the largest NPs increases slowly for t < 2 × 10 6 t 0 , which corresponds to nucleation and growth through monomer accretion, and then they grow in a step-wise manner through collisions with other NPs. Secondly, the crystalline structures emerge later at t > 3 × 10 6 t 0 , which corresponds to k B T < 0.4 , and then plateau rapidly. The figure shows data corresponding to the alloy, Janus, and core-shell system at a quench duration of 5×10 6 t0.

L s

Not only the size of the largest NPs increases during the quench but also, and more importantly, the entire size distribution evolves. For further analysis, the NPs are grouped by size into five classes: 2-10, 11-100, 101-1,000, 1,001-10,000, and 10,001-100,000 atoms. For the quench duration of 5×10 6 t 0 , Fig. 8 shows the abundance of atoms within each of the size classes as a function of time in the case of the alloy, Janus, and core-shell systems. The fraction of atoms within the smallest NPs of up to 10 atoms decreases monotonically with time in the alloy system, which is consistent with results obtained with pure germanium NPs [74,75]. In the case of the Janus and core-shell systems, the curves present local maxima at about the time when significant amounts of atoms of species B condense [See also Fig. 6]. The amount of atoms within the next larger class of NPs (11-100 atoms) exhibits a maximum at a simulation time of ≈ 2.0×10 6 t 0 . While, in the case of the alloy system, this is the only maximum, a secondary lower maximum, which corresponds to the formation of droplets mostly made up of atoms of species B, is found at ≈ 4.0×10 6 t 0 for the Janus NPs. This secondary maximum is not observed for the core-shell NPs either, but the fraction of atoms in NPs of this type and size class increases again at the end of the simulation. Independent of the system type, during a short time window (2.0-2.5×10 6 t 0 ), most of the atoms are bound in NPs of 101-1,000 atoms. This particular intermediate size corresponds to the critical nucleus for the gas to liquid transition as observed at approximately 300 atoms with monodispersed Lennard-Jones particles using umbrella sampling [89]. Later during the simulations, this fraction decreases again at the expense of larger NPs. The amount of atoms within the two largest classes of NPs increases monotonically with time in each of the systems, albeit none of the largest NPs are formed in the case of the alloy system. FIG. 8. Fraction of atoms within NPs of different sizes as a function of time in case of the simulation with a quench duration of 5×10 6 t0 for the alloy, Janus, and core-shell systems.

C. Tuning structural parameters via the quenching rate

In what follows, we explore how the quenching rate influences the probability distributions of the properties of the NPs. We focus on NP size, shape, crystallinity, and composition at the end of the quenching simulations. The fraction of atoms within NPs increases with quench duration [See Fig. 9]. The same tendency was also noted in simulations where the inert gas pressure was varied [72,74,75]. The differences between the three system types -alloy, Janus, and core-shell -are related to their different condensation temperatures as well as kinetic effects. While in the alloy case, atoms of each type may condense at the NP surface, condensation is limited to one side of the NP for either species in the case of the Janus NPs and one species in the case of the core-shell NPs. In either case, however, most (> 90%) of the atoms of the system atoms are bound in NPs. As can be appreciated from Fig. 9(b), not only the fraction of atoms in NPs, also the average size of the NPs increases naturally with the quench duration. This increase in size is particularly pronounced with the alloy and Janus systems. As stated earlier, the size distribution of the Janus and core-shell NPs is bimodal, which translates here in these cases to slightly lower average gyration radii.

FIG. 9. Fraction of atoms bound in NPs as well as the average gyration radii of NPs populations as a function of quench duration for the alloy, Janus and core-shell system at the end of the quenching simulation.

The abundance of NPs within the different size classes as a function of quench duration for the three system types considered in this work is shown in Fig. 10. The fraction of atoms found within the three classes of smallest NP size (up to 1,000 atoms) decreases exponentially (mostly straight lines on the logarithmic plot) with quench duration. While the rate of decrease becomes smaller with NP size in the case of the alloy and Janus systems, it becomes more important for the larger coreshell NPs. Considering the fact that, per NP, the second and third smallest class of NPs contain respectively one and two order of magnitude more atoms than the smallest class of NPs, less of these larger NPs are generated in most cases. The majority of atoms are part of NPs of medium size (1,001-10,000 atoms) at the end of almost all the simulations. Considering the logarithmic scale of Fig. 10, the number of atoms within these NPs reaches a clear maximum at a quench duration of about 4 × 10 6 t 0 , where they essentially make up the entire system. At larger quench durations, the number of atoms in the largest NPs of up to 1 × 10 5 atoms increases exponentially and becomes comparable to the atoms within the second largest class of NPs at the largest quench durations (0.9 -1.0 × 10 7 t 0 ) studied here. Since the simulations used here contain 1 × 10 5 atoms, we cannot observe larger NPs than that. We speculate, however, that their number would be very limited in the range of quench durations discussed here. Their abundance would become more important at larger quench durations at the expense of the two largest classes of NPs shown in Fig. 10. FIG. 10. Fraction of atoms within NPs of different sizes as a function of the quench duration for the alloy, Janus, and core-shell systems.

In the following, we characterize the NPs in terms of their shape, crystallinity, and composition. It turns out that the asphericity of NPs does not depend very much on the quench duration (data not shown), or on the type of system (alloy, Janus, or core-shell). Fig. 11 shows the average asphericity over all quenching rates as a function of NP size. While the smallest and largest NPs have higher asphericity indexes on average, the most spherical NPs have an intermediate size of 101-1,000 atoms. This size range was already identified as a critical size of the gas-liquid transition when studying the temporal evolution. As such, it suggests that the nucleation core should be spherical in these systems. It is, therefore, instructive to consider the internal structure of the NPs: Fig. 11 shows the fraction of crystalline atoms within the NPs, averaged over all quenching rates, as well. Surprisingly, the crystallinity is also little dependent on the quench duration (not shown). However, NP size has a significant effect again: The largest NPs have a higher amount of crystalline atoms. This is to be expected, as larger NPs have relatively fewer surface atoms, which are identified as non-crystalline by the a-CNA method.

The chemical composition of the NPs at the end of the quenching simulation is given in Fig. 12. In almost all our simulations (except for the smallest alloy NPs at quench durations in the range of 2-3×10 6 t 0 ), the amount of atoms of species A increases with NP size and decreases with quench duration. The alloy ratio of the two largest NP populations reaches about unity in all three system types at the highest quench durations studied here. In the alloy system, all NP populations have an alloy ratio between 0.65 and 1.15, meaning they always contain a non-negligible number of atoms of either species. In the case of the Janus and core-shell systems, this ratio tends to zero for NPs smaller than 1,000 atoms. In other words, these NPs are composed mostly of atoms of species B. Atoms of species B condense at lower temperatures and thus later in the quenching simulations (compare Fig. 6). At that time, atoms of species A are already absorbed into other NPs that are then larger at the end of the simulation. Therefore, close to the end of the simulations, there is mostly gas of atoms B left that finally condenses into small, almost monoatomic NPs of species B. Inversely, atoms of species A condense at higher temperatures, and, in the case of the core-shell NPs form the core of the NPs.

IV. CONCLUSIONS

In summary, molecular dynamics simulations were employed in order to study the formation mechanisms oc-curring in the physical synthesis of binary NPs. First, we confirm that with simple Lennard-Jones interactions, one can obtain NPs with different chemical ordering: alloy, Janus, and core-shell. Analyzing the probability distribution, we show that while alloy NPs display narrow size distributions, Janus and core-shell NPs exhibit bimodal distributions. Our simulations reveal that the nucleation and growth mechanisms depend qualitatively on the system type, with two-step processes occurring in Janus and core-shell systems. A critical cluster size emerges at intermediate time which is reminiscent of a classical nucleation core in equilibrium conditions. Studying the temporal evolution, we go beyond equilibrium simulations. In our case, the phase transitions occur at low free energy barriers and kinetic considerations related to quench durations are more likely to explain the formation mechanisms. The quench duration does not affect the crystallinity nor the asphericity of the NPs but remains an essential parameter to tune NP size and the chemical ordering, especially for Janus and core-shell systems. While our modeling follows the simplest possible approach, it allowed for studying numerous features of the physical synthesis of binary NPs. As such, we hope it will encourage further studies with the same modeling, but even more complex structural characterization using for instance novel machine-learning methodologies [90].

FIG. 1 .

 1 FIG. 1. Temperature evolution during the quench in the (a) alloy, (b) Janus and (c) core-shell cases. Continuous and dotted lines represent the temperature of the reactive and inert gas, respectively.

FIG. 2 .

 2 FIG.2. Temporal evolution during the temperature quench with ∆t = 5×10 6 t0. Time and temperature are given in units of 10 6 t0 and of , respectively. In the snapshots, blue and red dots represent A and B atoms. The inert gas atoms are not shown for better clarity.

FIG. 3 .

 3 FIG. 3. (a) Illustration of the order parameters p1 and p2. (b) Chemical ordering for the final NP population obtained at the quench duration of 5 × 10 6 t0.

FIG. 5 .

 5 FIG.5. Probability distribution of NPs from the final configuration of simulations with a quench duration of 5×10 6 t0. Panel (a-c) are concerned with the alloy, Janus, and core-shell systems.

FIG. 6 .

 6 FIG.6. Species-resolved fraction of atoms within NPs obtained at a quench duration of 5 × 10 6 t0 in the (a) alloy, (b) Janus and (c) core-shell cases.

FIG. 7 .

 7 FIG. 7. Number of atoms within the largest NP and overall fraction of crystalline atoms within NPs as a function of time.The figure shows data corresponding to the alloy, Janus, and core-shell system at a quench duration of 5×10 6 t0.

FIG. 11 .

 11 FIG. 11. Asphericity and crystallinity of NPs of different size classes at the end of the quenching simulations.

FIG. 12 .

 12 FIG.12. Ratio of atoms of species A and B within the final NPs grouped into 5 size classes as a function of quench duration.

  support of the Fonds de la Recherche Scientifique -FNRS. Computational resources have been provided by the Consortium des Equipements de Calcul Intensif (CECI) and by the Fédération Lyonnaise de Modélisation et Sciences Numriques (FLMSN).

TABLE I .

 I Lennard-Jones parameters for the interaction between reactive NPs in the three different cases.

	80 1.00 2.5	1.00 1.00 2.5	1.00 1.00 2.5
	B-B 0.40 0.88 2.2	0.62 1.00 2.5	0.50 1.00 2.5
	A-B 1.20 0.80 2.0	0.30 1.00 2.5	0.70 1.00 2.5


Djamel Selloum
email: selloumdjamel@gmail.com

Vincent Techer

Abdellah Henni

Sophie Tingry

Marc Cretin

Christophe Innocent

Performance comparison with different methods for Ethanol/O2 biofuel cell based on NAD + cofactor immobilized and activated by two type of carbon nanoparticles

Keywords: bioelectrochemistry, enzymatic biofuel cell, carbon particles, alcohol dehydrogenase, diaphorase

  Performance comparison with different methods for ethanol/O2 biofuel cell based on NAD+ cofactor immobilized and activated by two types of carbon nanoparticles

Introduction

The significantly environmental impact in the use of unrecoverable energy source will shortly force us to found a substitute way of energy production [START_REF] Surendra | Current status of renewable energy in Nepal: Opportunities and challenges[END_REF]. Fuel cells usually generate electricity through electrochemical conversion of fuels directly into electricity [START_REF] Escalona-Villalpando | Performance comparison of different configurations of Glucose/O2 microfluidic biofuel cell stack[END_REF]. Biofuel cells use biocatalysts as enzymes, in regards to biological cells, to convert chemical energy into electrical energy with a low power source at room temperature and under physiological conditions. The advantage of enzymatic biofuel cells is to exploit the selectivity of enzymes and total biodegradability of the catalysts.

Biofuel cells are an attractive alternative to conventional fuel cells, because they use biological catalysts which eliminate the dependence on precious metal catalysts. Biofuel cells reported in the literature since the 1960's [START_REF] Young | The theoretical aspects of biochemical fuel cells[END_REF][START_REF] Yahiro | Biochimica et Biophysica Acta (BBA) -Specialized Section on Biophysical Subjects[END_REF].

Enzymatic biofuel cells employ enzymes to catalyze chemical reactions, thereby replacing traditional electrocatalysts [START_REF] Kayali | Transport properties of alternative fuel microemulsions based on sugar surfactant[END_REF]. These systems generate electricity under mild conditions through the oxidation of renewable energy sources [START_REF] Patrick | Handbook of fuel cells. Fundamentals technology and applications[END_REF]. The advantages of biocatalysts are reactant selectivity, activity in physiological conditions, and manufacturability. However, biofuel cells show relatively low power densities and limited stability associated with the lifetime of the enzymes. The most studied biofuel cell operates with glucose as fuel and oxygen as oxidant [START_REF] Bullen | Biofuel cells and their development[END_REF][START_REF] Kendall | Fuel cell technology: A sweeter fuel[END_REF][START_REF] Service | Shrinking fuel cells promise power in your pocket[END_REF][START_REF] Selloum | A highly efficient gold/electrospun PAN fiber material for improved laccase biocathodes for biofuel cell applications[END_REF][START_REF] Selloum | Optimized electrode arrangement and activation of bioelectrodes activity by carbon nanoparticles for efficient ethanol microfluidic biofuel cells[END_REF][START_REF] Ma | Effect of immobilization supports on the amperometric response of silver nanoparticles enhanced glucose oxidase electrodes[END_REF][START_REF] Niiyama | High-performance enzymatic biofuel cell based on flexible carbon cloth modified with MgO-templated porous carbon[END_REF], in the presence of redox mediators to achieve efficient catalyst utilization.

Enzymatic biofuel cell systems offer main advantage to traditional fuel cells. Enzyme catalysts allow for a wide a variety of fuels to be utilized, because enzymes typically do not react with or get passivated by impurities in the fuel [START_REF] Mccormick | Photosynthetic biofilms in pure culture harness solar energy in a mediatorless bio-photovoltaic cell (BPV) system[END_REF].

Researchers have been attempting to develop direct ethanol fuel cells (DEFCs), but there have been problems because traditional precious metal catalysts (Pt-based catalysts) are unable to efficiently catalyze the oxidation ethanol and maintain an electrode with minimal fouling at low temperatures. However, living organisms are capable of efficiently catalyzing the oxidation of ethanol. Living organisms, such as Pseudomonas aeruginosa [START_REF] Veer Raghavulu | Comparative bioelectrochemical analysis of Pseudomonas aeruginosa and Escherichia coli with anaerobic consortia as anodic biocatalyst for biofuel cell application[END_REF], acetobacter [START_REF] Karthikeyan | Bioelectrocatalysis of acetobacter aceti and Gluconobacter roseus for current generation[END_REF], and gluconobacter [START_REF] Reshetilov | Testing of bacteria gluconobacter oxydans and electron transport mediators composition for application in biofuel cell[END_REF], contain enzymes that can oxidize a variety of alcohols, including ethanol.

An oxidized nicotinamide adenine dinucleotide (NAD + )-dependent dehydrogenase, such as alcohol dehydrogenase, has commonly been utilized in an enzymatic anode in biofuel cells [START_REF] Moore | Improving the environment for immobilized dehydrogenase enzymes by modifying nafion with tetraalkylammonium bromides[END_REF]. Such enzymes oxidize the corresponding fuels with the accompanying reduction of cofactor NAD + to nicotinamide adenine dinucleotide (NADH). In order to collect electricity and recover NAD + for a continuous reaction, the enzymatically produced NADH needs to be oxidized at the supporting electrode. Since the oxidation of NADH requires a large overvoltage on a conventional electrode, typically diaphorase (an enzyme for NADH oxidation) [START_REF] Antiochia | Kinetic and thermodynamic aspects of NAD-related enzyme-linked mediated bioelectrocatalysis[END_REF][START_REF] Ogino | Reactions between diaphorase and quinone compounds in bioelectrocatalytic redox reactions of NADH and NAD+[END_REF]. This paper describes the development of an ethanol biofuel cell based on NAD + cofactor immobilized and activated by two types of carbon nanoparticles. In particular, the work focused on the influence of carbon nanoparticles type on the electron transfer in attempt to improve efficiency of the bio-electrocatalytic process and in the first time, the performances of the devices were evaluated using two different methods.

Experimental Section

Materials

Laccase from Trametes Versicolor (20 U mg -1 solid), Diaphorase 3-20 U mg -1 , Alcohol Dehydrogenase 300 U mg -1 , β-nicotinamide adenine dinucleotide sodium salt (NAD + ), 2methyl-1,4-naphthoquinone (VK3), Aceton, polyethylenimine (PEI), 2,2'-azinobis(3ethylbenzothiazoline-6-sulfonate) diammonium salt (ABTS), Nafion  solution (5wt%), sodium phosphate dibasic dihydrate (Na2HPO4.2H2O) and sodium phosphate monobasic monohydrate (NaH2PO4.H2O) were purchased from Sigma-Aldrich and used without further purification. The phosphate buffer was prepared with Na2HPO4.2H2O and NaH2PO4.H2O pH 5, 7 or pH 9, 0.1 M. The carbon nanoparticles powder as Super P  and KS6 were purchased from TIMCAL.

Bioelectrodes preparation

The biocathode to be employed in the electroreduction of oxygen was prepared by adsorption of enzymes and mediators on the surface of the electrodes by drop casting. 333 µL of laccase 15 mg mL -1 and Super P  15 mg mL -1 in phosphate buffer 0.1 M (pH 5) solution was mixed on a vortex mixer. Sequentially, 100 µL of the solution was mixed with ABTS 5.4 mg mL -1 and 10 µL Nafion  . Then, 6 µL of the preparation was coated onto Au electrode and left to dry at room temperature before keeping in a low humidity environment.

The bioanode to be employed in the oxidation of ethanol was prepared by adsorption of successive coatings separated by a dried step at room temperature. 167 µL of ADH 30 mg mL - 1 and KS6 15 mg mL -1 in phosphate buffer 0.1 M (pH 7) solution was mixed on a vortex mixer and 6 µL of the preparation was pipetted onto the electrode and dried at room temperature. The same procedure was conducted for the immobilization of NAD + 30 mg mL -1 and then diaphorase 20 mg mL -1 . The last coating on the electrode consisted in pipetting 10 µL of VK3 60 mg mL -1 , 190 µL acetone and 10 µL PEI, followed by drying.

Physical characterization

Carbon powders were analyzed by scanning electron microscopy (SEM, S-4800, Hitachi), and X-Ray diffraction (PANAlyticalXpertPRO diffractometer equipped with a X'celerator detector using Ni-filtered Cu-radiation). Raman spectra were acquired using a Raman microspectrometer (Jobin-Yvon S.A.S., Horiba, France) equipped with an optical microscope under ambient conditions with a 647.1 nm laser and the acquisition range was 1000-2000 cm -1 .

Electrochemical measurements

The bioelectrodes were characterized separately by polarization curves performed in dioxygen-saturated phosphate buffer pH 5 0.1 M for the biocathode, or in phosphate solution pH 9, 0.1 M with 160 µL ethanol for the bioanode, after stabilization of the open circuit potential. Electrochemical measurements were performed on a potentiostat Autolab (Eco chemie, Netherlands) at 25°C in phosphate buffer, with a conventional three-electrodes system composed of a stainless-steel auxiliary electrode, a calomel saturated reference electrode and the electrode material as working electrode.

The BFC was characterized by the same electrochemical equipment. The cathode was connected to the working electrode, and the counter and reference electrodes were both connected to the anode. The catholyte solution consisted of phosphate buffer solution 0.1 M (pH 5) saturated with O2. The anolyte solution contained ethanol 160 µL in phosphate buffer 0.1 M (pH 9). Both solutions are separated by ion exchange membrane (Nafion from Dupont USA) [START_REF] Zebda | Hybrid layered double hydroxides-polypyrrole composites for construction of glucose / O2 biofuel cell[END_REF] The open circuit potential (Voc) was measured when no current flows. An equilibration time was maintained during 5min before data collection to reach maximum open circuit potential.

Polarization curves, known as I-V curves, were obtained by measuring the current while slowly varying the potential between the Voc and 0 V, or by applying external resistances from 100 Ω to 1 MΩ and measuring the resulting current and voltage. Current and power densities were determined using the geometric area of the electrodes.

Results and Discussion

Influence of carbon powder.

At first the carbon KS6 and Super P ® were characterized by SEM microscopy (Fig. 1A and1B). KS6 has a sheet structure composed of particle size dimensions of the order of 1 µm and by an intense X-ray diffraction peak at 2θ=26° typical of graphitic carbon [START_REF] Szabó | Composite graphitic nanolayers prepared by selfassembly between finely dispersed graphite oxide and a cationic polymer[END_REF] (Fig. 1C).

At contrary, the carbon Super P ® presents a porous structure composed of smaller nanoparticules of ten nanometers (Fig. 1B) and no characteristic peak on the XRD diffraction pattern (Fig. 1C). Their impact has been studied on the electrochemical behavior of the bioelectrodes based on immobilization of laccase as cathode and diaphorase (DI) and alcohol deshydrogenase (ADH) as anode. The Raman spectrum of KS6 carbon powder presented in Figure 1.D displays a peak at 1582 cm -1 , which correspond to the crystal graphite [START_REF] Wang | Raman Spectroscopy of Carbon Materials: Structural Basis of Observed Spectra[END_REF].This result is in agreement with the data provided by the supplier (TIMCAL).

Characterization of bioanode.

The working of bioanode based on alcohol deshydrogenase and diaphorase immobilization is presented in figure 2. Ethanol is oxidized to aldehyde by ADH and its cofactor NAD + . The regeneration of the cofactor is well achieved using an additional redox protein diaphorase which operates efficiently in the presence of the mediator VK3 [START_REF] Sakai | A high-power glucose/oxygen biofuel cell operating under quiescent conditions[END_REF]. The electrochemical behavior of the bioanode has been studied in the presence of the different carbon black powders. Fig. 3 shows the catalytic activity of the bioanode in phosphate buffer solution pH 9 in the presence of ethanol.

The bioanodes exhibit activity for ethanol oxidation with a Voc around -0.32 V/SCE and -0.30 V/ECS for KS6 and Super P respectively. (Fig. 3). Comparison between carbon powders KS6 and Super P ® . Scan rate 3.33 mV s -1 .

The high value of the peak is due to the large amount of VK3 reduced at the electrode when the substrate reacts with enzymatic chain. A decrease of the signal is observed during the second scan. The resulting values after stabilization gave the same results for electrode with Super P and KS6. Thus, electrodes with KS6 have been chosen for the elaboration of bioanode.

A poor stability has been observed by voltamperometric characterization with a loss of 31% in current density after 3 scans (Table 1). A decrease of 90% of current has been observed after one day of storage at 4°C. These weak performances can be explained by the denaturation of enzyme and the loss of mediator (particularly loss of NAD + ) in the solution. On the other hand, the enzymatic cascade with alcohol deshydrogenase and diaphorase and their mediators generates a local pH modification when the VK3 is not oxidized at the electrode (open circuit).

However, the chronoamperometric characterization presented in figure 4 showed stabilization of current after 1500 s. This result demonstrates the efficiency in operational conditions of the enzymatic catalysis on the electrode. At first a decrease of current is recorded followed by a steady state which corresponds to the stationary operation of the bioanode. Without enzyme no current is observed and also without VK3 (see fig 4). The presence of enzyme and VK3 is necessary to activate the catalytic cycle of enzymatic reaction according to the scheme presented in figure 2.

According to these results we can conclude that bioanode is able to work in the operating condition but a large denaturation is observed in the storage condition. This result can be explained by the reaction association on the electrode. When no potential is applied, or no connection with another electrode (cathode in the case of biofuel cell) substrate residues unremoved by rinsing, can react with enzyme at the electrode and activate the catalytic chain until the formation of reduced vitamin K3. In this condition, a high local concentration of reduced vitamin K3 occurs on the electrode. Moreover, when the catalytic chain stops by lack of substrate, the ratio of reduced and oxidized compound in the film and the variation of local pH can affect the enzymes integrity. Accordingly the storage stability is very poor for this kind of the electrode.

To keep the integrity of the electrode, operating is necessary to evacuate electron produced by the enzymatic reactions.

Characterization of biocathode

The cathode based on immobilization of ABTS and laccase has been studied in the same way. The electrochemical reactions are depicted in figure 5. Reduction of oxygen is catalyzed by enzyme laccase and ABTS as mediator.

Figure 5.Schematic representation of the biocathode system

Fig. 6 shows the electrochemical response of the biocathode towards O2electrocatalytic reduction in phosphate buffer solution pH 5 saturated with O2. From the polarization curves, the oxygen reduction current begins at 0.6 V/ECS, with overpotential of 50 mV, and the current densities feature a semi-plateau that indicates the control of the electrocatalytic reaction by diffusion of the oxygen to the electrode surface. The biocathode based on Super P ® ( ▲ ) shows improved kinetics of the electronic transfer and higher performance increase 500 µA cm -2 in 2folds compared to biocathode based on KS6.

13

The first outcome to be observed from polarization curves is that the addition of carbon particles to the enzyme layers enhances bioelectrocatalytic processes. The large surface area of the particles enhances the reactive species loading and thus the available active sites for electron generation. Besides, the electronic conductivity of the carbon nanoparticles can enhance electron transfer rate to the electrodes. However, it is evident that bioelectrodes based on carbon Super P ® display the best performances. The benefit of Super P ® may be explained by its high porous structure that offers a close proximity to the reactive species and improves diffusion of oxygen within the enzyme films. The storage stability is reported in figure 7. On the opposite of the bioanode, biocathode shows a small decrease of activity and 20 % of initial current is still recorded after 28 days. This result demonstrates the good stability of the biocathode, probably due to the stability of laccase immobilized [START_REF] Brunel | Oxygen transport through laccase biocathodes for a membrane-less glucose/O2 biofuel cell[END_REF][START_REF] Merle | Long-term activity of covalent grafted biocatalysts during intermittent use of a glucose/O2 biofuel cell[END_REF] in our experimental conditions even if 20 % of initial value of current is recorded after 28 days [START_REF] Manoand | O2 Reduction in Enzymatic Biofuel Cells[END_REF]. A decrease of 50 to 60 % of the initial value has been already observed after 12 h [START_REF] Franco | Product Analysis of Operating an Ethanol/O2 Biofuel Cell Shows the Synergy between Enzymes within an Enzymatic Cascade[END_REF]. Stability for a long time is a real limitation of this kind of system.

The operational stability has been studied by chronoamperometric curves characterization at 0.5 V/ECS. Figure 8 shows a current evolution with time. The initial value increases rapidly and reaches a stable value at 5000 s. This characterization confirms the stable response of the electrode in operating conditions.

Biofuel cell characterization

The good operational stability of bioanode based on carbon KS6 and biocathode based on carbon Super P allow application of these electrodes for biofuel cell construction.

The polarization curve and power output are depicted in figure 9. The initial open circuit potential of the biofuel cell was 1.02 V and the maximum power density was 300 µW.cm -2 . These results are in agreement with performances obtained with other ethanol biofuel cell [START_REF] Topcagic | Development of a membraneless ethanol/oxygen biofuel cell[END_REF][START_REF] Deng | A membraneless biofuel cell powered by ethanol and alcoholic beverage[END_REF]. In our case, the open circuit potential value is higher than the value previously reported. The immobilization of mediator, specifically NAD + , on the surface of electrode could explain this large open circuit potential. Indeed, all reacts are immobilized on the electrode and the potential is defined by amount of redox compound on the surface. A large amount of NAD + is located close to the enzyme and could react quickly without diffusion limitation.

The electrochemical characterization has been compared with power output given by connection of biofuel on the electrical resistance. In this last case, no potential and no current are applied on the electrodes. The values of potential and current are determined by the steady state of the both electrodes. An overshoot is observed in figure 9 corresponding to power density plot double back to the lower current after the maximum value is achieved. This kind of overshoot is frequently reported in the literature for the bioelectrochemical system based on the biological entities connected to the electrode [START_REF] Angioni | Novel composite polybenzimidazole-based proton exchange membranes as efficient and sustainable separators for microbial fuel cells[END_REF]. It can be interpreted as the limitation of electrochemical circuit in presence of high current [START_REF] Pasternak | Comprehensive Study on Ceramic Membranes for Low-Cost[END_REF].

Figure 10 presents the evolution of current with time for two resistances: 400 and 40 000 Ω and electrical scheme used to connect the biofuel cell.

Surprisingly, the value of current with the larger resistance was highest. This result was thus not in accordance with the omic law which impose that when the resistance increases, the current decreases at constant potential.

This result suggests thus, that the value of potential can be affected by the current delivered by the biofuel cell. This point is clearly showed on figure 9, with polarization curve and power output. The both curves present a decrease of intensity for high values of potential and power.

On the polarization curve, the current decrease for a potential below to 0.55 V. The changing of shape of the curve can be explained by the modification of internal resistance. In other hand, if the resistance is maintained constant, modification of potential can be observed. immobilization. Moreover, good performances have been recorded in terms of power density (300 µW.cm -2 ). However, in the case of multi enzymatic reaction it is necessary to consider with great discernment the polarization measurement because of effect of potential applied on the enzymes. The current production with the high resistance gave good performances; 90 µA.cm -2 during 20 hours.

Figure 1 .

 1 Figure 1.SEM images of carbon powders KS6 (A) and Super P ® (B), scale bar is 1.2 µm.

Figure 2 .

 2 Figure 2.Schematic representation of the bioanode system.

  the presence of Super P  , a sharp increase of the current density and an oxidation peak at -0.15 V/ECS are observed, corresponding to the oxidation of the mediator VK3 on the electrode and also the oxidation of ethanol in the solution (see fig2). In absence of substrate (ethanol) the value of current is very low, demonstrating the effect of ethanol oxidized by enzymatic reaction.

Figure 3 .

 3 Figure 3. Polarization curve (first scan) of an ADH/NAD + /Di/VK3 modified electrode in

Figure 4 .

 4 Figure 4. Chronoamperometric response at -0.2 V/SCE for a bioanode with carbon KS6 in

Figure 6 .

 6 Figure 6.Polarization curve of an ABTS/laccase modified electrode in phosphate buffer solution (pH 5, 0.1 M) under O2. Comparison between carbon powders KS6 and Super P ® . Scan rate 3.33 mV s -1 .

Figure 7 .

 7 Figure 7.Biocathode stability, results based on polarization curves. Electrodes stored at 4°C

Figure 8 .

 8 Figure 8.Chronoamperometric response at 0.5 V/SCE for a biocathode with carbon Super P ®

Figure 9 .

 9 Figure 9. Polarization and power density curves of the biofuel cell using carbon KS6 at the anode and carbon Super P ® at the cathode. Scan rate = 3.33 mV s -1 , with (black curve) and

Figure 10 .

 10 Figure 10.Evolution of the current density while the biofuel cell is connected to different

Table 1 .

 1 Stability of a KS6 modified bioanode

			KS6 anode
		V 0 (mV)	i (A.cm -2 ) -0.2	Current density lost
			V
	1 st scan	-415	680	-
	2 st scan	-420	614	-10%
	3 st scan	-411	470	-31%
	Day+1	-198	06	-90%


L_E

L_H

Keywords: Combined Cycle, cogeneration, security constrained unit commitment, MILP. I

This paper investigates an optimal scheduling method for the operation of combined cycle gas turbines (CCGT). The objective is to minimize the CO2 emissions while supplying both electrical and thermal loads. The paper adopts a detailed model of the units in order to relate the heat and power outputs. The grid constraints as well as system losses are considered for both the electrical and thermal systems. Finally, the optimal power dispatch lies on the hybridization of a Mixed Integer Linear Programing (MILP) scheduling with a greedy search method. Different sets of simulations are run for a small 5-bus test case and a larger model of Jurong Island in Singapore. Several load levels are considered for the heat demand and the impact of the steam pipe capacities is highlighted.

Combined heat and power (CHP) cogeneration systems are widely used assets for both residential and industrial applications due to their high global efficiency as well as economic and environmental benefits. Operation of such units has been extensively addressed in the literature, mainly targeting total revenue maximization and/or emissions minimization. Optimal scheduling problems are commonly investigated to supply both heat and electrical loads in the presence of storage devices [START_REF] Hashemi | A developed offline model for optimal operation of combined heating and cooling and power systems[END_REF]- [START_REF] Kia | Optimal day ahead scheduling of combined heat and power units with electrical and thermal storage considering security constraint of power system[END_REF] or renewable energy sources [START_REF] Shao | Modeling and Integration of Flexible Demand in Heat and Electricity Integrated Energy System[END_REF]. Another class of related studies refers to planning problems with the sizing of CHP generators [START_REF] Shao | Modeling and Integration of Flexible Demand in Heat and Electricity Integrated Energy System[END_REF] or the design of waste heat recovery networks [START_REF] Li | Transmission-Constrained Unit Commitment Considering Combined Electricity and District Heating Networks[END_REF]. Most of those optimization frameworks assume coarse models for the generators. Usually, an operating region for the units power and heat outputs is assumed resulting in linear constraints for the scheduling problems [START_REF] Hashemi | A developed offline model for optimal operation of combined heating and cooling and power systems[END_REF]- [START_REF] Kia | Optimal day ahead scheduling of combined heat and power units with electrical and thermal storage considering security constraint of power system[END_REF]- [START_REF] Li | Transmission-Constrained Unit Commitment Considering Combined Electricity and District Heating Networks[END_REF]. Another shortcoming often encountered in the literature lies in the lack of representation of both thermal and electrical networks. Recent studies have pointed out that coupling of thermal and electricity grids in the cogeneration system would bring additional flexibility, as well as complexity, to the system operation regime, making the optimization of such cogeneration systems difficult [START_REF] Li | Transmission-Constrained Unit Commitment Considering Combined Electricity and District Heating Networks[END_REF]- [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF]. This paper aims at Rémy Rigo-Mariani, Chuan Zhang, Alessandro Romagnoli, Markus Kraft, Keck Voon Ling, Jan Maciejowski Cambridge Centre for Advanced Research and Education in Singapore Ltd, 1 Create Way, 138602, Singapore Nanyang Technological University, 50 Nanyang Ave, 639798 Singapore A Combined Cycle Gas Turbine Model for Heat and Power Dispatch Subject to Grid Constraints addressing the two aforementioned points. First, a realistic model of combined cycle gas turbine (CCGT) is considered. In particular, that representation taken from previous work [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF] estimates the CO2 emissions of the units all along their operating range, by the computation of the chemical reaction within the combustion chamber. Secondly, the optimal scheduling of the units takes both power and heat network constraints into account, as well as the losses within the cables and the pipes. The last contribution of the paper lies in the hybridization of a Mixed Integer Linear Programming (MILP) scheduling with a greedy search method. The objective is to avoid prohibitive computational times when large numbers of units are considered. Also note that with recent developments in cyber-physical systems for energy system management (e.g. J-Park Simulator [START_REF] Zhang | Knowledge management of eco-industrial park for efficient energy utilization through ontology-based approach[END_REF]) , the model proposed in the paper could be integrated with real data from smart meters to unleash the potential of real-time control of CHP systems. The rest of the paper is organized as follows. Section III presents the model for the CCGT as well as the MILP formulation that allows one to compute the power/heat outputs. Section IV integrates that model in an optimal dispatch problem whose objective is the minimization of the CO2 emissions over a representative day (no cost of generation are considered here). Then, Section V introduces the representation of the power grid and heat network. Finally, sets of simulations are run for a small 5-bus test case and a model of the Jurong Island in Singapore.

III. CCGT MODEL

A. CCGT Operation for Heat and Power Generation

The work presented in this paper has been done in the context of the C4T project whose objective is carbon reduction in industrial/chemical activities. Jurong Island in Singapore is considered as a case study. That island hosts half the power generation capacity in Singapore. Most of the production is consumed on site by heavy customers (oil refineries and chemical facilities) while the surplus is exported to the Singaporean mainland. Previous work focused on system modeling with specific attention given to the representation of the combined CCGT on site [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF]. The operation of such units relies on the combination of two thermodynamic cycles. In the top Brayton cycle, the input air-mass flow (ma in kg/s) at ambient conditions is compressed by a ratio rc before entering the combustion chamber and being heated to the turbine inlet temperature (TIT in O C). Power is then produced by the expansion of the hot gases in the gas turbine (GT) whose shaft is connected to the generator. In the bottom Rankine cycle, the hot gases at the turbine-exhaust temperature (TET in O C) are used in a heat recovery steam generator (HRSG) to produce superheated steam at high temperature for a steam turbine (ST) with an efficiency ηst. Thus, the combined-cycle (CC) units allow exploitation of the available heat in the exhaust gas of the gas turbine and enhancement of the power generation efficiency from 40% (ηgt) to 60% (ηcc) [START_REF] Boyce | Handbook for cogeneration and combined cycle power plants[END_REF]. A conventional way to operate these units consists of the control of the inlet guide vane (IGV) of the compressor. The IGV is controlled up to a 30% air-flow reduction with a constant TET followed by a control of the maximum air flow [START_REF] Kim | Comparative analysis on the part load performance of combined cycle plants considering design performance and power control strategy[END_REF] (Fig. 1a). A model was developed in [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF] to estimate the performance and the CO2 emissions of different technologies over their whole operating range (typically from 30 % to 100 % of the nominal electrical output). Fig. 1b displays a sample of the obtained results, with increased efficiencies at the designed operating point (i.e. 100 %). The developed model allow one to identify the share between heat and power in the units all along the two cycles (Fig. 2a) and for different operating points (Fig. 2b). Note that the estimated amount of heat in the GT exhaust gas is significant. This justifies the use of a combined cycle to enhance the efficiency of the energy generation, with an overall electrical output computed as the sum of both GT and ST powers.

The previous work only considered the electrical output and performed an environmental unit commitment (UC) in order to supply a given power load profile over a representative day while minimizing the corresponding CO2 emissions. In this paper the CCGT units are considered as cogeneration assets with the ability to provide both heat and power, with the waste heat transferred into a thermal network. In addition, the amount of transferred heat is adjustable by varying the proportion of GT exhaust gases dedicated to the power generation at the ST stage. A power/heat ratio 𝛼 i,t ph (between 0 and 1) is then introduced considering the operation of unit i at time t. Thus for a given power output of the gas turbine P i,t gt (in MW) different operating conditions can be computed and correspond to distinct values for the overall electrical and heat outputs of the unitrespectively P i,t e and P i,t h in MW. Fig. 3 shows the results obtained while running the model for a 400 MW CCGT unit with a 240 MW GT and a 160 MW ST. As expected, both electrical and heat outputs increase with the operating point of the gas turbine. A case with 𝛼 i,t ph =1 corresponds to the nominal operation described in Fig. 1 with the maximum electrical output reached for a gas turbine operating at 100 % (Fig. 3a). That power generation significantly decreases when the power/heat ratio becomes lower, leading to more waste heat (Fig. 3b). In such a case the priority is given to the thermal generation. Only a small part of the nominal waste heat coming from the GT is available. Consequently, the output of the steam turbine is assumed to be null and the CCGT power only corresponds to the gas turbine output.

B. MILP Formulation for CCGT Operation

As will be presented in the next section, the objective of the proposed heat and power dispatch is to supply both thermal and electrical loads with minimal emissions of CO2. A MILP approach is considered in order to ensure reasonable computational times, as is commonly done in the literature when the commitment of a large number of units with temporal constraints is investigated Erreur ! Source du renvoi introuvable.- [START_REF] Li | Transmission-Constrained Unit Commitment Considering Combined Electricity and District Heating Networks[END_REF]. The gas turbine output and the power/heat ratio of the units are considered as controls. Additional variables and constraints need to be introduced in order to model both electrical and thermal outputs for a given set of controls. In other words, the concern here is to represent the surfaces plotted in Fig. 3 in the MILP formulation. Thus linearization methods for functions of two variables are considered [START_REF] D'ambrosio | Piecewise linear approximation of functions of two variables in MILP models[END_REF]. Traditional triangular or rectangular approaches offer the best performance but require a great number of variables (continuous and binary) that could lead to prohibitive computational times. For the sake of simplicity, a 1D method is considered here. With functions of two variables, the idea is to keep one continuous variable (P i,t gt in this case) while the other is discretized. Thus 𝛼 i,t ph is substituted by a binary variable 𝛼 i,t,r ph_b with r ∈ R the set of possible values for the power/heat ratio. Arbitrarily, a set of eleven different values is chosen and corresponds to the discretization along 𝛼 i,t ph displayed in Fig. 3 (i.e. r = 1 for 𝛼 i,t ph = 0 and r = 11 for 𝛼 i,t ph = 1). Thus 𝛼 i,t,r ph_b =1 allows one to identify at which ratio r a unit i is operated at time t. The electrical and heat outputs are then computed by defining a special ordered set over P i,t gt [START_REF] D'ambrosio | Piecewise linear approximation of functions of two variables in MILP models[END_REF] with a two-block piecewise linearization (set S). The following equations describe the set of variables and constraints to model the electrical power output (subscript e) -the thermal components (subscript h) are treated similarly. First, an additional binary variable w 𝑖,t,r,s 𝑒 is introduced and is equal to 1 if unit i is operated with ratio r in block s at time t. Constraints (1) ensure that only one operating block is identified provided that the unit is on at time t (i.e. ui,t = 1). A continuous variable 𝛽 𝑖,t,r,k 𝑒 represents the weight coefficients attached to the k breakpoints that border each block s (K=S+1). With constraints in (2), only the weights for the two breakpoints around the operating block are non-null. sx-1 and sx+1 denote the two segments around a breakpoint kx with w 𝑖,t,r,s x 𝑒 = 0 for the "extreme" breakpoints (i.e. x=0 and x=K+1). As the sum of those weights is equal to one by (3), the GT power is computed following (4) with the breakpoints G i,t,r,k gt_e entered as parameters. Fig. 4 shows an example for the GT power computation of a unit operating in a specific block s1 for a given heat/power ratio. Constraints ( 5) are introduced over the whole set R in order to compute the electrical output corresponding to the operating power/heat ratio and the calculated 𝑃 𝑖,𝑡 g𝑡_𝑒 . With an appropriate "Big M" value (typically 10 6 ) only the constraints referring to 𝛼 i,t,r ph_b =1 will be active [START_REF] D'ambrosio | Piecewise linear approximation of functions of two variables in MILP models[END_REF]. Equation ( 6) ensures that only a single value of the power/heat ratio is considered for every unit at each time step. 𝑃 𝑖,𝑡 g𝑡_𝑒 is then computed similarly to the GT output with the weight coefficients for the breakpoints G i,t,r,k e entered as parameters.

   S s t i e s r t i u w , , , ,             e s r t i e s r t i e k r t i e k r t i x x x w w 1 1 , , , , , , , , , , , , 0  
     K k t i e k r t i u , , , ,    e gt k r i K k e k r t i gt t i G P _ , , , , , , ,                                 b ph r t i e k r i K k e k r t i e t i b ph r t i e k r i K k e k r t i e t i M G P M G P _ , , , , , , , , _ , , , , , , , , 1 1 
        R r b ph r t i 1 _ , ,   
e k r i G 1 1 , , e gt k r i G _ , , 2 1 e gt k r i G _ , , 3 1 e k r i G 2 1 , , e k r i G 3 1 , , otherwise 0 , , ,  e s r t i w operating point 2 , , , 1 , , , 2 1 1 1
,

      e k r t i e k r t i e gt k r i e gt k r i gt t i G G P _ , , 2 _ , , 1 , 2 1 1 1       e k r i e k r i e t i G G P 2 1 1 1 , , 2 , , 1 ,       β 1 β 2 1 2 1 1 1 , , , , , ,   e k r t i e k r t i   breakpoint block s 1

IV. HEAT AND POWER DISPATCH

A. Classical Unit Commitment MILP Formulation

The model developed in the previous section is now included in an environmental unit commitment for both heat and power dispatch. As already mentioned, the objective is to supply power (P t L_e ) and heat (P t L_h ) load profiles with a set of different CCGT units and with minimal CO2 emissions. A classical twoblock (c ∈ C) linearization of the carbon cost is considered [START_REF] Carrión | A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem[END_REF] regarding the GT output. The objective over the time horizon is computed following [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF] with an hourly time step. That function depends on the operating block for the GT of unit I at time (p i,t,c ) and considering the corresponding block slope Ai,c, the base cost C0i and the start-up cot SUiwith vi,t =1 when unit i starts up at time t. An additional boiler with CO2 coefficient Abl (in kg/MWh [START_REF] Zhang | Life Cycle Emissions and Cost of Producing Electricity from Coal, Natural Gas, and Wood Pellets in Ontario, Canada[END_REF]- [START_REF] Eggleston | IPCC guidelines for national greenhouse gas inventories[END_REF]) is considered to supply the heat load if the units' heat outputs are not enoughsurplus denoted as P t bl (in MW). In order to ensure the convergence a thermal damp load P t dp also needs to be introduced if the units generate too much heat. Thus, implicitly priority is given to the supply of the electrical load in the current implementation. Note that both P t bl and P t dp are unbounded positive variables.

bl t bl T t I i t i i C c gt c t i c i t i 0i p v u P A v SU p A u C gt c t i t i t i                     , , , , , , , min ,

 

Constraints ( 8)- [START_REF] Boyce | Handbook for cogeneration and combined cycle power plants[END_REF] ensure that all the units properly work in the linear block identified in the objective function. Typical operating constraints in UC refer to minimum up and down times for the units, ramping limits or shut down cost [START_REF] Van Den Bergh | A Mixed-Integer Linear Formulation of the Unit Commitment Problem[END_REF]. For clarity they do not appear here and only the logical constraints [START_REF] Kim | Comparative analysis on the part load performance of combined cycle plants considering design performance and power control strategy[END_REF] are considered. Finally, ( 12) and ( 13) allow one to fulfill both power and heat balances at each time step. The problem is formulated in MATLAB using YALMIP [START_REF] Lofberg | YALMIP: a toolbox for modelling and optimization in MATLAB[END_REF] and solved using CPLEX 12.7.1 (16 threads in parallel, 16 GB RAM, 3.2 GHz processor).

B. Results for a Two Units Dispatch

Initially, simulations are performed for a dispatch with two units. Unit 1 is a high efficiency 400 MW CCGT with specific CO2 emissions of 390 kg/MWh (at its nominal point). Unit 2 has higher specific CO2 emissions with 450 kg/MWh (at its nominal point) for a maximum capacity of 250 MW. More information on the considered technologies and the modeling aspects can be found in [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF]. Three different strategies are investigated. In S0 the two units are optimally dispatched to feed the electrical load while the heat demand is supplied by the boiler. The same dispatch strategy is considered in S1 but the CCGT waste heat is injected into a thermal network. The boiler and damp load allow to adjust the amount of heat provided. Finally, S2 denotes the power/heat ratio management previously described. For S0 and S1 the optimal results correspond to the cleaner unit working at its maximum electrical output while Unit 2 provides the surplus of energy during the peak period (Fig. 5a). The ability to control the power/heat ratio tends to lower the power output of Unit 1 in order to favor its heat generation (Fig. 5b). At the same time the operating point of Unit 2 (in terms of GT power) is increased with higher power and heat outputs in S2. Note that with a significant electrical load compared to the heat demand, the power/heat ratio cannot reach values below 0.8 to ensure convergence (Fig. 6). In S2 more waste heat from the CCGT units is transferred to the thermal network which lowers the need for the additional boiler. Consequently, the overall amount of CO2 generated while supplying heat and power is reduced.from 5214 tons to 5195 tons in the previous simulations with specific emissions of 300 kg/MWh for the boiler. Obviously, the CO2 reduction becomes more significant with 'dirtier' boilers as shown in Table I. The same observation can be made when the level of the heat load increases regarding the power demand with 𝛿 = P t L_h /P t L_e (simulation in Fig. 5 corresponds to δ = 75 %). For higher heat load levels the improvement provided by the management of the power/heat ratios is more significant regardless of the boiler emissions (Fig. 7). C. Constant Power/Heat Ratios The computational time for the previous simulation with two units was very low at around 18 s. However, preliminary tests with five units showed no convergence after more than four hours of computation. Indeed, although the solver can easily handle the 5376 binary variables in the two-CCGT problem, this is not the case with five units (13440 binary variables) or thirteen units (34944 binary variables) as in the Jurong Island model [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF]. The complexity of the developed model mainly lies in the introduction of multiple possible values for 𝛼 i,t,r ph_b . A first simplification could consist of reducing the size of the set R with fewer values for the power/heat ratios (e.g. {0, 0.3, 0.6, 1}). Instead, a constant ratio over the time horizon is considered here. The set R remains the same as previously (i.e. 11 possible values) and constraint ( 6) is rewritten as follows with r i ref the constant power/heat ratio of unit i over the day (with r an integer in {1,..,11}).

          otherwise 0 if 1 _ , , _ , , b ph r t i ref i b ph r t i T t r r    
A new set of simulations is performed while varying the references for the power/heat ratios of the two units. Also, the impact of the power load level is investigated for three cases: low (L at 50%), middle (M at 75 %) and high (H at 100 %). The heat load profile remains the same as previously and A bl = 500 kg/MWh. An exhaustive search is performed to find the best set of constant ratios (r 1 ref and r 2 ref ) that minimizes the emissions. That approach is denoted as S3. As already observed in the previous subsection, results obtained with S2 are better than in a case where the power/heat ratio of the unit remains at one (i.e. S1) (Table II). The improvements increase when the heat load is more significant compared to the electrical demand (i.e. moving from H to L here). The strategy S3 with constant ratios displays intermediate results that become closer to the performances of S2 with lower power loads. Fig. 8 plots the results obtained while varying ratios r 1 ref and r 2 ref for the different power load levels. For higher electrical demands the dispatch strategy cannot converge with low ratios. In such cases the combined power output of the two units is not enough to supply the load (Fig. 8b,c). Considering case L, the convergence is obtained whatever the chosen reference ratios are (Fig. 8a). Each of the plotted surfaces displays a global minimum that corresponds to the optimal set of power/heat ratios with the lowest CO2 emissions. The exhaustive search performed by S3 is possible in a case with two units and only requires 30 min with 121 different sets of ratios estimated (i.e. 11 2 ). However, computational times might become prohibitive when more generators are considered (e.g. 11 13 possible sets for the Jurong Island model). To overcome that difficulty a greedy search method is implemented and the MILP dispatch is successively run for different sets of power/heat ratios. The idea of the method is to iteratively decrease the ratios of the units until no improvement is possible in the objective (obj * ). At each iteration the method identifies the unit whose power/heat ratio should be decreased in order to reduce the emissions as much as possible (Table III). That is done by independently decreasing the ratio of each unit i and computing the corresponding emissions (obji). For the three different power load levels (i.e. low L, medium M and high H) the developed algorithm returns the same solution as the exhaustive search S3. However, it is obvious that the greedy approach cannot guarantee a global optimum in every case. The main advantage lies in a limited computational time with a maximum number of evaluations (i.e. runs of the MILP dispatch) directly linked to the number of units and the size of the set Rmaximum of R×I ×I evaluations. The previous subsection introduced an optimal power/heat dispatch in the presence of cogeneration units. A security constrained unit commitment is now considered with the inclusion of models for power and thermal networks. The objective is to take account of possible constraints regarding the flows within electrical cables and steam pipes.

1) Electrical Grid Model

For the power system a traditional DC power flow model is considered with the computation of Shift Factors (SF L E B ) [START_REF] Bardulescu | Software tool for power transfer distribution factor (PTDF) computing within the power systems[END_REF]. Then the power flow within a line l e at time t is computed by considering both positive and negative components F l e ,t e+ , F l e ,t e-as well as the power injection at each bus. A linear coefficient 𝛿 l e e for the branch losses is introduced together with the matrices MBI and M B𝐋 𝐄 e that represent the grid topology -Mbi = 1 if unit i is connected to bus b; M b𝑙 𝑒 e = 1 if line l e starts/ends at bus b. Finally, line power flow is computed according to [START_REF] Eggleston | IPCC guidelines for national greenhouse gas inventories[END_REF]. As in the method developed in [START_REF] Li | Fully reference-independent LMP decomposition using referenceindependent loss factors[END_REF] the losses within a line are equally distributed at the start and end buses as additional loads. The power balance constraint is modified following [START_REF] Van Den Bergh | A Mixed-Integer Linear Formulation of the Unit Commitment Problem[END_REF] to consider the balance at each bus b for every time step t. Finally, constraints (17) ensure that the line power flow remains below the specified limit F l e e_M .

 

                        

2) Thermal Grid Model

The thermal grid displays two main differences compared to the power system. Firstly, the steam flows F l h ,t h in the pipe l h are strictly unidirectional because of the irreversibility of turbomachinery [START_REF] Lund | 4th Generation District Heating (4GDH): Integrating smart thermal grids into future sustainable energy systems[END_REF]. Secondly, those steam flows can only take values below a predefined capacity F l h h_des within a certain range θ (typically 25 %) [START_REF] Wang | A novel model for steam transportation considering drainage loss in pipeline networks[END_REF] (( 18)). In order to ensure the convergence of the scheduling problem, controllable damp loads P b,t dp and boilers P b,t bl should be considered at each bus. In particular, they allow more flexibility if the units' operating conditions and the heat flow limits do not allow the supply of loads P b𝑙 ℎ h or if there is a local excess of heat generation. Finally the heat balance at each bus is expressed using [START_REF] Li | Fully reference-independent LMP decomposition using referenceindependent loss factors[END_REF] with M B𝐋 𝐇 h+ and M B𝐋 𝐇 h-B×L H matrixes that map the heat network -M b𝑙 ℎ h+ = 1 if line l h ending at bus b and M b𝑙 ℎ h = 1 if l h starting from bus b.

The equation also considers linear loss coefficients δ l h h that depends on the pipe length.

des h l h t l des h l h h h F F F _ , _ 100 1             
                   B. Obtained Results.

1) 5-bus Test Case

The power/heat dispatch problem with grid constraints is implemented for the 5-bus system displayed in Fig. 9. Parameters for power lines and maximum outputs of the generators are derived from [START_REF] Li | Small test systems for power system economic studies[END_REF] with technical parameters corresponding to the CCGT units in Jurong Island (five units A, B, C, D and E) . A thermal network is designed with the topology depicted in Fig. 9b with specified directions for the steam flows. The pipes capacities as well as the heat loss coefficients are given in Table IV. Three daily profiles are considered for the electrical loads at buses 2,3 and 4. The heat profiles are assumed to follow the same patterns, and different levels of power/heat demands can be investigated. First, an optimal scheduling is considered with a heat load equals to 150 % of the power demand profile and the boilers specific CO2 emission is set at 500 kg/MWh. A preliminary test is performed while running the optimal dispatch with all the power/heat ratios set to 1 for every unit (S1) and it corresponds to 9639 tons of CO2 emitted for the representative day. Then the hybrid Greedy/MILP dispatch (SGD) presented in the previous section is applied to minimize the CO2 emissions with adjusted r i ref (integers in {1,11} for power/heat ratio in [0,1] ). Emissions are significantly reduced to 8123 tons (16 % reduction) with the CCGT units generating more heat than in S1. As displayed in Fig. 10a the increase of the waste heat generation in SGD allows extensive use of the boilers to be avoids, resulting in reduced CO2 emissions. For low heat load levels, the CCGT heat production is even greater than the demand. The use of damp loads is then required, which is not necessary in S1. During peak load, hours the heat generation remains slightly greater than the load as it feeds the system heat losses (Fig. 10a). The convergence of the greedy method is reached after 45 min of computation with 34 iterations, which corresponds to 169 evaluations of the objective function (i.e. runs of the MILP dispatch with constant ratios). The optimal solution returned displays relatively low values for the power/heat ratios compared to the previous case study with two units (unit A: 0.9, unit B: 0.1, unit C: 0, unit D: 0.6 and unit E: 0.1). That is explained by the electrical power load being significantly lower than the installed capacity: the peak load is 620 MW for 1.2 GW of total capacity. There is no problem for the CCGT units to operate at lower power/heat ratios while supplying the electrical demand. Thus those ratios decrease along the iterations of the greedy search method as well as the CO2 emissions until no improvement is possible (Fig. 10b). Fig. 11 displays the line flows normalized by the maximum capacities. The powers within the electrical lines are not important (below 50 % of load) as the electrical demand level is low compared to the system capacity (Fig. 11a). Regarding the steam pipes, the lower bound for the heat flows appears to be a binding constraint (Fig. 11b). Indeed, except for line 1-2 all the flows remain close to the minimum value which then requires the generation of additional heat (either with the boilers of the CCGT units).

An additional set of simulations is performed, varying the heat load level (in % of power demand) as well as relaxing the steam flow constraint with greater values for θ (which is not realistic). The CO2 emissions are computed for the two scenarios S1 and SGD and the obtained results are compiled in Table V. As already observed in Section IV the improvements provided by an optimal dispatch compared to S1 tend to increase with greater levels of heat demand. The relaxation of the thermal flow constraints also allows greater CO2 reductions in all cases. The obtained results confirm that the constraint is a binding limit for all the considered load levels and while applying both dispatch S1 and SGD. Thus the design capacity of the pipe should be appropriately determined in the planning phase when dealing with cogeneration problems [START_REF] Li | Transmission-Constrained Unit Commitment Considering Combined Electricity and District Heating Networks[END_REF]- [START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF]. On the one hand, undersized pipes would limit the use of the heat generation capacity of the units. On the other, oversized lines might require additional boilers to operate under low load condition to maintain a minimal flow in the branches. Final set of simulations refers to the model of Jurong Island. The power grid is represented by a 208 buses, 219 branch model and the thermal network displays 14 nodes and 14 steam pipes. The distinction is made between purely electrical customers and nodes with both thermal and electrical loads (Fig. 12a). Heavy customers (e.g. oil refineries) are assumed to get their own steam supply while the other thermal loads are aggregated into mutualized heating districts (Fig. 12b). Losses are computed based on the cable/pipe characteristics and the estimated distances between the buses. Simulations are run for A bl = 500 kg/MWh and θ = 25 %. Three levels of heat demand are investigated (i.e. low L, medium M and high H). Similar to Section IV, S1 corresponds to an optimal dispatch with all the power/heat ratios set to 1 for every unit and SGD to the application of the hybrid Greedy/MILP dispatch. Results obtained in the various cases are displayed in Table VI. The computational time to run SGD remains significant even with the use of the greedy search from 3 h for the case L to more than 5 hours for case H, for the representative day. As already observed the improvements are more significant for higher heat load levels. Following the remarks of the previous subsection, better results could be obtained by investigating different steam pipe capacity configurations. Additionally, not all the units might have to be considered in the power/heat dispatch. Indeed, for the three investigated load levels, only three to four CCGT actively participate in the steam supply with optimal power/heat ratios lower than 1.

VI. CONCLUSIONS

This paper successfully implemented a CCGT model for power and heat cogeneration as well as an efficient power dispatch strategy based on the minimization of CO2 emissions. The obtained results show significant improvements compared to a case where cogeneration units are not controllable (i.e. priority given to electrical output). Those improvements become greater when higher heat load levels or dirtier additional boilers are considered. The observations also showed the necessity to appropriately size the steam pipes as lower bounds for heat flows are binding constraints in most cases. Further work may investigate design strategies. Also, other cost components should be included in the dispatch objective, in addition the CO2 emissions that have been considered so far.

VII. ACKNOWLEDGMENT

  power/heat load at bus b at time t Ai,c slope for carbon cost of unit i in block c C0i, SUi base cost and start-up cost of unit i p i,c gt_M upper bound of block c for GT of unit i lines loss coefficients θ range of heat flow in steam pipes (in %) SF L E B L E ×B matrix with generation shift factors M B𝐋 𝐄 e B×L E matrix that maps the power network M BI B×I matrix that maps the units in the grid M B𝐋 𝐇 h+ , M B𝐋 𝐇 h-B×L H matrixes that map the heat network II. INTRODUCTION

Fig. 1 :Fig. 2 :

 12 Fig.1: CCGT normal operation -a) operating variables -b) performances[START_REF] Zhang | Optimal Allocation of CHP-Based Distributed Generation on Urban Energy Distribution Networks[END_REF]

Fig. 3 :

 3 Fig. 3: CCGT power/heat operationa) power outputb) heat output

Fig. 4 :Fig. 5 :

 45 Fig.4: Example of unit i operating in block s1 with ratio r1 at time t

Fig. 7 :

 7 Fig. 7: CO2 Vs the heat load levela) A bl =300 -b) A bl =500 -c) A bl =700

Fig. 8 :

 8 Fig. 8: Obtained results with S3a) level Lb) level Mc) level H

  of Losses and Grid Constraints.

Fig. 9 :

 9 Fig. 9: 5 bus systema) electrical networkb) thermal network

Fig. 10 :

 10 Fig. 10: Obtained resultsa) heat load supplyb) greedy search method

Fig. 11 :

 11 Fig. 11: Branch flowsa) power linesb) steam pipes

Fig. 12 :

 12 Fig. 12: Jurong Islanda) power gridb) thermal network

TABLE VI

 VI

	CO2 emissions for different heat load levels
	Heat load level	S1	SGD
	L	33385	32201
	M	41757	39415
	H	56558	53801

This work was supported by the Singapore National Research Foundation (NRF) under its Campus for Research Excellence and Technological Enterprise (CREATE) programme, specially the Cambridge Centre for Advanced Research in Energy Efficiency in Singapore (CARES).


Laurent Larrieu

The " Groupement Forestier des montagnes particulières de Hèches ": a case from the mountains of southwestern France

The « Groupement Forestier des montagnes particulières de Hèches »: a case from the mountains of southwestern France Laurent Larrieu An property owned by 100 people and collectively managed

INRAE

  livestock issues • 1 for timber issues A chairman elected by the owners A 10-members board which meets four times a year Once a year, a meeting gathers all the members of the ownership group and the main strategic guidelines are discussed and then approved.

About 6000 m3/y of timber harvested + 120 m3/y of fuelwood for the owners