V Hamm

C Kervévan
email: c.kervevan@brgm.fr

D Thiéry

CO 2 -DISSOLVED: a Novel Concept Coupling Geological Storage of Dissolved CO 2 and Geothermal Heat Recovery -Part 4: Preliminary Thermo-Hydrodynamic Simulations to Assess the CO 2 Storage Efficiency

Keywords: CO2 storage, dissolved CO2, geothermal energy, saline aquifer, thermo-hydrodynamic modeling, storage efficiency

The CO 2 -DISSOLVED project aims at assessing the technical-economic feasibility of coupling CO 2 storage in a saline aquifer and geothermal heat recovery. The proposed infrastructure basically relies on a standard geothermal doublet where CO 2 would be injected after having been entirely dissolved in the cooled brine. The objective of this preliminary thermo-hydrodynamic modeling study was to quantify the expected CO 2 storage lifetime and efficiency. The results first confirmed that CO 2 will inevitably be produced in the extracted brine after 2-15 years of continuous injection, depending on the operating parameters. However, mass balance calculations evidenced that after a 30 year injection period, 37-85% of the total CO 2 injected should remain stored in the aquifer.

Introduction

This study was conducted in the framework of the CO 2 -DISSOLVED project [START_REF] Kervévan | CO2-DISSOLVED: a novel approach to combining CCS and geothermal heat recovery[END_REF], funded by the ANR (French National Research Agency). This project proposes to assess the feasibility of a novel CO 2 injection strategy in deep saline aquifers, combining injection of dissolved CO 2 (instead of supercritical CO 2 ) and recovery of the geothermal heat from the extracted brine. This approach relies on the closed-loop geothermal doublet technology where the warm water is extracted at the production well and the cooled brine is re-injected in the same aquifer via a second well (injection well).

The objective of the work presented here was to provide first quantitative elements to assess the feasibility of the concept in terms of storage lifetime and efficiency. For that purpose, we performed hydrodynamic modeling of several test-cases relying either on the conventional doublet technology (deviated injection and production wells), or on a set of three wells (one injector / two producers, or conversely). A third case was also investigated using a horizontal well architecture instead of a typical deviated well for the injection. The selected hydrodynamic parameters are typical of the deep geothermal resources of the Paris basin (Dogger aquifer). Several distances between wells and production flow rates were simulated. In those preliminary simulations, the dissolved CO 2 was considered as a non-reactive tracer with no impact on the CO 2 laden brine properties (viscosity and density), and no chemical reactivity. We also performed some simulations accounting for coupling with heat transfer (dependency of brine viscosity and density upon temperature).

The first part of this paper describes the conceptual model on which our simulations were based on. The second part of the paper deals with the results achieved for the various scenarios simulated. The results presented focus more specifically on both the general hydrodynamic behavior of the system, and the storage efficiency defined here as a mass balance ratio between the CO 2 actually stored in the aquifer and the total CO 2 injected at the end of the injection period.

Modeling approach

Nomenclature

Variables

V th thermal front velocity (L.T -1 ) ω porosity (-) ρ A C A Aquifer global heat capacity (J.L -3 .K -1 ) ρ F C F fluid heat capacity (J.L -3 .K -1 )
Superscripts ___ rate of the corresponding variable, with respect to time (e.g. mass rate) Subscripts inj at the injection well prod at the production well stored stored in the aquifer

Conceptual hydrodynamic model

The conceptual model relies on a 2D single-layer geometry for the flow and CO 2 non-reactive tracer transport calculations, and on a pseudo 3D multilayer model for the coupled flow, CO 2 and heat transport calculations. The finite volume code MARTHE, developed at BRGM [START_REF] Thiéry | Groundwater Flow Modeling in Porous Media Using MARTHE[END_REF][START_REF] Thiéry | Presentation of the MARTHE 7.4 code developed[END_REF], has been used for the simulations. The reservoir model is based on homogeneous hydraulic and thermal properties. The horizontal dimensions of the domain are 20 x 10 km divided in rectangular cells of 100 x 100 m, with a horizontal refinement of 20 x 20 m in the vicinity of the wells (Fig. 1). The layer thickness is 20 m, which is representative of the cumulated average thickness of the productive layers in the Dogger aquifer. In the case of heat transport calculation with dependency of fluid viscosity and density upon temperature, the aquifer layer is refined and both base and cap-rock layers were added. Heat conduction in base and cap-rock is solved in z direction using a semi-analytical solution developed by Vinsome and Westerfeld [START_REF] Vinsome | A simple method for predicting cap and base rock heat losses in thermal reservoir simulators[END_REF] added to 3D numerical heat calculation (conduction plus advection) in the aquifer layers. The parameters used for the thermo-hydrodynamic simulations are summarized table 1. They are typical measurements as flow meter log, well testing, plug measurements representative of the Dogger aquifer.

CO 2 mass balance calculation

In order to quantify the amount of CO 2 effectively stored in the aquifer at the end of a simulated injection period, we first applied a general mass balance equation according to which:

(

Introducing more familiar variables such as concentration ( ) and volumetric flow rate ( ), equation (1) can be rewritten as: [START_REF] Thiéry | Groundwater Flow Modeling in Porous Media Using MARTHE[END_REF] In these simulations, we considered that no chemical reaction occurred. As a consequence, it can be reasonably assumed that the source/sink term in the above equations can be omitted ( ⁄ ). Hypothesizing a constant flow rate all along the injection period ( ), we can then easily infer from equation (2) the total mass of CO 2 accumulated in the aquifer ( ):

∫ ∫ ( ) (3)
C inj is considered as a constant term in our simulations, and C prod is an output of the hydrodynamic model that is tabulated at every time step, so that the integral can be easily calculated using classical integration algorithms included in general purpose numerical computation software (we used a short Scilab script to calculate the integral in equation ( 3)).

Results and discussion

All the simulations presented here assumed no chemical interaction between the brine and the porous matrix of the aquifer. Similarly, dissolved CO 2 was assumed to behave as a non-reactive tracer; in particular, the effect on brine density and brine viscosity were not accounted for. Only the temperature effect on water properties (viscosity, density) was accounted for in the simulations considering heat transfer.

Consequently, in all the calculations, the CO 2 concentration in the injected water was assumed to be constant all along the injection period; the absolute value of this concentration was then meaningless, that is why most of the results involve a dimensionless relative concentration comprised between 0 and 1 (ratio of the CO 2 concentration in water to the CO 2 concentration in the injected water). For the mass balance calculations however, we assumed a CO 2 concentration in the injection water of 50 g.L -1 (approximately 1 mol.L -1 ) which is a reasonable order of magnitude for CO 2 solubility in the typical (P, T, S) conditions of the exploited geothermal aquifers of the Paris basin (e.g., see fig. 4 in [START_REF] Kervévan | Improvement of the Calculation Accuracy of Acid Gas Solubility in Deep Reservoir Brines: Application to the Geological Storage of CO2[END_REF]).

Single doublet, without heat transfer

This first series of simulations aimed at quantifying the effect of the key operating parameters of the doublet system, i.e. the distance between the well shoes, and the injection flow rate. Three typical distances between wells (1, 1.5, and 2 km) and flow rates (100, 200, and 300 m 3 .h -1 ) were investigated, all being typical of actual operating parameters of the geothermal doublets of the Paris basin. One of the primary questions to answer was about the quantification of the CO 2 migration time in the aquifer between the injection and the production wells. In a doublet system, the fact that an injected product is going to be present in the extracted brine after a while (at increasing concentration with time in case of continuous injection), is inevitable. In our simulations where injection is continuous, CO 2 is considered as a tracer. Consequently, it is expected to be present at the production well after a certain period to be determined by our calculations.

The results provide quantitative estimates of the concentration of the CO 2 produced at the extraction well and on the timeframe of the mass transport process. Unsurprisingly, the curves on Fig. 2 clearly show, at a given distance between wells (1.5 km in this case), that the higher the flow rate, the shorter the CO 2 breakthrough time at the production well, and the higher the CO 2 concentration in the extracted brine. It takes approximately 2 to 8 years in this example to get a 0.1 relative concentration. On Fig. 3, the effect of the distance between wells is evidenced: the higher it is, the higher the breakthrough time is and the lower the CO 2 concentration is. For a distance of 4 km between the wells, the threshold concentration value of 0.1 is not reached at the end of a 30 year continuous injection period.

Consequently, for an efficient CO 2 storage using a doublet system, it seems at first glance that we should have to minimize the flow rate and maximize the distance between wells. However, simple technical reasons prevent us from choosing the distance between wells we would like to have: for instance, in the geothermal reservoirs of the Paris basin, the depth of the productive layers of the Dogger aquifer is, on average, between 1,500 and 2,000 m. With the constraint of a typical distance between the well heads of a few meters (for minimal surface footprint and heat loss), and a limit of 55° in the deviation angle, the maximum distance between the well shoes cannot exceed 1.5 to 2 km. So, the distance of 4 km between wells presented on Fig. 3 is actually an ideal case but such a doublet could not be implemented in the context of the currently exploited geothermal aquifer of the Paris basin. On the other hand, decreasing the flow rate would decrease the CO 2 mass injected as well as the quantity of geothermal energy recoverable. However, the most important point to consider when studying a system basically dedicated to CO 2 storage is to determine how much CO 2 will be effectively stored at the end of an injection period. Based on the breakthrough curves presented previously, it is demonstrated that CO 2 is going to be produced at the extraction well after a few years of injection. Even though this produced CO 2 is not going to be released to the atmosphere because the injection-production wells work in closed-loop, these results strongly suggests that part of the injected CO 2 is not going to remain stored in the aquifer. In terms of mass balance however, equation (3) demonstrates that since C prod remains lower than C inj (i.e. relative concentration remains strictly lower than 1), the mass of CO 2 stored is then necessarily greater than 0; so, CO 2 storage is going to be effective.

In order to quantify more precisely the amount of stored CO 2 , we first reworked the data of Fig. 2 to express them in terms of storage rate (Fig. 4). The results confirm the general trend observed on Fig. 2 according to which the lower flow rate provides a less and later significant decrease in the storage rate. However, it is important to remark that, even though a very significant decrease in the storage rate is observed relatively early (after 2-3 years of injection) for the highest flow rate of injection (300 m 3 .h -1 ), this storage rate always remains greater than for the lower flow rates. Now, performing a numerical integration of equation ( 3) for all the operating parameters, we were able to calculate the cumulated masses of stored CO 2 over the 30 year injection period considered in these simulations. The results are synthesized on Fig. 5. They clearly illustrate the impact of both the distance and the flow rate on the accumulated mass of CO 2 effectively stored in the aquifer: for a given flow rate, the accumulated mass increases with the distance; for a given distance, the stored mass increases with the flow rate.

For a better interpretation of these results, it is interesting here to introduce the notion of storage efficiency that we define, for a given injection period, as the percentage of the mass of CO 2 effectively stored in the aquifer with respect to the mass injected. On Fig. 6, the results expressed in terms of storage efficiency provide a new insight on the understanding of the storage system: to maximize the efficiency of the storage, it is necessary to work with the lowest flow rate (100 m 3 .h -1 ) and the highest distance between the wells shoes (2 km); however, with these parameters, the absolute value of the total mass of CO 2 stored after 30 years (1,118 kt) is surprisingly lower than for the least efficient case (Q = 300 m 3 .h -1 ; d = 1 km) where the mass of CO 2 stored reaches 1,447 kt.

At this stage however, it is too early to draw definitive conclusions on the best strategy to select: optimizing the storage efficiency or maximizing the total amount of stored CO 2 . Other aspects have to be considered first: optimization of the geothermal energy recovered as a function of the needs, expected lifetime of the doublet, configuration of the site (surface and subsurface), economic aspects, etc. This will be treated later in the project, when application test cases will have been selected.

Single doublet with heat transfer

In this second set of simulations, we simulated the coupling of heat transport and non-reactive tracer CO 2 mass transport by taking into account brine density and viscosity dependence on temperature (thermal coupling). The conditions used for these simulations are a constant flow rate of 200 m 3 .h -1 , a distance between wells of 1,500 m and a re-injected fluid temperature of 40 °C. The initial fluid temperature was set to 70°C which is the average temperature of the Dogger aquifer in the Paris basin.

Figure 7 shows the concentration of the CO 2 produced at the extraction well on the left axis versus the production temperature on the right axis. There is a retardation factor of about 4 between the thermal breakthrough and the tracer breakthrough. This retardation factor can be expressed by equation ( 4): [START_REF] Vinsome | A simple method for predicting cap and base rock heat losses in thermal reservoir simulators[END_REF] Consequently, after 30 years of geothermal exploitation, the production temperature drop is less than 0.5°C, which is not significant compared to the relative CO 2 concentration variation at the same date. Furthermore, the comparison between the growth of the CO 2 plume and of the thermal plume show a great difference of the domain impacted (Fig. 8). This process is due to the important role of the base and cap rock that slow down the propagation of the thermal front. The effect of temperature on brine properties and consequently on the tracer velocity is limited and results only in a slightly higher breakthrough time and lower relative concentrations at the production well. Indeed, as the doublet flow rate is maintained constant (prescribed flow rate), an increase in fluid viscosity around the injector due to the cooling implies a slight increase in pressure gradient between wells with constant velocity. The effect of temperature upon density (at constant salinity) is negligible for the thickness considered. In a later stage, the effect of dissolved CO 2 on fluid viscosity and density should be considered. 9 and 10, based on the same principle as in §3.1, show respectively the evolution of the stored CO 2 rate with time, and the total mass of CO 2 stored and the storage efficiency after a period of 30 years for both cases (with and without thermal coupling). These figures confirm that the change in temperature does not change much the stored CO 2 rate and, as a consequence, the total mass stored after 30 years of injection of a cold CO 2 -rich solution . In these examples, it then appears that the temperature of the injected fluid has less influence for the CO 2 storage than for the geothermal heat exploitation.

Other wells configurations

In this last section, with an objective of improving the storage efficiency, we tested some other wells configurations such as a set of three wells (triplet) instead of a doublet or the use of a horizontal or sub-horizontal well for injection. The horizontal well technology is already used in oil or gas industry extraction and a first case is envisaged for both production and injection in the Dogger aquifer for geothermal purpose. In the next figures, the notation used are "IP" for a classical doublet with one injection and one production well, "PIP" or "IPI" in the case of an alignment of three wells with one central injector and two producers, or two injectors surrounding a central producer, respectively. The "IIP" configuration, with one producer and two injectors grouped together was also tested hereafter. The conditions for the simulations are identical to the §3.1 case, i.e. no thermal effect, a constant total flow rate of 200 m 3 .h -1 (divided evenly between the two injectors or producers), and a distance between each wells of 1,500 m. This distance is maximised with respect to the wells configuration (one central vertical well and two deviated wells). Nevertheless, the drilling of recent new wells in the Dogger aquifer is currently using two deviation angles at different depths leading to an increase in the distance between well shoes. In the case of a horizontal drain for the injection, the length is considered to be 1,000 m.

Figure 11 compares the CO 2 relative concentration at the production well for a classical doublet, three wells configurations, and 2 wells with one horizontal injection well. In the PIP or IPI cases, the CO 2 breakthrough time is delayed compared to the case of a doublet, but the concentration that reaches the production well(s) becomes higher than the one observed for a doublet after 10 years of production. In the IIP case, the CO 2 breakthrough time is also greater than for a doublet but the CO 2 relative concentration at the production well remains always significantly smaller than the one calculated for a doublet. Finally, in the case of a horizontal injection well, the CO 2 breakthrough time is the most delayed, but after 30 years of production, the CO 2 relative concentration remains smaller than for a classical doublet, but higher than for the IIP scheme. 12 to 14 show the corresponding stored CO 2 rates and the total stored CO 2 and storage efficiency over an injection period of 30 years. They indicate that the IPI or PIP schemes does not improve the CO 2 storage efficiency and even slightly degrade it. In contrast, both the IIP and the horizontal injector schemes significantly improve the CO 2 storage efficiency (respectively 68.4% or 67.1% instead of 55.9% for a classical doublet).

Nevertheless, we have to keep in mind that we can also reach those rates with a classical doublet by decreasing the flow rate and/or increasing the wells distance (fig. 6). The benefit of those technical solutions using three wells or horizontal wells must also be balanced with their higher initial capital cost and the actual energy needs, in order to define an optimal injection scenario.

Conclusion

This series of preliminary thermal-hydrodynamic simulations provides a first overview of the influence of the operating parameters and the wells configuration on the potential lifetime and efficiency of the CO 2 storage system. The results acquired so far demonstrate that:  Dissolved CO 2 inevitably reaches the production well after a while. Typical breakthrough times obtained here range between 2 and 15 years for a technically achievable doublet in the context of the Dogger aquifer of the Paris basin. However, it has to be noticed that this CO 2 will not be released to the atmosphere since the production/heat exchanger/injection circuit is a closed loop. Increase in CO 2 concentration in the produced brine will then result in a lower quantity of gaseous CO 2 that could be dissolved at the injection well.

 The storage efficiency, defined here as the mass ratio between stored and injected CO 2 , decreases with the flow rate, while it increases with the distance between the wells.

 For all the injection scenarios considered in our simulations, even the less favorable, a significant ratio of the injected mass of CO 2 , varying from 37% to 85%, remains stored in the aquifer at the end of an injection period of 30 years (typical duration of an operating permit for a geothermal doublet in the Paris basin).

 Using a combination of two injection wells and one production well (triplet) or two wells with one horizontal injection well instead of a classical doublet, improves the storage efficiency (all equivalent operational parameters being equal). However, the required infrastructure would then be obviously more expensive; economic calculations applied to a well-defined test-case are necessary to conclude on the viability of those options.

In the next tasks of the project, more accurate simulations will be performed, accounting for both the effect of dissolved CO 2 on fluid viscosity and density, and the water/rock interaction processes. The combined effects of those chemical-physical processes is expected to slow down the migration velocity of dissolved CO 2 toward the production well of the doublet (or triplet). If this behavior is confirmed by numerical modeling, better storage efficiency will then be obtained, corroborating the conclusions of this study on the viability of the CO 2 -DISSOLVED concept for achieving dissolved CO 2 storage in a saline aquifer exploited by a geothermal doublet (or triplet).

  C concentration of dissolved CO 2 in the aquifer (M.L -3 ) d distance between the wells of the doublet, bottom hole (L) dm/dt accumulation term in the dissolved CO 2 mass balance equation (M.T -1 ) dr/dt source/sink term in the dissolved CO 2 mass balance equation (M.T -1 ) m mass of dissolved CO 2 in the aquifer (M) P pressure (M.L -1 .T -2 ) Q CO 2 laden brine volumetric flow rate of injection (L 3 .T -1 ) S salinity (M.L -3 ) T temperature (K)  injection period (T) R thermal retardation factor (-) V e effective pore velocity (L.T -1 )

Fig 1 .

 1 Fig 1. Top view of the meshing and different wells locations (I: injection well, P: production well at different distance from the injector (P1: 1 km, P2, 1.5 km or P3: 2 km

Fig. 2 .

 2 Fig. 2. CO2 breakthrough at the production well for continuous injection at various flow rates. The injection and production well shoes are located 1.5 km apart (P2 in Fig 1).

Fig. 3 .

 3 Fig. 3. CO2 breakthrough at the production well for continuous injection with various distances between the injection and the production well shoes (P1, P2 or P3). The injection flow rate is constant at 200 m 3 .h -1 .

Fig. 4 .

 4 Fig. 4. Evolution of the stored CO 2 rate all along a 30 year injection period, as a function of the injection flow rate (Q). The injection and production well shoes are located 1.5 km apart.

Fig. 5 .

 5 Fig. 5. Comparison between the cumulated masses (in kt) of stored and injected CO2 over a period of 30 years of continuous injection, as a function of the distance between wells (d) and the injection flow rate (Q).

Fig. 6 .

 6 Fig. 6. Storage efficiency over a period of 30 years of continuous injection, as a function of the distance between wells (d) and the injection flow rate (Q).

Fig 7 .

 7 Fig 7. CO2 and thermal breakthrough at the production well for continuous injection

Fig. 9 .

 9 Fig. 9. Evolution of the stored CO2 rate all along a 30 year injection period

Fig. 11 .

 11 Fig. 11. CO2 breakthrough at the production well for continuous injection and different wells configurations

Fig.

  Fig.12to 14 show the corresponding stored CO 2 rates and the total stored CO 2 and storage efficiency over an injection period of 30 years. They indicate that the IPI or PIP schemes does not improve the CO 2 storage efficiency and even slightly degrade it. In contrast, both the IIP and the horizontal injector schemes significantly improve the CO 2 storage efficiency (respectively 68.4% or 67.1% instead of 55.9% for a classical doublet).Nevertheless, we have to keep in mind that we can also reach those rates with a classical doublet by decreasing the flow rate and/or increasing the wells distance (fig.6). The benefit of those technical solutions using three wells or horizontal wells must also be balanced with their higher initial capital cost and the actual energy needs, in order to define an optimal injection scenario.

Fig. 12 .

 12 Fig. 12. Stored CO2 rate for different wells configurations

Table 1

 1

	: Model parameters
	Parameters	Values
	Rock porosity	0.15 [-]
	Specific storativity	10 -6 [m -1 ]
	Average aquifer transmissivity	40 [Dm]
	Average aquifer temperature	70 [°C]
	Fluid density	0.4 [mPa.s]
	Average fluid salinity	20 [g/l]
	Fluid heat conductivity	0.6 [W/m/°C]
	Fluid heat capacity	4.18 [MJ/m 3 /°C]
	Rock heat conductivity	2.5 [W/m/°C]
	Rock heat capacity	2.2 [MJ/m 3 /°C]
	Longitudinal heat diffusivity	20 [m]
	Transversal heat diffusivity	10 [m]

Acknowledgements

We would like to thank both the ANR (French National Research Agency) and the Research Directorate of BRGM for funding this study in the framework of the CO 2 -DISSOLVED project. We are also grateful to the project team members for allowing us to publish these results.


Juliette Heintz

Jean-Pierre Belaud

Nishant Pandya

Moises Teles

Dos Santos

Vincent Gerbaud
email: vincent.gerbaud@ensiacet.fr

Computer aided product design tool for sustainable product development

Keywords:

Molecular graph Bio-based molecule Sustainable product design

A computer aided product design (CAPD) tool is proposed that finds mixtures matching target properties. Genetic algorithm crossover and mutation operators are completed with insertion or deletion operators adapted for side branches. A new substitution operator is devised for cyclic molecules. The mixture fitness is evaluated by a weighted sum of property performances. Molecules are represented by molecular graphs. They are split into molecular fragments which are built from polyatomic groups. Molecules or molecular fragments can be fixed, constrained or left free for building a new molecule. Building blocks are chemical functional groups or bio-sourced synthons. A specific coding of hydrogen-suppressed atoms is devised that can be used with various property estimation models where atom connectivity information is required. Illustration is provided through three case studies to find levulinic, glycerol and bio-based derivatives as substitute for chlorinated paraffin, methyl p-coumarate ester solvent and blanket wash solvent, respectively.

Introduction

The chemical industries are on the frontline of sustainable development due to the potential impact on the environment, health and safety of its product and process activities. Regulations such as the European REACH (REACH, 2006) and VOC [START_REF] Voc | VOC Solvents Emissions Directive[END_REF] directives or the keen interest of consumers for eco-labelled products push the chemical industries to reconsider the products which they use and produce.

In Europe, the cost of registering chemicals to comply with REACH could exceed D 2.1 billion, based on about 30,000 substances (ECHA, 2012). Therefore there is a strong incentive to find substitute molecules and chemical products. New products need to obey environmental, health and safety constraints in addition to usual product and process requirements. Economists have argued that a doubly green chemistry perspective prevails among chemical industry engaged in green activities: one green for the reduction of their impacts on environment and one green for the use of renewable raw materials [START_REF] Garnier | The emergence of doubly green chemistry, a narrative approach[END_REF]. The first perspective is a direct transcript of the definition of sustainable growth in the founding Brundtland 1987 report. The second is the seventh principle of green chemistry [START_REF] Anastas | Green chemistry theory and practice[END_REF]. As it should allow sustainable issues like toxicity or degradability to be met more easily, the use of bio-sourced molecules or synthons is a major stimulus when looking for a new product.

For finding a substitution molecule, the usual 'trial and error' approach seems inefficient unless high throughput screening is used. Instead, reverse engineering approaches, like Computer Aided Molecular Design (CAMD) are fit to handle several properties and to propose molecular structures matching the target values of these properties. In some cases, the problem of substituting a molecule may result in proposing a mixture. This further brings forth the challenge of computing mixture properties which may not always obey a linear mixing rule.

This paper presents a Computer Aided Molecular Design tool and its tailoring for finding alternative bio-sourced molecules and mixtures, with the help of model driven engineering (MDE) concepts. The CAPD tool follows the general methodology of CAMD tool with several modifications. By using a genetic algorithm, this tool simultaneously optimizes the molecular structure of the components and their compositions in the mixture in order to best fit the desired properties at normal operating conditions set by the user.

After a section devoted to present background information related to CAMD, we describe the data structures and methods. They concern molecular representation, atom coding, fragment builder along with specific genetic operators to build or delete side chemical branches and to enhance changes in aromatic rings while keeping their aromaticity. Their implementation into a three software-component tool is then presented using MDE concepts. Three case studies are presented to illustrate some of the features of the tool: mixture search (case 1), search of a molecule with predefined bio-sourced synthons (case 2), two level search (case 3).

Background

Computer Aided Molecular Design (CAMD) aims at finding molecules that satisfy a set of property targets defined in advance [START_REF] Achenie | Computer aided molecular design: theory and practice[END_REF]. CAMD relies upon four main concepts, namely, a molecular representation model, a set of property calculation models, a solving method and a performance criterion. Candidate molecules can be searched in a database or built from chemical groups. Their fitness is evaluated thanks to property estimation models by comparing the values of estimated property and the target property. Then they are discriminated according to their performance and either modified, kept as is or rejected, with the help of the solving algorithm. During the problem setting, in addition to the initial definition of the property target values, chemical blocks are pre-selected to be used in the molecular construction.

The CAMD problem solving method has often been tailored to a specific representation model. The early "generate and test" method was developed for a set of chemical groups that were also used by the group-contribution property estimation method [START_REF] Gani | A group contribution approach to computer-aided molecular design[END_REF][START_REF] Constantinou | Computer aided product design: problem formulations, methodology and applications[END_REF]. A vector of groups and their occurrences described candidates. However, a single vector may correspond to several isomer molecules and in this case a final step is required to generate the true molecules. To overcome this, some representation describing explicitly the group interconnections have been used: a genetic algorithm with adapted operators was used to generate polymers with a symbol string encoding [START_REF] Venkatasubramanian | Computer-aided molecular design using genetic algorithms[END_REF], a binary representation of atom connectivity in molecules was used with a MILNP method [START_REF] Churi | Novel mathematical programming model for computer aided molecular design[END_REF], an adjacency matrix was used with a simulated annealing [START_REF] Ourique | Computer-aided molecular design with simulated annealing and molecular graphs[END_REF], a graph representation was used with TABU search [START_REF] Lin | Computer-aided molecular design using Tabu search[END_REF] and recently a graphbased representation issued from signature descriptors was used with a genetic algorithm [START_REF] Herring | De novo molecular design using a graph-based genetic algorithm approach[END_REF]. These explicit representations of molecule are fit for many kinds of property estimation methods once a routine for finding the groups or descriptors of the corresponding estimation method is provided.

Regarding the fitness of a candidate molecule, the differences between the predicted and target values of all properties are aggregated in a global objective function through either an arithmetic mean [START_REF] Vaidyanathan | Computer-aided synthesis of polymers and blends with target properties[END_REF] or a geometric mean [START_REF] Castillo | Modified desirability functions for multiple response optimization[END_REF]. The geometric mean penalizes severely the fitness when an individual property prediction/estimation method is far from target. In that way it is more discriminant than the arithmetic mean.

The evaluation of the performance of each molecule relies upon the calculation of property values that have been classified as product-properties, process-related properties and usage-related properties [START_REF] Costa | Chemical product engineering: an emerging paradigm within chemical engineering[END_REF]. Product attributes found desirable or undesirable by consumers belong to the latter class. For the CAMD problem, product requirements have to be translated into target property values, which have been done by using problem templates (Mattei et al., 2014a,b). Most product and process properties are usually described by group contribution methods [START_REF] Joback | Estimation of pure-component properties from groupcontributions[END_REF][START_REF] Constantinou | New group contribution method for estimating properties of pure compounds[END_REF][START_REF] Martin | Prediction of the acute toxicity (96-h LC50) of organic compounds to the fathead minnow (Pimephales promelas) using a group contribution method[END_REF]Marrero andGani, 2001, 2002;[START_REF] Nannoolal | Estimation of pure component properties: Part 1. Estimation of the normal boiling point of non-electrolyte organic compounds via group contributions and group interactions[END_REF][START_REF] Nannoolal | Estimation of pure component properties: Part 2. Estimation of critical property data by group contribution[END_REF][START_REF] Hukkerikar | Group-contribution+ (GC+) based estimation of properties of pure components: improved property estimation and uncertainty analysis[END_REF] or QSAR/TI topological index/QSPR methods [START_REF] Veith | Structure-activity correlations in studies of toxicity and bioconcentration with aquatic organisms[END_REF][START_REF] Karelson | Quantum-chemical descriptors in QSAR/QSPR studies[END_REF][START_REF] Gani | Automatic creation of missing groups through connectivity index for pure-component property prediction[END_REF][START_REF] Chemmangattuvalappil | A novel methodology for property-based molecular design using multiple topological indices[END_REF]. Some environmental, health and safety (EHS) properties like R-phrase or CMR classification are described by similarity methods, relying upon the finding of specific molecular patterns in molecules [START_REF] Gallenos | Mini-review on chemical similarity and prediction of toxicity[END_REF].

The problem of designing a mixture is referred to as Computer Aided Product Design (CAPD) where individual molecules within the mixture and their composition must be found. Some CAMD methods have been extended to CAPD with an additional composition search [START_REF] Klein | Computer aided mixture design with specified property constraints[END_REF][START_REF] Gani | Computer aided molecular and mixture design with specified property constraints[END_REF][START_REF] Vaidyanathan | Computer-aided synthesis of polymers and blends with target properties[END_REF][START_REF] Duvedi | On the design of environmentally benign refrigerant mixtures: a mathematical programming approach[END_REF][START_REF] Churi | The optimal design of refrigerant mixtures for a twoevaporator refrigeration system[END_REF][START_REF] Sinha | A systematic method for integrating product attributes and molecular synthesis[END_REF]. Overall, CAPD raises new issues compared to CAMD: firstly, more properties have to be matched, including more usage-related product properties or the mixture stability. Secondly, several mixture property models such as boiling point and flash point, exhibit nonlinear mixing rules and need to be solved with built-in routines, which may increase the computation time. Thirdly, some usagerelated properties may not be described by any suitable prediction model.

Several approaches have been taken to solve CAPD problem: some have performed a sequential search of each mixture components individually, before checking mixture properties, stability and composition [START_REF] Gani | Chemical product design: challenges and opportunities[END_REF]Conte et al., 2011;[START_REF] Papadopoulos | Molecular design of working fluid mixtures for organic Rankine cycles[END_REF]Mattei et al., 2014a,b), some others have done decomposition of the problem into a set of subproblems [START_REF] Karunanithi | A new decomposition-based computer-aided molecular/mixture design methodology for the design of optimal solvents and solvent mixtures[END_REF], while some have solved the problem globally for a given application, for example polymer blends [START_REF] Vaidyanathan | Computer-aided synthesis of polymers and blends with target properties[END_REF]. As part of a methodology for the design of formulated products, Gani and co-workers (Conte and Gani, 2011;Conte et al., 2011;Mattei et al., 2014a,b) have conceived the Virtual Product-Process Design Laboratory. They propose to run sequentially a design scenario within a computer aided stage: select a problem template and translate product needs into properties (Mattei et al., 2014a,b), choose an active ingredient of the product from the database, then design the solvents with their MIXD algorithm either from a pre-defined list or generated with a CAMD tool [START_REF] Conte | Innovation in Integrated Chemical product-process Design: Development through a Model-based Systems Approach[END_REF] and then add additives from another list and finally end up with the optimization of composition. To escape the computer-aided stage, a verification scenario is run with more accurate models, possibly involving model developments. An ultimate experimental validation ends the design activity. For overcoming the problem of consumer attributes not described by models, Solvason et al. (2009) combined an enumerating CAMD technique and MDOE (mixture design of experiments) technique. Illustrated with the formulation of a refrigerant mixture, they first solve a reverse formulation problem to find property relations that match user-defined attributes. Those relations are then used as target of a reverse problem aiming at finding the suitable mixture.

Methods and data structures

We have developed a CAPD tool, named as IBSS (Integrated Bio Sourced Search). It follows the general methodology of CAMD tools and is aimed at finding mixtures in which some molecule may bear bio-sourced fragments. The problem of finding a single molecule is handled as a mixture with one element. The methods and data structures developed to cope with that tailoring are now presented.

Optimization problem

The CAPD problem is multi-objective since several properties must be matched. It is transformed into a single-objective problem, aiming at maximizing a global performance, GloPerf, described by an objective function OF, subject to k equality constraints and i inequality constraints on property targets P. It can be modelled as follows:

OF = max(GloPerf (MG i , z i , cond j )) s.t. P k (MG i , z i , cond j ) = P k,fixed P l,lowerbound ≤ P l (MG i , z i , cond j ) ≤ P l,upperbound s.t. constraints on MG i , z i , cond j

(1)

The optimization variables are the molecular graph structure MG i of the individual i mixture components, the mixture composition z i and j conditions cond j . The conditions, cond j , affect the performance calculation by imposing conditions under which the properties are calculated.

The optimization variables can be constrained to allow the user to tailor the problem: the composition of any molecule, z i and condition cond j can be fixed, bounded or free. For example, the user can impose mole fraction of an ingredient, specify a physical state of the molecule or define the range of operating conditions. Any molecule MG i of the mixture can be fixed (ex. an active ingredient), sourced from a list of molecule (ex. a list of additives or solvents) or left free for optimization. In that latter case, one or more chemical fragments can be fixed or taken from a list of fragments to design the molecule (ex. to impose a renewable material derivative fragment in the molecule).

The global performance, GloPerf, is formulated as the product of a penalty function and of a weighted sum of np individual performance PropPerf p with weight w p with respect to each property target.

GloPerf (MG i , z i , cond j ) = min ır =1 (ı r • (1 -Penal r )) • np p=1 w p • PropPerf p (MG i , z i , cond j ) np p=1 w p (2)
The penalty function min (ı r •(1 -Penal r )) is related to user defined rules. Each rule r contains data related to a molecular pattern described as an opened molecular graph and is assigned a penalty percentage Penal r . ı r is equal to 1 if the rth rule is violated, 0 otherwise. Typical rules describe unrealistic structures from the chemical synthesis point of view, or molecular patterns that are correlated with toxicity.

Each individual performance PropPerf p for the property p, compares the predicted value x with the targeted value P. The user can select among mathematical functions F(x) as shown in Table 1: Gaussian [START_REF] Venkatasubramanian | Computer-aided molecular design using genetic algorithms[END_REF], desirability functions [START_REF] Castillo | Modified desirability functions for multiple response optimization[END_REF] or straight functions.

CAMD solution robustness is shattered by the property model prediction uncertainty. Solutions have been proposed in the literature, like the use of fuzzy logic operators to define upper and lower bounded property ranges associated to degrees of satisfaction [START_REF] Ng | Optimal chemical product design via fuzzy optimisation based inverse design techniques[END_REF], as can be done here with the straight function representation. Alternatively, the knowledge of property model uncertainty for some group contribution methods [START_REF] Hukkerikar | Group-contribution+ (GC+) based estimation of properties of pure components: improved property estimation and uncertainty analysis[END_REF] can be used to define the Tol parameter in the Gaussian function representation.

The search algorithm

The search algorithm selected is a genetic algorithm with elitism policy as earlier proposed by [START_REF] Venkatasubramanian | Computer-aided molecular design using genetic algorithms[END_REF] in CAMD. Modification operators are added to alter the mixture composition, conditions and molecules and to perform a multilevel search. The population size, the elitism value, the number of level and all the probabilities of operators are defined by the user.

The initial population of individuals is generated randomly within the predefined constraints on the optimization variables related to MG i , z i , cond j . The method for building fragments from chemical building blocks is described later.

The CAPD search can be performed in several sequential levels [START_REF] Harper | Computer-aided molecular design with combined molecular modeling and group contribution[END_REF][START_REF] Korichi | Computer aided aroma design I -molecular knowledge framework[END_REF]. At low level, simple and/or fast-computing property prediction models are used over a large population. Then as the level increments, more complex and/or time-consuming models are used over a smaller population originated from the fittest individuals of the previous level population. At the next level, the same set of building blocks and molecular structures is kept. In the meantime, the objective function can be modified according to the user's initial choices: property estimation models can be dropped, added or substituted by more complex ones.

Mixture representation data

Mixture representation

The mixture structure is customisable as presented in Fig. 1. Each mixture is an assembly of items and conditions. Each item contains one molecule and one mole fraction value. Each molecule is further split into interconnected fragments. The fragments are further built from basic or complex functional groups.

Initially, the user defines the mixture structure: the number of molecules, their type (fixed, list or free) and composition constraints. For each free molecule, he sets the number of fragments, fragment type (fixed, list or free) and fragment interconnections. For each free fragments, he defines the building groups to be used and their maximum number Different building block list can be used for different fragments. A molecule may contain a single free fragment. In that case the fragment has no external connections.

Molecular representation

We have selected molecular graph for the molecular representation which is described by an adjacency matrix (Achenie et al., 2003). The diagonal elements are either a hydrogen-suppressed atom basic group, or a complex group or a fragment. The offdiagonal elements are bond type connections. A matrix with a diagonal that contains basic groups exclusively is called an extended molecular graph hereafter.

A molecule graph is the aggregation of its fragment graphs. Fig. 2 describes the acetoin molecule (3-hydroxybutanone) built as a structure of three interconnected fragments. Fragment 1 is built from two basic groups and one complex group. Fragments 2 and 3 are just one basic group.

Basic groups (BGs) represent a hydrogen-suppressed atom, like O, OH, CH 3 , CH 2 , >CH 2 , NH 3 , N, etc. BGs are often similar to some first order groups in group contribution methods. They are assigned a "BG.ID" attribute. It is an 'elementary group' integer EG = P 1 P 2 P 3 P 4 that is displayed in the extended molecular graph diagonal. P 1 refers to the atomic number preceded with a 1 (106 for C, 107 for N, 108 for O, 117 for Cl. . .), P 2 refers to the highest bond order, P 3 to the type of the atom attached to, including its occurrence in an aromatic or non-aromatic cycle, and P 4 to the number of implicit hydrogen atoms bonded to the atom. BGs also bear a "bondvector" attribute that represents the number of single, double and triple bonds. e.g. for BG.short formula = "=C<", BG.ID = "106200" and BG.bondvector = "[2;1;0]".

Complex functional groups (CGs) are multi-atom groups. They are useful for a compact description of multi-atom chemical functions, like carboxyl groups, R COOH, nitrite R O N O, nitro R NO 2 , peroxy ROOR ′ , ester RCOOR ′ , acetals RCH(OR ′ )(OR ′′ ), sulfenyl RSOR ′ . . . A non-exhaustive list of BGs and CGs along with BG.ID is provided as Supplementary material. As they are kept intact when applying modification operators, complex groups are suitable to describe bio-sourced molecule derivative or synthons and to keep them in the molecule candidates.

CGs inherit from the basic group attributes, but the "ID" attribute has no specific meaning and is assigned a unique incremental value 2xxxxx defined by the user. Additional CGs attributes are a molecular graph "CG.graph" describing the complex group in terms of basic groups and a "CG.connectionVector" integer attribute that represents the number and location of the external connections of the group in the molecular graph (see Fig. 2).

Fragments "Fgt" are described as an adjacency matrix and an external connection vector. The adjacency matrix contains basic or complex functional group information (see Fig. 2). The vector represents the single, double or triple bond external connections, and their location on the functional groups. Similar bond type external connections on one group are distinguished by a letter.

Group vector

The classification of functional groups by [START_REF] Constantinou | Computer aided product design: problem formulations, methodology and applications[END_REF] is adapted into group vectors "GV" to provide the list of basic and complex functional groups authorized by the user to build a free fragment where k groups are allowed. A group vector GV is represented by the following way:

GV = {N 1 , N 2 , . . ., N n } (3)
Where N i is the number of groups in the fragment that have i connections, from 1 to n. Basic functional groups have up to 4 connections. Some groups with sulphur and phosphorous atoms where atom valence can be 6 and 5 respectively, though they are described with 4 connections (see a list of basic groups in Supplementary material). For complex groups, the number of connections n can be higher than 4. Authors have used some groups with 6 external connections, especially those that are sourced from renewable synthons.

In addition, the following chemical feasibility rules coming from the octet rule hold:

n j=1 N j = k (4) n j=1 N j (2 -j) = 2m -extconnections (5)
where m is a number equal to 1 minus the maximum number of cycles allowed by the user and extconnections is the number of external connections of the fragment.

Fragment creation

The method for building a free fragment from a preselected list of basic or complex functional groups is developed to ensure a diversity of structures, in particular heterocycles and aromatic cycles which may often disappear during modification by a genetic algorithm.

The procedure for constructing a fragment is described as an activity diagram in Fig. 3 which is explained below:

(1) Initialization of the user parameters: minimum k min and maximum k max number of BGs and/or CGs functional groups in a fragment, the maximum number of cycles m max, the number of external connections for the fragment extconnections and the preselected BGs and CGs ordered into N j lists.

(2) Choice of the fragment parameters: a value of k between k min and k max and a value of m between 1 and 1 -m max are randomly selected. Using those values, the possible group vectors are determined, e.g. assuming that groups with up to 4 external connections exist (n = 4) and k = 6, five group vectors are possible GV 1 = {1,5,0,0}, GV 2 = {2,3,1,0}, GV 3 = {3,1,2,0}, GV 4 = {3,2,0,1}, GV 5 = {4,0,1,1}. A GV i is picked and used as a basis for the fragment construction. Each group vector has the same probability to be chosen to ensure a higher diversity of the generated structures.

(3) Addition of a new element to the fragment: At this point either an acyclic group or an entire cycle structure can be added. The decision is made randomly considering the number of cycles yet to be constructed and the remaining number of elements that can be inserted. If an acyclic group is selected, a N j is selected randomly from GV i . Then one of the groups of the N j list is picked and is added.

Step 3 is repeated until all k groups have been added (go to step 6). If a cycle is selected go to step 4. (4) Cycle building: The cycle size and its aromaticity are decided before the cycle construction starts. Then all the elements that form the cycle are inserted one after the other until the cycle is complete. Side branches to the cycle are inserted only once the cycle is closed to the cycle groups that still bear unsaturated external connections. This way the elements of a cycle are consecutive in the adjacency matrix to make the graph easier to handle. (5) Fused cycle building: At the end of step 4 it is possible to add a fused cycle. The decision is made randomly considering the number of cycles yet to be constructed and the remaining number of elements that can be inserted. Then size and aromaticity are decided. For a non-aromatic fused cycle, the attachment points of the fused cycle are searched on the last inserted cycle and its adjacent cycles. A couple of attachment points are randomly chosen and the shortest way between these two points is determined. This path is then considered as a part of the fused cycle to build. Then the number of elements still to be added is randomly chosen and the elements are added one by one. The process to construct an aromatic cycle is the same with supplementary constraints: the attachment points must be consecutive in the last cycle; they must be connected with a double bound and must concern aromatic groups. This process goes on until the fragment is complete (go to step 6). (6) Fragment complete: All the groups having been added, the corresponding molecular consistency is tested. If so, the possible complex groups are expanded into basic groups to create the expanded molecular graph of the fragment.

Molecule modification operators

The operators used to alter the molecular graph adjacency matrix are summarized in Fig. 4: mutation, crossover, insertion, deletion and substitution.

Regarding the insertion and deletion operators, we have added a specific branch constructor or destructor. We have also developed a new substitution operator to improve the chances of structural changes in aromatic rings without breaking their aromaticity.

• The mutation operator consists in the random replacement of a single group by a group from the N i reference list that bears the same external connections, e.g. >NH by >CH 2 both with two single bonds in Fig. 4. • The crossover operator consists in randomly choosing two identical non-cyclic bond types (single, double or triple bound) in two extended molecular graphs. This is checked with the P 3 value in the ID attribute of a group. The semi-graphs are then switched and recombined to form two new molecules. • The insertion operator consists in the random addition of a group in the graph. We have developed the possibility to insert a group that has more than two connections, like CH<, thus enabling the completion of the graph with a new branch (Fig. 4). • The deletion operator consists in the random removal of a group with at least two connections of the same type in the graph. To be consistent with the insertion operator, it is possible to delete a group that has more than two external connections. This may possibly induce the deletion of a side branch of the molecule. In Fig. 4, the group C< is randomly chosen for deletion. The extra branches are deleted and the two remaining branches are directly reconnected. The branch NH is deleted because it is not consistent with the remaining connection type. Thus another suitable group F is reconnected instead.

• The new substitution operator combines the principles of mutation and insertion and consists in the replacement of a group by a group that has more connections. It was developed because the other operators performed poorly in modifying aromatic cycles without destroying their aromaticity: a mutation operator alone would be ineffective because of the limited number of aromatic groups, here only enabling to replace C arom H by N arom and conversely. The crossover operator cannot be applied on rings.

The insertion and deletion operators would only allow adding or deleting the aromatic heteroatoms O and NH , to maintain the aromaticity.

Implementation

The methods described above have been implemented in the "IBSS" CAPD tool software prototype. The iterative process of the IBM-RUP software development method [START_REF] Kroll | The rational unified process made easy: a practitioner's guide to the RUP[END_REF] was applied. It is centred on the architecture and driven by the functional needs that should cover the CAPD tool. Those needs were expressed by the partners of the French ANR CDP2D 2009 project InBioSynSolv, aiming at designing biosolvents. Those needs highlighted that the sustainability of the candidate mixture may rise from the occurrence of bio-sourced fragments within the molecular structure and by the use of EHS property models to evaluate the performance of each candidate. To speedup the implementation process, Model Driven Engineering, MDE principles were followed with the help of UML 2.0 (Unified Modelling Language) and BPMN (Business Process Modelling Notation) diagrams. It produced architectural, behavioural, functional and structural UML 2.0 views that are now briefly presented.

Architectural view

The CAPD tool is developed as component-based software. Each of the three components, 'Man-Machine-Interface MMI'; 'SEARCH' and 'P3', is a software package that encapsulates a set of functions and data and communicates through interface with the other components (Fig. 5). Next to them, an XML-structured database contains the basic and complex groups.

In the spirit of MDE, effective coding was started after the CAPD tool architecture and needs definition were well advanced. The first MMI component is written in java and aims at providing an input XML file to the search component. The second component "SEARCH" is built around an object-oriented architecture and is written in C# within the Visual Studio ® environment. The third property calculation component "P3" is a Dynamic-Link Library written in VB.NET. It contains a library of property estimation models and automatic group finders or molecular descriptor routines used to translate the molecular graph information sent by the search component into suitable inputs for the property estimation models. As a standalone component, various interface methods enable the use of the P3 component with the search component or with other independent software. Currently, thirty properties can be estimated from twenty property estimation models which are listed in the Supplementary material.

Behavioural view

The behavioural view presents the different processes of the tool and highlights the components interoperability (Fig. 6).

The interoperability between the MMI and the Search components is asynchronous via an XML file as input and a text file for the results. The interoperability between the Search and the Property calculation components is synchronous and Windows-Library like. The XML file is generated through the MMI-component and loaded in the search component interface. Then, the Resolution package launches the search algorithm, a genetic algorithm with a single level in Fig. 6. First the initial population is generated. Then properties are evaluated by the property calculation component to calculate the performance of each mixture. The population is modified and evaluated again until a stop criterion is satisfied putting an end to the search. Results are then saved in a text format to the MMI package of the search component. The data in these files are afterwards displayed by the MMI component.

Functional view

The functional view highlights the potential users and the main functionalities of the software, like the function "launching a CAPD search" represented in Fig. 7. In compliance with the hierarchical decision making process for sustainable product design presented elsewhere [START_REF] Heintz | Chemical enterprise model and decisionmaking framework for sustainable chemical product design[END_REF], we distinguish the expert user from the basic user. The former can access all parameters. The role of basic user is intended for people with moderate expertise in property estimation and chemistry, typically business or technical managers. The basic user has access to following functionality of the tool: to select pools of chemical building blocks classified by raw material sources, to choose product requirements to be considered, to define the rough structure of the product mixture and to set generic product constraints, like a fixed ingredient within the mixture. Martin and Young ( 2001) Log(BCF) [START_REF] Veith | Structure-activity correlations in studies of toxicity and bioconcentration with aquatic organisms[END_REF] The mixture data should preferably be set before the objective function data because some property calculation models must be chosen for each mixture component which number must be known beforehand. The algorithm data can be specified at any time. Any data set can be saved and retrieved independently. When the data are correctly defined, the application offers the user to run the search algorithm. The results are displayed as a list of product candidate which the user can choose to save.

Structural view

The structural view concerns the model abstractions, object classes and their relationships.

• The mixture data structure implements the hierarchy described earlier in Figs. 1 and2: mixture > composition + molecules + conditions > fragments > building blocks; along with the methods like the fragment constructor displayed in Fig. 3. A routine is developed to use the basic and complex groups stored in the "Block" database.

• The search algorithm data structure contains the genetic algorithm parameters and the modification operator methods which are described earlier.

• The objective function structure contains all the data and methods related to property target values and performance function in addition to the property estimation model for pure compounds and mixtures. In addition, we allow the basic user to be logged in the MMI component to identify product requirements as needs [START_REF] Yunus | A systematic methodology for design of tailor-made blended products[END_REF] rather than as precise property names and models. Requirements are edited by an expert user into a set of calculable properties, with specific targets and property estimation models as done by others (Mattei et al., 2014a,b). Table 2 displays an example.

Case studies

5.1. Case study 1: blanket wash mixture substitution 5.1.1. Problem setting "Blanket wash" are needed to clean ink residues from rubber blankets in the lithographic printing process. In replacement of petroleum sourced solvents, [START_REF] Sinha | A systematic method for integrating product attributes and molecular synthesis[END_REF] used a CAPD approach on a pre-selection of seven water-soluble and low EHSimpacts solvents and solved a MINLP problem to find the optimal composition of the aqueous blend. [START_REF] Heintz | Chemical enterprise model and decisionmaking framework for sustainable chemical product design[END_REF] revisited the whole decision process leading to substitute blanket wash. It led to the specification of a tree of requirements translated here into property target values for a CAPD search with the IBSS tool. We search for an aqueous mixture and we optimize simultaneously the organic molecule and its molar composition.

A first requirement consists in solubilizing the phenolic resin ink which is evaluated by the Relative Energy Difference, RED, property given by Eq. ( 6). RED is computed by using Hansen solubility If RED < 1, the solvent dissolves the solute. For the best solubilization results, we set a target value of RED = 0.

Other requirements refer to the cleaning process conditions: the solvent should keep its fluidity over the rubber blanket surface, which is evaluated by setting target values for viscosity and surface tension. It should comply with Volatile Organic Compound, VOC limits, evaluated with the vapour pressure. Environmental Health and Safety issues are assessed with a five EHS indices model and a flammability class is evaluated with a flash point model. The property target values are displayed in Table 3 along with the property weights, operation conditions, linear or non-linear mixture models and pure compound model. A Gaussian function is used for each property performance in Eq. ( 2).

In the global performance Eq. ( 2), we set penalties, namely the occurrence of three consecutive oxygen atoms or of C-C-C ring (Penal r = 100% => forbidden).

The organic solvent structure is split into two fragments, one with a core synthon traceable from various renewable biomass material stocks and the other built from less than 10 groups among a pool of simple and complex groups displayed in Fig. 8.

The search is ran with the following set of parameters: number of generations (300), population size (100), elitism (30) and the probabilities of crossover, mutation, insertion and deletion (respectively 65, 15, 10 and 10). It is completed in less than 40 min. Cyclic compounds are excluded. The probability for composition change and molecule change are 0.7 and 0.3 respectively.

Results

The output file displays a list of hundred mixtures rated by their performance [START_REF] Heintz | An integrated framework for product formulation by computer aided mixture design[END_REF] and the best solution (Fig. 9) emerged after 40 generations. [START_REF] Sinha | A systematic method for integrating product attributes and molecular synthesis[END_REF] suggested a mixture of ␥butyrolactone and water (45/55 mol%) for which posterior evaluation of the performance with our criteria is 0.94. Our solution reaches a performance of 0.96 for a composition of 31 mol%. Confidential issues prevent us to display the formula of organic molecule which is different than ␥-butyrolactone. Fig. 9 shows how the organic molar fraction affects the RED property. The minimum is at 0.3, indicating that the genetic algorithm correctly finds the best solution. It also shows that adding water to the biosolvent improves its cleaning ability. a G(Tol,Val): Gaussian type function see Table 1.

Case study 2: substitution of chlorinated paraffins

Problem setting

We seek an alternative molecule to chlorinated paraffins (CPs). With a general formula C n H 2n+2-x Cl x , they have been used as additives in high-temperature lubricants and cutting fluids for metalworking, plasticizers and flame retardants in plastics, sealants and leather [START_REF] Bayen | Chlorinated paraffins: a review of analysis and environmental occurrence[END_REF]. Sourced from petroleum, their use and risks (toxicity, potential carcinogenic) have been assessed and are regulated (EU Report, 2008).

Table 4 summarizes the targeted physical properties. We seek a suitable alternative sourced from levulinic acid (LA), which is known as a biomass platform chemical with derivatives already in use as lubricants, coatings and printing/inks [START_REF] Bozell | Production of levulinic acid and use as a platform chemical for derived products[END_REF]. The structure of substitute molecule consists in two fragments, one is fixed as LA and the other is constructed as assembly of less than 10 groups among 19 basic groups and 8 complex groups as displayed in Fig. 10.

Two different searches are executed: one uses only basic and complex groups to generate candidates (set 1); the other uses a fixed fragment (levulinic acid) and generate LA derivatives using basic and complex groups (set 2).

The following parameters are used: number of generations (300), population size (100), elitism (10) and the probabilities of crossover, mutation, insertion and deletion (respectively 20, 50, 15 and 15). The formation of cyclic and aromatic compounds is not allowed. No penalty related to specific chemical sub-structure is set.

Results

The search with set 1 (no LA fragment) produces the molecule shown in Fig. 11 that was brought out by the genetic algorithm after the 100th generation with a performance of 0.9122. It has great similarities with known CPs alternatives and displayed in Fig. 11. Notice that viscosity is not computed because of no suitable group Introducing LA fragment in the molecule, we display the influence of k max (maximum number of building groups in fragment 2), on the performance of the genetic algorithm in Fig. 12. It shows that acceptable solutions are found if we use k max > 4. The maximum performance (1.000) is always found for k max ≥ 6, and found more rapidly with k max = 7.

With the set 2 (LA + basic + complex groups), the maximum performance (1.00) is achieved after 100 generations for k max = 7. The best 10 molecules of the 100th generation are given in Fig. 13, along with a general scheme for the present approach and their predicted values of various properties in Table 5.

For the set 2, the flash point of molecules 6 and 9 could not be estimated. Indeed, group contributions do not exist for the group >CO (molecules 1-8 and 10) and for the group CHCO (molecule 9). The Configuration Interaction "CI" correction of the Hukkerikar model was not used [START_REF] Hukkerikar | Group-contribution+ (GC+) based estimation of properties of pure components: improved property estimation and uncertainty analysis[END_REF]. For the molecules other than 6 and 9, the viscosity was not predicted because the Joback and Reid method does not handle the phosphate group. Besides, the method of Hukkerikar et al. does not have the contribution of the group aC-PO4 for the boiling point; therefore, the values of boiling point shown in Table 5 are only first approximations that do not take into account the aC-PO4 group. Such property model deficiencies for complex structures confirm that CAPD computer based approaches should necessarily be validated by laboratory experiments. As mentioned before, the global performance was renormalized.

Analysing the results for molecules 6 and 9, predicted melting points values are greater than the target values but remains below the reported average absolute error of the model (17.65 K). These molecules can be retained for further experimental analysis of their melting point. Nevertheless, they will likely be waived by chemists because they exhibit a double ketone sequence, C( O) C( O) , that is known to be unstable. Alternatively, this rule could be coded in the IBSS tool and a penalized performance could have been used instead.

Case study 3: find solvents for extraction of natural antioxidants

Problem setting

The objective is to find a solvent for methyl p-coumarate (MpCA), an ester of cinnamic acid isolated from plants. This class of compounds is known for their antibacterial, antifungal and/or antiviral activity [START_REF] Galanakis | A knowledge base for the recovery of natural phenols with different solvents[END_REF][START_REF] Sova | Antioxidant and antimicrobial activities of cinnamic acid derivatives[END_REF]. [START_REF] Panteli | Solubilities of cinnamic acid esters in organic solvents[END_REF] found that tert-pentanol was the best solvent for MpCA among tert-butanol, tert-pentanol, ethyl acetate and n-hexane.

Here we search for a glycerol based solvent. The molecular structure we look for consists in two fragments. Fragment Fgt1 is taken from a list of two glycerol derivatives with either sn-1 or sn-2 substituent. The other fragment is constructed as assembly of less than 10 groups among those displayed in Fig. 14. We use the IBSS tool with a two level search. At level 2, we keep the three best molecules from the list of candidates found at level 1 and improve the prediction of their solubility by using a solid-liquid equilibrium calculation instead of the RED property (see Eq. ( 6)) used at level 1. The solubility x of an ideal solid phase in a liquid solvent is given by Eq. ( 7) in first approximation [START_REF] Prausnitz | Molecular thermodynamics of fluid phase equilibria[END_REF]:

ln x = H m RT 1 - T T m -ln (7)
where is the activity coefficient (computed by the UNIFAC modified Dortmund 1993 model), H m and T m are the fusion enthalpy and the melting temperature of MpCA, respectively. Table 6 summarizes the physical properties to be satisfied by a potential candidate.

The following set of parameters are used: number of generations (300), population size (50), elitism (5), k max (4) and the probabilities of crossover, mutation, insertion and deletion (respectively 40, 50, 5, 5). The formation of cyclic and aromatic compounds is not allowed.

Results

The maximum performance (0.9802) is achieved in 60 generations. The best 10 molecules of the 60th generation are given in Fig. 15 and their properties are displayed in Table 7.

None of those molecules exhibit a sn-2 glycerol derivative. Indeed, the contribution to the melting point predictions of the sn-2 OCHCH2OH group is much larger (10.9421) than the sn-1 OCH2CHOH group contribution (1.7527) in the Hukkerikar's method [START_REF] Hukkerikar | Group-contribution+ (GC+) based estimation of properties of pure components: improved property estimation and uncertainty analysis[END_REF]. Thus, isomers can show a difference in melting point up to 100 K. Even though all the generated Fig. 13. Candidates for chlorinated paraffins substitution. Building blocks: levulinic acid + basic groups + complex groups. 100th generation, kmax = 7. molecules have predicted melting points greater than the target, this difference is within the reported average absolute error of the model (17.65 K), and these molecules can be retained for further experimental analysis. It can be also noted that the best molecule is a compromise, with a larger RED than other molecules but with an overall better performance.

The top 3 molecules proposed by level 1 are retained for further analysis of their solubility with the UNIFAC-SLE approach. The results of level 2 shift the order of molecules 1 and 2 in terms of performance compared to that of the level one by using the RED model: the molar fraction solubility of molecule 2 is predicted at 0.0787, lower than that of molecule 1 at 0.0896. The predicted solubility of molecule 3 is the highest: 0.0902. This value confirms with the level 1 RED prediction, in which molecule 3 is that with the lowest RED value.

Conclusion

We have described the methods, data structures and software implementation of IBSS, a computer aided product design (CAPD) tool that is aimed at finding mixtures. Based on CAMD concepts, it is able to optimize simultaneously the mixture components, compositions and mixture conditions.

With the help of a model driven engineering approach, the architecture and the functional needs of the CAPD tool were devised, and specifically aimed at finding mixtures with molecules that can be sourced from pools of renewable materials. This is done by setting constraints on the mixture molecules or on fragments within molecules, like bio-sourced synthons.

A molecular representation based on an adjacency matrix was selected. It was made flexible to represent both atom-based structures and fragment-based structures. Diagonal hydrogensuppressed atomic elements of the matrix were described with a novel coding of four indexes describing the atom number, the highest bond type, the atom neighbouring context and the number of hydrogen atoms. The coding was also devised for the identification of chemical groups or descriptors necessary to evaluate properties with models from various authors. A complex group describing polyatomic structures was created to represent polyatomic chemical functions and synthons. It enabled to keep intact bio-sourced synthons during the search and promote their occurrence in the final solution if their performance was deemed high enough.

A genetic algorithm was selected to optimize simultaneously the mixture elements, its composition and additional operating conditions. It was made capable to perform a multilevel search with different objective functions. Modification operators were adapted to cope with the mixture search context and with the aforementioned molecular representation. Classical crossover and mutation operators were used, while insertion and deletion operators were adapted to insert or delete side branch. A new substitution operator was proposed to maintain cyclic molecules through the genetic algorithm generations. The building of fragments from basic or complex groups was made more efficient by using a vector group classification inventorying building groups on the basis of their external connections.

The mixture performance is calculated through a sum of weighted property performance that can be penalized if specific molecular patterns occur, like those found in toxic molecules.

The tool was implemented as a set of three independent software components, namely a MMI component, a CAPD search component and a property library component; associated to a database of basic and complex functional groups. To cope with the difficulty for some users of expressing requirements in terms of numerical target values, we distinguished product requirement (better suitable to basic user) vs calculable properties (for expert users).

Future work is ongoing to set the CAPD tool within a virtual laboratory decision-making framework as described by [START_REF] Heintz | Chemical enterprise model and decisionmaking framework for sustainable chemical product design[END_REF]. We also intend to benefit from the flexible tool architecture, first by adding alternative solving methods other than genetic algorithm, and by testing the suitability of the molecular graph representation for novel property models based on higher dimensionality which are useful to handle conformation dependent properties.

Fig. 1 .

 1 Fig. 1. Overview of the mixture structure and its substructures.

Fig. 2 .

 2 Fig.2. Molecular graph representation of 3-hydroxybutanone as three fragments.

Fig. 3 .

 3 Fig. 3. Activity diagram of the creation of a free fragment object.

Fig. 4 .

 4 Fig. 4. Molecule modification operators.

Fig. 5 .

 5 Fig. 5. CAPD tool UML component diagram.

Fig. 6 .

 6 Fig. 6. BPMN diagram of the three IBSS components behaviour.

Fig. 7 •

 7 Fig. 7 describes the use case diagram of the 'launch a search' functionality. It gives access to the definition of the problem data through three data subsets and describes actions available to the user.

Fig. 7 .

 7 Fig. 7. "Launch a search" UML use case diagram.

Fig. 8 .

 8 Fig.8. Organic molecular structure and building blocks in the blanket wash mixture search.

Fig. 9 .

 9 Fig. 9. Performance and property values for the best biomass derivative mixture and influence of its fraction on RED property.

Fig. 10 .

 10 Fig. 10. A chlorinated paraffin example and building blocks used to generate alternative molecules.

Fig. 11 .

 11 Fig. 11. Candidate and known alternatives for chlorinated paraffins (CP) substitution.

Fig. 12 .

 12 Fig. 12. Influence of the maximum number of building blocks on the genetic algorithm performance.

Fig. 14 .

 14 Fig. 14. Methyl p-coumarate solubilization: target molecule and building blocks.

Fig. 15 .

 15 Fig. 15. Candidates for methyl p-coumarate solubilization. Glycerol and basic groups as building blocks; 60th generation.

Table 1

 1 Property performance functions.

Table 2

 2 Examples of product requirements and associated calculable properties.

	Product requirement	Calculable	Default pure compound
		property	calculation model
	Fluidity	Viscosity Molecular	Conte et al. (2008) Atomic weight summation
		weight
	Volatility	Boiling point Vapor pressure	Marrero and Gani (2001) Riedel (1954)
		Log(Kow)	Marrero and Gani (2002)
	Toxicity	-Log(LC50)

Table 3

 3 Blanket wash substitution calculable properties, target, model and parameters.

	Property name	Weight	Target	Performance function a	Mixture model	Pure cpnd model	Operating conditions
	Molecular weight	1	<200 g/mol	G(20,0.8)	-
	Flash point	1	>323.15 K	G(5,0.8)	Linear	Catoire et al. (2006)
	Vapor pressure	1	<0.00267 bar	G(10 -4 ,0.8)	Linear	Riedel (1954)	T = 298.15 K
	RED	4	=0	G(1,0.9)	Volumic fraction	HSPiP	See Eq. (6)
	Env. waste	0.2	>8	G(1,0.8)	Linear	Weis and Visco (2010)
	Env. impact	0.2	>8	G(1,0.8)	Linear	Weis and Visco (2010)
	Health	0.2	>8	G(1,0.8)	Linear	Weis and Visco (2010)
	Safety	0.2	>8	G(1,0.8)	Linear	Weis and Visco (2010)
	LCA	0.2	>8	G(1,0.8)	Linear	Weis and Visco (2010)
	Viscosity	1	∈ [0.8;1.4] cP	G(0.1,0.8)	Non linear; Tamura	Conte et al. (2008)	T = 298.15 K
					and Kurata (1952)
	Surface tension	1	∈ [30;45] dyn/cm 2	G(5,0.8)	Non linear; Rice	Conte et al. (2008)	T = 298.15 K
					and Teja (1982)
	Density	1	∈ [0.9;1.1]	G(0.05,0.8)	Linear	HSPiP 3.1, Model, 2010
	Log(Ws)	4	>4 mg/L	G(0.5,0.8)	Linear	Marrero and Gani (2002)

Table 4

 4 CPs substitution case: product requirements and calculable properties with corresponding CAMD parameters.

	Product requirement	Calculable property	Weight	Target	Performance function a	Pure cpnd model	Operating conditions
	Liquid at usage temperature	Melting point Boiling point	1 1	<283.15 K >523.15 K	G(10,0.8) G(5,0.8)	Hukkerikar et al. (2012) Hukkerikar et al. (2012)
	Non flammable at high temperatures	Flash point	1	>505.15 K	G(5,0.8)	Hukkerikar et al. (2012)
	Fluidity	Viscosity	1	>37 cP	G(10,0.8)	Joback and Reid (1987)	T = 298.15 K
	Low potential to bioaccumulation	Log(Kow)	1	<3	G(2,0.8)	Hukkerikar et al. (2012)
	a G(Tol,Val): Gaussian type function see Table 1.

Table 5

 5 Predicted properties for the best 10 candidates for chlorinated paraffins substitution (100th generation, kmax = 7, levulinic + basic groups + complex groups).

		Performance	Melting point (K)	Boiling point (K)	Flash point (K)	Viscosity (cP)	Log(Kow)
	Molecule 1	1.0000	274.80	358.67	518.95	-	0.31
	Molecule 2	1.0000	266.34	360.24	507.30	-	0.30
	Molecule 3	0.9999	283.30	362.01	520.23	-	0.36
	Molecule 4	0.9986	284.72	361.60	506.98	-	-0.97
	Molecule 5	0.9982	284.94	362.59	527.74	-	-0.82
	Molecule 6	0.9976	284.14	317.97	-	35.19	-0.72
	Molecule 7	0.9964	285.71	353.64	512.73	-	-1.45
	Molecule 8	0.9902	287.39	369.86	537.69	-	0.37
	Molecule 9	0.9884	280.08	307.51	-	32.38	-0.96
	Molecule 10	0.9656	291.30	364.25	522.11	-	-1.33

Table 6

 6 MpCA solvent extraction case: product requirements and calculable properties with corresponding CAMD parameters. The Hukkerikar model is used to compute the Hansen solubility parameters, which are used to compute RED.

	Product requirement	Calculable property	Weight	Target	Performance	Pure cpnd model	Pure cpnd model
					function a	Level 1	Level 2
	Must dissolve MpCA	RED	2	=0	G(1,0.95)	Hukkerikar et al. (2012) b	UNIFAC + SLE at 25 • C
	Liquid at room temperature	Melting point Boiling point	1 1	<283.15 >373.15 K	G(10,0.8) G(5,0.8)	Hukkerikar et al. (2012) Hukkerikar et al. (2012)	--
	Non flammability	Flash point	1	>343.15 K	G(1050.8)	Hukkerikar et al. (2012)	-
	Low potential to bioaccumulation	Log(Kow)	1	<3	G(2,0.8)	Hukkerikar et al. (2012)	-

a G(Tol,Val): Gaussian type function see Table

1

.

b

Table 7

 7 Predicted properties for the best 10 candidates for methyl p-coumarate solubilization (60th generation, kmax = 4, glycerol + basic groups).

	Level 1

Acknowledgment

This scientific work was supported by the French National Research Agency (InBioSynSolv ANR-CP2D-2009-08) in partnership with Rhodia-Solvay Company, ENSCL and LCA-INPT.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compchemeng. 2014.09.009.