Emilie Bosc

Romuald Pépion

Patrick Le Callet

Martin Köppel

Patrick Ndjiki-Nya

Luce Morin

Muriel Pressigout

Perceived quality of DIBR-based synthesized views

Keywords: Virtual view synthesis, multi-view video, 3DTV, quality assessment, quality metrics

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Any new device or method needs to be assessed through adapted and reliable tools. 3D Video technology has brought out new challenges such as the question of DIBR-based views evaluation. 3D Video applications, such as 3D Television (3DTV) or Free viewpoint Television (FTV) provide the users a novel type of experience thanks to the added value (depth, or immersion), compared to conventional 2D video. However, the success of these applications depends on their ability to offer high-quality visual content. FTV and 3DTV require virtual view generation, which is ensured by the use of DIBR algorithms. This process induces new types of artifacts and consequently impacts on quality that has to be identified considering various contexts of use. Artifacts from 2D video compression have been widely studied, but DIBR deals with new types of distortions that are mainly geometric distortions.

Yet, synthesized views, generated either from decoded and distorted data or from original data, need to be assessed. The best assessment tool remains the human judgment as long as the right protocol is used. Subjective quality assessment is still delicate while addressing new types of conditions because one has to define the optimal way to get reliable data. Subjective tests are time-consuming consequently one should draw big lines on how to conduct such experiment to save time and observers. Since DIBR is introducing new conditions, the right protocols to assess the visual quality with observers are still an open question. Protocols might vary according to the purpose (impact of compression, DIBR techniques comparison …).

Objective metrics are requested because subjective quality assessment tests are cumbersome. They are meant to predict human judgment and to facilitate the quality evaluation. However, the evaluation of their performances relies on their closeness to subjective assessment results. As, the way to conduct the subjective quality assessment protocols is already questionable, reliable objective quality metrics among existing ones that could be useful in DIBR context are all the more uncertain.

Yet, trustworthy working groups partly base their future specifications, concerning new strategies for 3D video, on the outcome of objective metrics. Considering the test conditions may rely on usual subjective and objective protocols (because of their availability), the outcome of wrong choices could result to a poor quality of experience for users, or worse, the failure of 3D Video for industrials. This study is motivated by two main questions: are usual requirements for subjective evaluation protocols still appropriate for assessing 3D synthesized views? And second, are usual metrics sufficient for assessing 3D synthesized views?

Previous studies ( 1 , 2 ) already proposed to answer our second questioning, regarding the reliability of usual objective metrics. However, often, objective metrics are applied on single distorted views and these measurements are compared to subjective assessments on 3D displays whose synthesis process is blind. In 2 , artifacts from both synthesis and compression are assessed without distinction.

This study aims at answering the two questions above, but considering the complexity of 3D Video, and in order to precisely keep in check what is being measured, a first approach is to be independent from the 3D, that is to say both the stereopsis and the 3D display whose technology is still a major factor of quality degradation. Thus, the study presented in this chapter focuses on the quality evaluation of DIBR-based synthesized views in 2D conditions. Besides, these conditions are plausible in a Free Viewpoint Video (FVV) application.

The problem we consider has not been treated in such a way in the literature. In our context, the results of our study confirm the common idea that 3D video is not just the extension of 2D video and that commonly used 2D video metrics are not sufficient to assess the synthesized views quality. The synthesized images quality is measured through pixelbased or perceptual-like metrics. The images are synthesized with seven different DIBR algorithms. Correlation with human perception is estimated thanks to subjective assessments.

The document is organized as follows: Section 2 highlights the challenges brought by 3D in the context of quality assessment. Section 3 presents the material used for the experiments. Section 4 considers the first question regarding the subjective assessment protocols while Section 5 studies the reliability of the objective tools. Section 6 concludes the paper.

NEW CHALLENGES IN TERMS OF QUALITY ASSESSMENT

Recently, the complexity of 3D video has been considered. It is now admitted that 3D Video is not the trivial extension of 2D to 3D, thanks to studies such as Seuntiens' 3 . In the latter, Seuntiens claims that not only the image quality should be addressed in a quality evaluation, but as well new factors such as naturalness, visual experience, and presence. This suggests the use of a multidimensional measurement that does not exist yet as an objective tool, and that has not been deployed yet in subjective quality assessment of 3D media. Indeed, it is a tough task designing subjective quality assessment frameworks, when new concepts are introduced. Besides, Chen et al. [START_REF] Chen | New requirements of subjective video quality assessment methodologies for 3DTV[END_REF] proposed new requirements from this perspective. They proposed to reconsider several conditions in this context, such as the viewing conditions (viewing distance, monitor resolution), the test material (depth rendering according to the chosen 3D display), viewing duration, etc. This certainly reinforces our assumption concerning the fact that 3D Video is much more complex than 2D conventional video, and that tools usually thought of as appropriate in the context of 2D conventional video should be reconsidered.

As a matter of fact, the study of artifacts occurring when generating a novel view from DIBR algorithms reveals that the distortions are mainly located around strong discontinuities, transitions of depth. They are quite different from artifacts we are used to deal with in 2D video compression. This comes from the fact they arise from the very synthesis process, and are then geometric distortions.

These distortions are directly related to the quality of the incoming data values (e.g. quality of depth estimation) and the synthesis process strategies. The latter mainly aimed at dealing with the critical problem in DIBR, namely the disocclusion. This term refers to the discovered regions during the synthesis of a new viewpoint. Extrapolation techniques are meant to fill the disoccluded regions. In this section, we comment the artifacts from seven different synthesis process strategies. These seven DIBR algorithms are referenced from A1 to A7:

-A1: based on Fehn [START_REF] Fehn | Depth-image-based rendering (DIBR), compression and transmission for a new approach on 3D-TV[END_REF] , where the depth map is pre-processed by a low-pass filter. Borders are cropped, and then an interpolation is processed to reach the original size.

-A2: based on Fehn [START_REF] Fehn | Depth-image-based rendering (DIBR), compression and transmission for a new approach on 3D-TV[END_REF] . Borders are inpainted by the method proposed by Telea 6 .

-A3: Tanimoto et al. [START_REF] Tanimoto | Reference Softwares for Depth Estimationand View Synthesis[END_REF] , it is the recently adopted reference software for the experiments in the 3D Video group of MPEG.

-A4: Müller et al. [START_REF] Müller | View synthesis for advanced 3D video systems[END_REF] , proposed a hole filling method aided by depth information.

-A5: Ndjiki-Nya et al. [START_REF] Ndjiki-Nya | Depth image based rendering with advanced texture synthesis[END_REF] , the hole filling method is a patch-based texture synthesis.

-A6: Köppel et al. [START_REF] Köppel | Temporally consistent handling of disocclusions with texture synthesis for depth-image-based rendering[END_REF] , uses depth temporal information to improve the synthesis in the disoccluded areas.

-A7: corresponds to the unfilled sequences (i.e. with holes).

Commonly encountered artifacts, in DIBR-based synthesized images are described in the following. In most of the cases, these artifacts are located around strong discontinuities, but are more perceptible when the contrast between the background and the foreground is important.

Object shifting: a region may be slightly translated or resized, depending on the chosen extrapolation method (if the method chooses to assign the background values to the missing areas, object may be resized), or on the encoding method (blocking artifacts in depth data results in object shifting in syntheses).

Blurry edges: when data is pre-processed by low pass filters, it leads to blurry synthesized regions. Indeed, to avoid, the lack of information in disoccluded areas, it is tempting to low pass filter the per-pixel depth information, as proposed in A1. However, this smoothes out the discontinuities and results in blurry regions in synthesis. This is annoying especially around the edges that are expected to be sharp. It is obvious around the background/foreground transitions. These remarks are confirmed on Figure 1-A1, where we also note the curved stairway, and on Figure 1-A3 around the disoccluded areas.

Texture synthesis: inpainting methods can fail in filling complex textured areas. A2 and A3 use the inpainting method in [START_REF] Telea | An Image Inpainting Technique Based on the Fast Marching Method[END_REF] . However, this inpainting method fails to reconstruct complex textures. A4 also fails in rendering complex textures. To overcome these limitations, a hole filling approach based on patch-based texture synthesis is proposed in A5. However, candidate patches may not be perceptually close to the disoccluded region and then lead to artifacts. Furthermore, the rectangular patches can lead to blocky artifacts.

Flickering: when errors occur randomly in depth data along the sequence, pixels are wrongly projected: some pixels suffer slight changes of depth, which appears as flickers in the resulting synthesized pixels. To avoid this effect, proposed the use of a background sprite is proposed in A6. Background image information is stored in the sprite and updating by the original and the synthesized image information of previous frames to provide depth consistency. All along the frames, background texture is copied from the stored information to the disoccluded region. The copy relies on depth data, and is then still dependent on depth data quality Tiny distortions: in synthesized sequences, a large number of tiny geometric distortions and illumination differences are temporally constant and perceptually invisible. However, pixel-based metrics may penalize these distorted zones.

Figure 1: DIBR-based synthesized frames of the "Lovebird1" sequence.

EXPERIMENTAL MATERIAL

3D Video involves many factors that are not completely under control yet. Among them, the stereopsis phenomenon, or the relationship between depth rendering and the quality of experience over 3D-adapted displays, can be cited. To be independent of combination of uncontrolled parameters, the quality of images synthesized from DIBR algorithms is first assessed in 2D conventional conditions. This is meant to be a preliminary experiment, and it is also a plausible context for an FVV application.

A1 Original depth Original image

Warped depth map A7 A5 A3 A4

A6

Three different MVD sequences are used. The sequences are Book Arrival (1024x768, 16 cameras with 6.5cm spacing), Lovebird1 (1024x768, 12 cameras with 3.5 cm spacing) and Newspaper (1024x768, 9 cameras with 5 cm spacing). These sequences are also involved in the test benchmarks of 3DV Group of MPEG.

Seven DIBR algorithms processed the three MVD sequences to generate, for each one, four different viewpoints. Six of the DIBR algorithms, from A1 to A6, are cited in section 2, and the seventh correspond to the unfilled sequences (i.e. with holes) is labelled A7.

Subjective and objective evaluations were processed to evaluate the synthesized views quality.

As for the subjective quality evaluations, standardized methodologies for subjective multimedia quality assessment, such as Absolute Categorical Rating (ACR) [START_REF]ITU-T P.910 Subjective video quality assessment methods for multimedia applications[END_REF] , have proved their efficiency regarding the quality evaluation of 2D conventional images. ACR was then used to collect perceived quality scores. The tests were conducted in an ITU conforming test environment. The stimuli were displayed on a TVLogic LVM401W, and according to ITU-T BT.500 [START_REF] Bt | Methodology for the subjective assessment of the quality of television pictures[END_REF] .

As for the objective quality metrics, measurements were realized over the 84 synthesized views by the means of MetriX MuX Visual Quality Assessment Package [START_REF]MetriX MuX Home Page[END_REF] except for two metrics: VQM and VSSIM. VQM were available at [START_REF]Video Quality Research[END_REF] ; VSSIM was implemented by the authors, according to [START_REF] Wang | Video quality assessment based on structural distortion measurement[END_REF] . The reference was the original acquired video. Still image metrics were applied on each frames of the sequences and then averaged by the number of frames.

SUITABILITY OF SUBJECTIVE QUALITY ASSESSMENT PROTOCOLS

ACR is selected as an example of subjective assessment protocol to reconsider. ACR methodology consists in presenting test objects (i.e. images or sequences), one at a time and in a random order, to the observers. Observers score the test item according to a discrete category rating scale. This quality scale is depicted in Table 1: ACR quality scale. A semantic meaning is assigned to the five integers of the quality scale. The results of an ACR test are obtained by averaging observers' opinion scores for each stimulus, in other words, by computing mean opinion scores (MOS). ACR requires many observers to minimize the contextual effects (previously presented stimuli influence the observer opinion, i.e. presentation order influences opinion ratings). The accuracy increases with the number of participants. The efficiency of ACR methodology is proved as soon as at least 24 observers participate, according to VQEG.

ACR test was conducted over the 84 video sequences synthesized from the seven DIBR algorithms.

Table 2 shows the ranking of the algorithms from the obtained subjective scores. The obtained values indicate the quality of the presented sequences is perceived as "poor", in average. Although ACR scores allow the ranking of the algorithms, they do not directly provide knowledge on the statistical equivalence of the results. Student's t-test has been performed over the scores to meet this need.

Table 1 depicts the results of the Student's t-test processed with the values. The number in parentheses indicates the required minimum number of observers that allow the statistical distinction (VQEG recommends 24 participants as a minimum [START_REF]VQEG 3DTV Test Plan for Crosstalk Influences on user Quality of Experience[END_REF] , values in bold are higher than 24 in the table). It seems that in some cases, 24 participants are not sufficient to provide a sharp distinction. This may be due to the fact that some of the algorithms induce the same type of distortions. Then, the perceived quality is similar to observers.

Table 3: Results of Student's t-test with ACR results. Legend: ↑: superior, ↓: inferior, O: statistically equivalent. Reading: Line"1" is statistically superior to column "2". Distinction is stable when "less than 7 observers participate.

A1 A2 A3 A4 A5 A6 A7 A1 ↑(7) ↑(3) ↑(3) ↑(2) ↑(3) ↑(1) A2 ↓(7) ↑(2) ↑(2) ↑(1) ↑(2) ↑(1) A3 ↓(3) ↓(2) O (>32) ↑(9) O (>32) ↑(1) A4 ↓(3) ↓(2) O (>32) O (>32) O (>32) ↑(1) A5 ↓ (2)
↓( 1) ↓( 9)

O (>32) ↓(15) ↑(1) A6 ↓(3) ↓(2) O (>32) O (>32) ↑(15) ↑(1) A7 ↓(1) ↓(1) ↓(1) ↓(1) ↓(1) ↓(1)

RELIABILITY OF OBJECTIVE METRICS

Different objective metrics are used in the presented experiments:

-Peak Signal to Noise Ratio (PSNR) measures the signal fidelity of a distorted image compared to a reference. It is based on the measure of the Mean Squared Error (MSE).

-As an alternative to pixel-based methods, Universal Quality Index UQI [START_REF] Wang | A universal image quality index[END_REF] is a perceptual-like metric. It models the image distortion by a combination of three factors: loss of correlation, luminance distortion, and contrast distortion.

-PSNR-HVS [START_REF] Egiazarian | New full-reference quality metrics based on HVS[END_REF] , based on PSNR and UQI, is meant to take into account the Human Visual System (HVS) properties.

-PSNR-HVSM [START_REF] Ponomarenko | On between-coefficient contrast masking of DCT basis functions[END_REF] is based on PSNR but takes into account Contrast Sensitivity Function (CSF) and betweencoefficient contrast masking of DCT basis functions.

-Single-scale Structural SIMilarity (SSIM) [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] combines image structural information: mean, variance, covariance of pixels, for a single local patch. Multi-scale SSIM (MSSIM) is the average SSIM scores of all patches of the image.

-Video Structural Similarity Measure (V-SSIM) [START_REF] Wang | Video quality assessment based on structural distortion measurement[END_REF] is a FR video quality metric which uses structural distortion as an estimate of perceived visual distortion. At the patch level, SSIM score is a weighted function of SSIM of the different component of the image (i.e. luminance, and chromas). At the frame level, SSIM score is a weighted function of patches' SSIM scores (based on the darkness of the patch). Finally at the sequence level, VSSIM score is a weighted function of frames' SSIM scores (based on the motion).

-Visual Signal to Noise Ratio (VSNR) [START_REF] Chandler | VSNR: A wavelet-based visual signal-to-noise ratio for natural images[END_REF] is also a perceptual-like metric: it is based on a visual detection of distortion criterion, helped by CSF.

-Weighted Signal to Noise Ratio (WSNR) that uses a weighting function adapted to HVS denotes a weighted Signal to Noise Ratio as applied in [START_REF] Damera-Venkata | Image quality assessment based on a degradation model[END_REF] and expressed in [START_REF] Lee | Foveated video quality assessment[END_REF] .

-Visual Information Fidelity (VIF) [START_REF] Sheikh | Image information and visual quality[END_REF] uses a statistical model and an image distortion model to assess an image. VIFP is a pixel-based version of VIF.

-Information Fidelity Criterion (IFC) [START_REF] Sheikh | An information fidelity criterion for image quality assessment using natural scene statistics[END_REF] uses a distortion model to evaluate the information shared between the reference image and the degraded image.

-Noise quality measure (NQM) [START_REF] Damera-Venkata | Image quality assessment based on a degradation model[END_REF] quantifies the injected noise in the tested image.

-Video Quality Metric (VQM) was proposed by Pinson and Wolf in [START_REF] Webster | An objective video quality assessment system based on human perception[END_REF] . VQM combines different perceptual calculated parameters (motion, blurring, global noise, color distortion) to compute the final quality measure.

Table 4 shows the rankings obtained from the objective scores (and the subjective scores as a reminder). Subjectivebased ranking points A1 as the best out of the seven algorithms. But objective-based rankings judge it as the worst. As well, A6 obtains the best results for numerous objective metrics, but subjective assessments judge the synthesized views as "poor". From the observation of this table, our assumption is that objective metrics detect non-annoying artifacts. Note, that objective metrics are coherent to one another. They are correlated.

The consistency between objective and subjective measures is evaluated by computing the correlation coefficients for the whole fitted measured points.

Table 5 shows the results. None of the tested metric is 50% close to human judgment when assessing synthesized views.

According to MOS and PC scores, the most correlated metrics are VSNR, VIFP and IFC (47.3%, 46.9%, and 45.6%).

In order to consider further the consistency, Figure 2 shows the ACR values over the VSNR scores. The figure shows that a difference of nearly 3dB (A3 and A6) leads to nearly the same ACR score. However the range of confidence intervals is quite important for cases such as A3's. This can be explained by the fact that the results presented in Figure 2 are obtained from the twelve synthesized sequences of each algorithm whose quality scores varied depending on the baseline distance between reference cameras and target viewpoints, and also depending on the sequences. These results show inconsistencies between subjective and objective scores. This suggests that objective metrics are not able to reliably predict human judgment, in the tested conditions.

CONCLUSION

This paper addresses the issue of evaluating virtual synthesized views with tools usually appropriate in the context of 2D conventional video quality assessment. The assessments of the seven test algorithms by objective measurements and subjective ratings show that among all tested objective metrics, VSNR VIFP and IFC are the most correlated with perceptual evaluation provided by MOS scores, with low correlation scores though. As well, the results also show metric's inability to predict human experience, in the tested conditions. The objective metrics certainly indicated the presence of distortions, but they were not able to predict annoying distortions in the sense of human opinion. New methods are then required for assessing virtual synthesized views as pixel-based and perceptual-based metrics fail.

Because view synthesis induces geometric distortions, a new model should take into account depth data and particularly the preservation of coherent edges after the synthesis.

Figure 2 :

 2 Figure 2: Correlation between MOS and VSNR.

Table 2 :

 2 Ranking of algorithms according to subjective scores.

		A1	A2	A3	A4	A5	A6	A7
	ACR	2.70	2.41	2.14	2.03	1.96	2.13	1.28
	Rank order	1	2	3	5	6	4	7

Table 4 :

 4 Rankings according to measurements.

		A1	A2	A3	A4	A5	A6	A7
	ACR	2.70	2.41	2.14	2.03	1.96	2.13	1.28
	Rank order	1	2	3	5	6	4	7
	PSNR	19.02 24.99 23.227 25.994 26.035 26.04 20.89
	Rank order	7	4	5	3	2	1	6
	SSIM	0.648 0.844 0.786 0.859 0.859 0.859 0.824
	Rank order	7	4	6	1	1	1	5
	MSSIM	0.664 0.932 0.825 0.948 0.948 0.948 0.888
	Rank order	7	4	6	1	1	1	5
	VSNR	13.14 20.41 18.75 21.786 21.965 21.968 20.73
	Rank order	7	5	6	3	2	1	4
	VIF	0.129 0.393 0.313 0.423 0.423 0.424 0.396
	Rank order	7	5	6	2	2	1	4
	VIFP	0.153 0.415 0.342 0.446 0.446 0.446 0.419
	Rank order	7	5	6	1	1	1	4
	UQI	0.359 0.664	0.58	0.598 0.598 0.598 0.667
	Rank order	7	5	6	3	3	3	1
	IFC	0.779 2.399 1.926 2.562 2.562 2.564 2.404
	Rank order	7	5	6	2	2	1	4
	NQM	8.66 15.933 13.415 16.635 16.739 16.739 10.63
	Rank order	7	4	5	3	1	1	6
	WSNR	14.41 20.85 18.853 21.76 21.839 21.844 16.46
	Rank order	7	4	5	3	2	1	6
	PSNR HSVM 13.99 19.37 18.361 21.278 21.318 21.326 16.23
	Rank order	7	4	5	3	2	1	6
	PSNR HSV 13.74 19.52 17.958 20.795 20.823 20.833 15.91
	Rank order	7	4	5	3	2	1	6
	VSSIM	0,662 0,879 0,809 0,899 0,898 0,893 0,854
	Rank order	7	4	6	1	2	3	5
	VQM	0.888 0.623 0.581 0.572 0.556 0.557 0.652
	Rank order	7	5	4	3	1	2	6

Table 5 :

 5 Correlation coefficient between objective and subjective scores.

	PSNR SSIM MSSIM VSNR	VIF VIFP UQI	IFC NQM WSNR PSNR	PSNR	VSSI	VQM
									HVSM	HVS	M
	ACR 34.5	45.2	27	47.3	43.9 46.9	20.2 45.6 36.6	32.9	34.5	33.9	33	33.6

ACKNOWLEDGEMENTS

We would like to acknowledge the Interactive Visual Media Group of Microsoft Research for providing the Breakdancers data set, the MPEG Korea Forum for providing the Lovebird1 data set, the GIST for providing the Newspaper data set, and HHI for providing Book Arrival. This work is supported by the French National Research Agency as part of PERSEE project (ANR-09-BLAN-0170) and CAIMAN project (ANR-08-VERS-002).


Jean-Charles Manceau
email: jc.manceau@brgm.fr

Jeremy Rohmer

Jérémy Rohmer

Analytical solution incorporating historydependent processes for quick assessment of capillary trapping during CO 2 geological storage

Keywords: CO 2 geological storage, safety, quick assessment, residual trapping, hysteresis

Analytical solution incorporating history-dependent processes for quick assessment of capillary trapping during CO2 geological storage

Introduction

In addition to non-emitting and renewable energy production among others, carbon dioxide capture and storage (CCS) is considered as a promising option to ensure the necessary decrease of CO 2 anthropogenic emissions (IPCC 2005).

Regarding the storage capacity needed to develop geological sequestration at a large scale, saline aquifers offer good potentials [START_REF] Bachu | Sequestration of CO 2 in geological media in response to climate change: road map for site selection using the transform of the geological space into CO 2 phase space[END_REF]. Its industrial development in those formations is above all conditioned by its safety demonstration (IEAGHG 2007). A comprehensive strategy should then be adopted, made up of a risk assessment, a monitoring plan and a corrective measures plan as stated by the European Directive on the geological storage of carbon dioxide (EC 2009).

Several modes of trapping take part in the confinement within saline aquifers: as a gas phase, the buoyant CO 2 can be trapped below an impermeable cap rock layer (structural trapping, e.g. [START_REF] Bachu | Aquifer disposal ofCO 2 : hydrodynamic and mineral trapping[END_REF] or disconnected from the main cluster and immobilized due to capillary forces (capillary or residual trapping, e.g. [START_REF] Juanes | Impact of relative permeability hysteresis on geological CO 2 storage[END_REF]; dissolved in brine, it can be trapped in the aqueous phase (solubility trapping, e.g. Ennis-King and Paterson 2005) before reacting with rocks and precipitating (mineral trapping, e.g. [START_REF] Bachu | Aquifer disposal ofCO 2 : hydrodynamic and mineral trapping[END_REF][START_REF] Gunter | Aquifer disposal of CO 2 -rich greenhouse gases: Extension of the time scale of experiment for CO 2 -sequestering reactions by geochemical modeling[END_REF]). These trapping mechanisms occur under specific conditions (notably geological, hydraulic and chemical) and over different time scales. Structural trapping alone is not sufficient to prevent migration of CO 2 in case of cap-rock anomaly for instance [START_REF] Nghiem | Risk mitigation trough the optimization of residual gas and solubility trapping for CO 2 Storage in saline aquifers[END_REF]. The evaluation of safer modes (residual, solubility and mineral trapping) is then essential; this study is focused on residual trapping as it has been shown as the fastest mode of trapping (e.g. [START_REF] Qi | Design of carbon dioxide storage in aquifers[END_REF]Juanes et al. 2010).

When injecting supercritical CO 2 (gas-like since it is less dense and viscous than the native brine -we use the term gas in the following) into saline aquifers (liquid phase), the native brine is drained by CO 2 (the liquid saturation decreases). The formation minerals being naturally water-wet, the liquid phase can be considered as the wetting fluid (Juanes et al. 2005). When the injection stops, according to the driving forces, some parts of the initial plume will experience a wetting stage (gas saturation decreases) while others will continue being in drainage stage.

Wetting (also called imbibition) may either occur naturally in response to regional groundwater flow or to gravity currents (i.e. "natural" imbibition), or artificially in response to liquid injection through wells in depth (i.e. "artificial" imbibition).

Artificial imbibition aims at increasing the storage efficiency (e.g. [START_REF] Qi | Design of carbon dioxide storage in aquifers[END_REF] or at forcing the immobilization of the injected CO 2 in case of significant "irregularities" detected at the reservoir level [START_REF] Manceau | Forcing gaseous CO 2 trapping as a corrective technique in the case of abnormal behavior of a deep saline aquifer storage[END_REF] or at the level of a potentially impacted groundwater aquifer [START_REF] Esposito | Remediation of Possible Leakage from Geologic CO 2 Storage Reservoirs into Groundwater Aquifers[END_REF].

When imbibition occurs, the gas saturation decreases until a residual gas saturation that corresponds to the amount of gas that is capillary trapped and hence immobilized. This Darcy's scale observation is explained (e.g. by [START_REF] Lenormand | Mechanisms of the displacement of one fluid by another in a network of capillary ducts[END_REF][START_REF] Suicmez | Effects of wettability and pore-level displacement on hydrocarbon trapping[END_REF]) by the displacement phenomena governing the interactions between the wetting and non-wetting phase at pore scale (i.e. piston-like displacement, cooperative pore-body filling and snap-off).

These phenomena lead to the immobilization of the wetting liquid during drainage (through by-passing) and of the non-wetting gas (through by-passing and snapoff) explaining the residual trapping observed. The trapped fraction depends on the initial saturation since, for instance, when the non-wetting phase initial concentration is lower, the number of potential sites where snap-off can occur is lower as well.

Moreover, a contact angle hysteresis modification between drainage and imbibition occurs due to the surface roughness and due to the wettability change once the non-wetting phase has been introduced [START_REF] Valvatne | Predictive pore-scale modeling of multiphase flow[END_REF]de Gennes et al. 2004). This hysteresis effect is taken into account at Darcy's scale with the dependence of the characteristic curves describing the two-phase flow (i.e. evolution of the relative permeability and capillary pressure functions) on the saturation history. In order to measure supercritical CO 2 -brine characteristic curves, several experimental technics exist as the steady and unsteady state usual approaches. However, the specificities of the CO 2 (notably the reactivity with rock) and of the reservoir rock (types and heterogeneity) as well as the large range of pressure-temperature conditions that have to be tested make difficult the full interpretation of the measurements [START_REF] Müller | Supercritical CO 2 -brine relative permeability experiments in reservoir rocks -Literature review and recommendations[END_REF]. Hence, the hysteresis phenomenon as well as its amplitude still remains little known in the field of geological CO 2 storage.

Aiming either at evaluating the plume behaviour under a natural groundwater flow or designing the injection parameters (e.g. rate, injection period) in order to create an "artificial imbibition", an estimation of the temporal trapped gas quantity evolution and of the spatial distribution of gas saturation (mobile and immobile) is of first importance.

In the present paper, we develop in this purpose an analytical solution applying the classical fractional flow theory, and including trapping and hysteresis models to incorporate history dependant processes. In this study, we focus on the "artificial" imbibition case.

Review of existing models

By its very objective, the developed solution should describe the evolution of the gas plume as well as the residual trapping that occurs as soon as brine is injected.

Both phenomena are described through different kinds of models from which our solution is derived. In this matter, we present a brief overview of the existing models.

Trapping and hysteresis models

Multi-phase flow modelling implies the use of characteristic curves describing, at Darcy's scale, the interactions between the different phases. They represent the evolutions of both capillary pressure and relative permeabilities (associated to the wetting and the non-wetting phase). In a first approach, the capillary pressure and relative permeabilities are functions of one phase saturation only (e.g. [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF][START_REF] Brooks | Hydraulic properties of porous media[END_REF].

However these models fail in representing the hysteretic phenomena explained in the previous section. In this view, the trapped fraction should be determined at first (with a trapping model) and, subsequently, the history-dependent relative permeabilities and capillary pressure functions should be assessed as a function of the saturation, of the trapped fraction and of the process in progress (with a hysteresis model). Please note that the trapped wetting phase fraction is commonly considered constant. As an illustration of the hysteresis model, Figure 1 depicts the evolution of the gas relative permeability during drainage and wetting processes.

One of the most classic non-wetting phase trapping models is Land's model (1968). This model being only suitable for water-wet media, [START_REF] Spiteri | A new model of trapping and relative permeability hysteresis for all wettability characteristics[END_REF] developed a new model for other wettablilities. Several hysteresis models have been built based on Land's trapping model, amongst whom are [START_REF] Carlson | Simulation of relative permeability hysteresis to the nonwetting phase[END_REF], and Parker and Lenhard's [START_REF] Parker | A model for hysteretic constitutive relations governing multiphase flow -1. Saturation-pressure relations[END_REF][START_REF] Lenhard | A model for hysteretic constitutive relations governing multiphase flow -2. Permeability-saturation relations[END_REF].

They are based on classical non-hysteretic relative permeability and capillary pressure models and give the characteristic curves during drainage and wetting (for the primary drainage-imbibition cycle and for the following cycles).

[Fig. 1]

Models describing gas plume dynamics

Safety (e.g. leakage concerns) and performance (e.g. injectivity issues) assessments in the field of CO 2 geological storage notably require plume extension models. Many analytical and numerical models have been developed, based on existing fluid mechanics theories. Analytical models allow quick and flexible computations while numerical models allow dealing with the complexity of the phenomena.

Two main families have emerged among analytical and semi-analytical models aiming at characterizing the evolution of one phase in another one [START_REF] Dentz | Abrupt-Interface Solution for Carbon Dioxide Injection into Porous Media[END_REF]. The first approach studies the shape of one phase plume assuming a marked delimitation between the two phases that do not mixed (at Darcy's scale sense) while the second is focused on the evolution of one phase in terms of its saturation profile.

The first family of models (called here sharp interface models) are based on the so-called sharp interface approximation stipulating that the CO 2 saturation gradient thickness is low compared to the distances concerned. Capillary forces are neglected in these models. After Bear's work (1972), some expressions have been derived specifically for CO 2 sequestration in deep saline aquifers, first for the injection period when gravity forces can be neglected [START_REF] Nordbotten | Injection and storage of CO 2 in deep saline aquifers: analytical solution for CO 2 plume evolution during injection[END_REF][START_REF] Nordbotten | Similarity solutions for fluid injection into confined aquifers[END_REF] and then for post-injection period (with gravity forces) [START_REF] Hesse | Gravity currents in horizontal porous layers: transition from early to late self-similarity[END_REF]. While these works were carried out for horizontal and undisturbed aquifers, others studied the effects of sloppy formations [START_REF] Hesse | Scaling analysis of the migration of CO 2 in saline aquifers[END_REF]) and of groundwater flow [START_REF] Juanes | Upscaling of capillary trapping under gravity override: application to CO 2 sequestration in aquifers[END_REF][START_REF] Juanes | The Footprint of the CO 2 Plume during Carbon Dioxide Storage in Saline Aquifers: Storage Efficiency for Capillary Trapping at the Basin Scale[END_REF][START_REF] Macminn | CO 2 migration in saline aquifers. Part 1: Capillary trapping under slope and groundwater flow[END_REF]. Some of the previously cited works [START_REF] Hesse | Scaling analysis of the migration of CO 2 in saline aquifers[END_REF][START_REF] Juanes | Upscaling of capillary trapping under gravity override: application to CO 2 sequestration in aquifers[END_REF][START_REF] Juanes | The Footprint of the CO 2 Plume during Carbon Dioxide Storage in Saline Aquifers: Storage Efficiency for Capillary Trapping at the Basin Scale[END_REF][START_REF] Macminn | CO 2 migration in saline aquifers. Part 1: Capillary trapping under slope and groundwater flow[END_REF]) also included capillary trapping after imbibition following [START_REF] Kochina | Groundwater mound damping[END_REF].

However, these models assumed a constant residual gas saturation after wetting, hence neglecting the history-dependent processes.

The second family of models (named here Buckley-Leverett models) uses the fractional flow theory and is based on Buckley-Leverett equation [START_REF] Buckley | Mechanism of fluid displacements in sands[END_REF]) that, neglecting gravity and capillary pressure, allows finding a solution for saturation variability. Concerning CO 2 geological sequestration, [START_REF] Noh | Implications of coupling fractional flow and geochemistry for 2 injection in aquifers[END_REF] and [START_REF] Zeidouni | Analytical solution to evaluate salt precipitation during CO 2 injection in saline aquifers[END_REF] used this approach to model drainage occurring during injection. They took into account solubility and partitioning of the several components (CO 2 , water) into both phases; [START_REF] Noh | Implications of coupling fractional flow and geochemistry for 2 injection in aquifers[END_REF] also considered capillary trapping during wetting with a constant residual gas saturation.

These two kinds of models are very similar since the underlying theory is alike.

The governing equation comes from continuity and Darcy's equations. In the sharp interface models, the purpose is to estimate the shape of the gas plume. The sharp interface and the vertical equilibrium (or Dupuit) approximations are made so that the gas saturation can be replaced by the height of gas (interface height).

The capillary pressure, which is neglected in the Buckley-Leverett models, is replaced, in the sharp interface models, by a "pseudo-capillary" pressure [START_REF] Gasda | Vertical equilibrium with sub-scale analytical methods for geological CO 2 sequestration[END_REF], which is the difference of both phases fluid-static pressures. Practically, the additional approximations made in the sharp-interface models allow using the Buckley-Leverett models in two dimensions [START_REF] Yortsos | A theoretical analysis of vertical flow equilibrium[END_REF].

In this paper, we propose a solution for the gas plume evolution taking into account variable residual gas saturation. As described by the existing trapping models, this value is highly influenced by the gas saturation reached before imbibition. Moreover, these history-dependent effects have been shown to be significant (this was numerically demonstrated for instance by [START_REF] Doughty | Modeling geologic storage of carbon dioxide: Comparison of non-hysteretic and hysteretic characteristic curves[END_REF].

Therefore, in our study, we prefer Buckley-Leverett models to the sharp interface ones since the saturation variability is of first importance. A similar choice has been made by [START_REF] Furati | Effects of relative permeability history dependence on two-phase flow in porous media[END_REF], followed by van [START_REF] Van Kats | A mathematical model for hysteretic two-phase flow in porous media[END_REF], who studied the effects of relative permeabilities hysteresis.

Mathematical model

The model should enable the representation of a CO 2 injection (at a constant rate ܳ ) followed by a forced imbibition (at a constant rate ܳ ) and the associated saturation variations in terms of trapped and mobile gas (see the conceptual layout on Figure 2). We focus on the artificial brine injection and use a 1D radially symmetric geometry.

[Fig. 2]

We apply the fractional flow theory [START_REF] Buckley | Mechanism of fluid displacements in sands[END_REF], considering the flow incompressible and immiscible.

The continuity equation gives:

(1)

డௌ డ௧ ொ ೕ గఝ డ డ మ ൌ 0,
where ݂ ൌ is the fraction of gas in the total flowing stream. Note that ݍ ௧௧ ൌ ܳ , ܳ being the volume injection rate (constant over the CO 2 and brine injection periods, but possibly different).

Developing ݂ with Darcy's equation and neglecting gravity (1D horizontal model) and capillary pressure, we obtain:

(2)

݂ ൌ ା ൌ ଵ ଵା ഊ ഊ
with the mobility ߣ ൌ ఓ . ݂ is then dependent on the saturation only (and on the history), which allows writing:

(3)

డௌ డ௧ ொ ೕ గఝ ௗ ௗௌ డௌ డ మ ൌ 0.
We scale equation (3) for the gas injection period and note ܶ the injection time at ܳ rate and ݎ the radius of the cylinder with a height ݄ and porosity ߮ containing the volume of gas injected during ܶ and we choose:

߬ ൌ ௧ ் ೕ and ߦ ൌ బ .
The governing equation for the gas injection period gets simpler:

(4)

డௌ డఛ ௗ ௗௌ డௌ డక మ ൌ 0
For the brine injection period, i.e. ݐ ܶ , we decide:

(5)

߬ ൌ 1 ொ ொ ௧ି் ೕ ் ೕ
This arbitrary choice leads to the same governing equation for both injection periods. The reformulation of the equation ( 4) leads to the so-called Buckley-Leverett equation:

(6) ቀ డక మ డఛ ቁ ௌ ൌ ௗ ௗௌ Classically, ߞ ൌ క మ
ఛ denotes the similarity variable, i.e. the suitable combination of primary variables corresponding to the velocity of the waves, which are solutions of the conservation law (see e.g. [START_REF] Lax | The formation and decay of shock waves[END_REF]). In other words, the equation ( 6) states that the increasing rate of the "dimensionless area" ߦ ଶ (in this case of an axisymmetric geometry) swept over by a specified saturation is given by

ௗ ௗௌ , i.e.
by the variations of the gas fractional flow when saturation changes.

Then, the Buckley-Leverett equation naturally implies studying fractional flow function. As mentioned previously, the fractional flow is a function of the dynamic viscosities, which can be taken constant under the incompressible, immiscible and isothermal flow assumptions and of the relative permeabilities that depends on the saturation as well as on the history of the saturation evolutions.

Because of hysteresis effects, separated works are done according to the flow process occurring (drainage or imbibition).

Phase I: gas injection

When injecting gas in native brine, both liquid and gas permeabilities will follow the so-called primary drainage curve (that depends on the model chosen).

The flow function is generally considered as S-shaped in literature (e.g. Furati [START_REF] Buckley | Mechanism of fluid displacements in sands[END_REF]. This multivalued situation is due to the no-capillary pressure approximation. The discontinuity must satisfy the jump condition [START_REF] Lax | The formation and decay of shock waves[END_REF][START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF]). This gives the velocity of the shock (denoted :)ܫ݄݇ܿݏ

Phase II: Brine injection

When injecting brine, the gas saturation decreases and one moves along the fractional flow curve in the opposite direction than for phase I. However, this function is different from the drainage one, since the relative permeability curves are prone to hysteresis effects.

The problem is the same though, since the curve is still S-shape meaning that a saturation discontinuity will occur at the leading edge of the injected brine, with a velocity given by the Rankine-Hugoniot condition. Contrary to the previous phase, gas is still present and flows upward as well as downward of the imbibition front. Moreover, the saturation on both sides of this front evolves over time: . This expression is clearly complex, as infinity of such imbibition gas fractional flow curves have to be computed to solve the entire problem.

We therefore decide to introduce the following approximation ܵ ,௦ூூ ି ൌ ܵ ሺܵ ,௦ூூ ା ሻ; we call it the "instantaneous trapping" approximation since we consider that, as soon as the imbibition front reaches a specific saturation, a jump towards the associated residual saturation is instantaneously done. Please note that this assumption is conservative in that the front velocity obtained is lower.

This choice leads to: This assumption is justified by the following reasoning. Gas fractional flow function depends on mobility ratio (ߣ ߣ ⁄ ሻ in that the lower this ratio, the higher the gas flow and thus the more left-shifted the inflection point of the S-shaped function. Given relative permeability functions, one can infer that the saturation interval where the gas fractional flow is convex is smaller when the viscosity ratio decreases (e.g. [START_REF] Mcworther | Exact integral solutions for two-phase flow[END_REF]). In the case of CO 2 geological storage, gas viscosity is much lower than liquid one and the viscosity ratio is particularly low.

The difference between the saturation directly downstream the imbibition front and the residual gas saturation is then likely to be slight (additional discussion on the validity of this approximation is provided in section 5.2).

By introducing this assumption, it is possible to establish for each brine injection time the saturation profile that is made of the profile corresponding to the drainage phase downward of the advancing brine front and, upward, of the residual saturation calculated from the saturation reached by the brine advancing front at earlier time (see Figure 3 for the illustration of the elaboration methodology of the gas saturation profile). Taking ܵ * ൌ ܵ ,௦ூ ି , one retrieves the complete trapping dimensionless time.

One can deduce the dimensionless radius at which this collision occurs:

(16) ߦ ௦,ௌ * ൌ ߦ ௌ * ሺ߬ * ሻ ൌ ට ௗ ,ೌೌ ௗௌ ሺܵ * ሻ • ߬ሺܵ * ሻ,
from which the dimensioned volume of mobile gas ܸ as well as the mobile gas mass fraction ݎ݂ can be inferred as a function of the saturation reached by the imbibition front:

(17) ܸ ሺܵ * ሻ ൌ ଶொ ் ೕ ఘ  ܵ • ௌ , ష ௌ * ߦ ௌ ሺܵ * ሻ • ௗక ೄ ௗௌ ሺܵ * ሻ݀ܵ (18) ݎ݂ ሺܵ * ሻ ൌ 2  ܵ • ௌ , ష ௌ * ߦ ௌ ሺܵ * ሻ • ௗక ೄ ௗௌ ሺܵ * ሻ݀ܵ.
Combining it with equation ( 15), we obtain this fraction as a function of the brine injection time.

Results on a reference test case

The following reference test case is considered: supercritical CO 2 is injected into a 20 m thick and 1200 m depth homogeneous aquifer (initially completely saturated with native water) through a vertical injection well that fully penetrates the reservoir. Properties and initial conditions of a deep aquifer likely to be suitable for CO 2 sequestration are provided in Table 1 mainly based on Pruess and Spycher paper (2007), with no salinity in our case.

[Table1]

The following scenario is considered: CO 2 is injected in the storage aquifer at a rate of 250 kt/y. An irreversible "abnormal behaviour" is assumed to occurr and to be detected after one year of injection (e.g. CO 2 escaping from the storage aquifer). To mitigate such a risk event, the following corrective action is implemented; water is injected at a rate of 160 kt/y (44 m 3 /h) in the view to enhance and accelerate gas trapping, i.e. to immobilize CO 2 by forcing the wetting process (see e.g. [START_REF] Manceau | Forcing gaseous CO 2 trapping as a corrective technique in the case of abnormal behavior of a deep saline aquifer storage[END_REF].

Please note that we considered Land's trapping model associated to Lenhard and Parker's hysteretic relative permeability functions [START_REF] Lenhard | A model for hysteretic constitutive relations governing multiphase flow -2. Permeability-saturation relations[END_REF].

Table 2 gives the values chosen for the trapped saturation and permeability functions computation.

[Table2]

Figure 4 depicts the analytical model results. It provides the analytical solutions for gas saturation profile after several water injection duration times and the temporal evolution of mobile gas fraction during water injection. The injected 250 kt of CO 2 are completely trapped after less than one year of water injection (Figure 4b).

[Fig. 4]

5 Discussion on the underlying assumptions

Fluid properties dependence on pressure

The flow is considered incompressible in Buckley-Leverett theory and the fluids density and viscosity are fixed. However, the pressure varies in a significant range (over time and space) when injecting CO 2 (see e.g. [START_REF] Zhou | A method for quick assessment of CO 2 storage capacity in closed and semi-closed saline formations[END_REF]) so that fluid parameters can change over a wide range of values during injection. Please note that the effects of CO 2 incompressibility approximation in analytical sharpinterface models have been recently explored by Vilarrasa et al. (2010). They found out that the error in the interface position due to this assumption is rather small when viscous forces dominate.

From a user perspective, only one couple density-viscosity is chosen in the analytical model. Several choices can be done:

the values for an estimated pressure that would correspond to a maximum or a mean value during both injections (CO 2 and water), but this requires the use of a pressure estimation model;

the values at initial pressure (we consider the temperature constant in that paper); this choice was made in section 4.2. These values can easily be computed from tables or formulas, but the viscosity and density may be underestimated, hence influencing the multiphase flow processes.

As the second choice seems more convenient, we explore in the following the influence on the results of such an assumption. To validate the analytical model on that point, i.e. to explore the effects of the spatial and temporal pressure variations, a comparison with numerical simulations that account for these variations is necessary. For this purpose, the multiphase fluid flow transport code TOUGH2 (Pruess et al. 1999) Lenhard and Parker's relative permeabilities [START_REF] Lenhard | A model for hysteretic constitutive relations governing multiphase flow -2. Permeability-saturation relations[END_REF].

For a proper comparison, please note that, in the numerical simulations:

the geometry and properties are the same than in the reference test case;

the solubility of components in both phases is not considered;

the gravity is neglected in the numerical simulations (1D horizontal mesh);

the salinity is set to zero;

the capillary pressure is neglected;

the flow is isothermal.

To investigate the influence of pressure variations (across space and time) a range of numerical simulations are carried out given different injection conditions, i.e. compressibility conditions: (1) 95 and 60 kt/y for CO 2 and water respectively ("low" injection rates); (2) 250 and 160 kt/y ("moderate" injection rates) and (3) 500 and 315 kt/y ("high" injection rates).

The saturation profiles after four months of water injection and the mobile gas mass fraction are compared to the solutions given by the analytical model for a density-viscosity couple corresponding to the initial pressure. The comparison shows very satisfactory match (see Figure 5 and Figure 6); the best fit is obtained logically at low rates when compressibility effects are the less significant (when the pressure is the closest to the initial one); the relatively good match is observed between the profiles even at high rates despite the deviations between densityviscosity values taken in the analytical model and the values computed in the numerical simulations (Figure 7 depicts these deviations just after the end of the CO 2 injection). More precisely, increasing the pressure does not result in strong viscous (see on Figure 8). In addition, Figure 8 shows larger trapped quantities for lower viscosity values.

According to this discussion, defining constant viscosity and density value at initial conditions may lead to an overestimation of the leading drainage front.

However, regarding the imbibition front and above all mobile gas mass fraction, viscosity and density-related processes act to compensate each other leading to little differences even for relatively large pressure differences.

[Fig. 8]

Influence of the capillary pressure

Including capillary effects in equation ( 1) amounts to add a diffusion-like term (see e.g. [START_REF] Pinder | Subsurface Hydrology[END_REF]. Practically, this will induce changes in the flow but also a smoothing of the sharp fronts (i.e. CO 2 and injected water fronts).

Integrating the capillary term in the fractional flow formula gives:

(19) ݂ ൌ ଵି ೖഊ ಲ • ങು ങ ଵା ഊ ഊ
, with ܲ ൌ ܲ െ ܲ . , is negative in the drainage part of the plume and positive in the imbibition zone, which means that capillary effects will lead to a gas flow increase in the drainage part and a reduction in the wetting one. In literature, these effects are mostly neglected when advective effects are enough significant so that the difference in capillary pressure during saturation evolution is usually considered low compared to pressure gradients. Errors due to this simplification will be maximized for high saturation gradient [START_REF] Juanes | Multiscale finite element methods for miscible and immiscible flow in porous media[END_REF]: at both fronts level, there is a capillary fringe where the differences in saturation profile with or without capillary effects get more significant; the fronts are then smoother with capillary effects. In particular, the smoothing of the imbibition front renders questionable the instantaneous trapping hypothesis explained in section 3.2 that leads to a more sharpened front.

To examine the effects of neglecting capillary pressure and of the additional instantaneous trapping approximation, we run numerical simulations of the application case including capillary pressure and for different injecting rates (the same than for the previous section). In the code used, the capillary pressure evolution is derived from van Genuchten's capillary pressure function in order to take into account hysteresis effects. The chosen capillary strength is 1 bar with a van Genuchten ݉ equal to 0.457.

Figure 9 and Figure 10 confirm the more significant influence of capillary effects for lower rates. The gas saturation profiles show that the analytical model underestimates the position of the leading CO 2 and water advancement fronts.

Quantitatively, the no-capillary pressure and instantaneous trapping approximations lead to more significant differences than the incompressible flow one. However, CO 2 injection operations and possible remediation measures in case of abnormal behaviour of the storage complex mostly imply moderate to high rates (i.e. industrial rates of the same order of magnitude of other underground activities such as deep geothermal activities) and hence, the simplifications made do not affect considerably the results, especially the trapped quantity evolution (Figure 10).

Summary and further works

An analytical Buckley-Leverett model has been set-up for CO 2 injection followed by a brine injection in a radially symmetric geometry. An analytical solution for saturation profile has been derived, which renders possible the assessment of the mobile gas remaining in the system as a function of the brine injection time. It then provides the quantity and the location of the mobile gas, which are both of primary importance when dealing with CCS safety.

This assessment method is flexible to the chosen trapping and hysteresis models.

It is rather simple and quick since, concerning the hysteresis model, it only requires the formula of the first drainage relative permeabilities functions. The solutions indeed do not depend on capillary pressure curves neither on relative permeabilities for imbibition and for the other cycles. This simplicity allows the implementation in spreadsheet type software.

The analytical solution is based on a range of simplifications and assumptions whose influence has been explored by means of the comparison with numerical simplification are low especially concerning the evolution of the mobile gas quantity; the capillary pressure influence logically decreases when injection rates increase and becomes negligible as soon as industrial rates are concerned.

Moreover, the instantaneous trapping approximation does not have major consequences on the results.

A second trapping mechanism plays an important role, namely the solubility trapping (CO 2 dissolution) and acts as an additional safety guaranty against gaseous CO 2 leakage (e.g. Ennis-King and Paterson 2005). Further developments of the analytical solution will then be focused on the integration of mutual solubility between phases.

Acknowledgements

The authors want to thank Arnaud Réveillère for his help to this work. This study has been The dissolution process is removed in the numerical simulator The original version is available at www.springerlink.com

  1997 and[START_REF] Medeiros | Hysteresis in two-phase flow: a simple mathematical model[END_REF], i.e. convex on part of the saturation interval and concave on the complementary one. It implies a physical impossibility in Buckley-Leverett equation, the derivative of the flux function having the same value for two different saturations. The problem can nevertheless be solved by means of a shock, thus introducing a discontinuity in the saturation profile

  zero (downstream of the advancing gas front, only native brine is present).[START_REF] Welge | A simplified method for computing oil recovery by gas or water drive[END_REF] sets-up a practical and graphical technique to find out the shock saturation as well as its velocity. As a matter of fact, the discontinuity (or advancing gas front) saturation can then be retrieved graphically plotting the line passing through the origin and tangent to the fractional flow curve; the front velocity is its slope. The rest of the gas saturation profile, obtained through Buckley-Leverett equation, then starts, at the vicinity of the injection point, from the residual liquid saturation and decreases progressively to reach the saturation front saturation.

  downward to the front, drainage is still going on and the saturation can be obtained by means of the drainage fractional flow curves; wetting is occurring upward and the saturation is obtained through imbibition fractional flow curves.Please note that the imbibition fractional curves depend on the saturation reached when the flow reversal occurs (corresponding to the downward saturation), which means that the velocity of the imbibition front depends on the reversal flow saturation as well: imbibition gas fractional flow function associated to a reversal saturation equal to ܵ ,௦ூூ ା

  imbibition curves which considerably simplifies the problem.

[

  Let us consider a saturation value ܵ. The first step consists in determining the brine injection time (߬), at which the brine front intercepts the saturation ܵ. The "dimensionless area" reached by the injected brine at ߬ is:(10) ߦ ሺ߬ሻ ଶ ൌ  ߞ ௦ூூ ሺݐሻ த ଵ.ݐ݀The dimensionless area reached by the gas saturation ܵ at ߬ is: a chosen time ߬ * , corresponding respectively to collision saturations 1 െ ܵ and ܵ * , we get the time at which the brine front reaches one specific saturation (here ܵ * ):

  including the EOS module ECO 2 N (Pruess 2005) and a hysteresis module (see Doughty 2009) is used. In this last module, Land's trapping model is considered with hysteretic characteristic functions derived from van Genuchten's capillary pressure function (van Genuchten 1980) and based on

[

FiguresFigure 1 :Figure 2 :Figure 3 :Figure 5 :Figure 7 :Figure 8 :Figure 9 :

 1235789 Figures

Figure 10 :

 10 Figure 10: Analytical and numerical (with capillary effects) results concerning the mobile gas mass fraction evolution for different injection rates. The dissolution process is removed in the numerical simulator

Table 1 :

 1 Initial conditions chosen in the simulations

	Initial pressure

differences for the imbibition front positions and mobile gas mass fraction while the gas advancement front seems more sensitive to those pressure effects.

[ Fig.5] [Fig.6] [Fig.7] In order to explain the results obtained, we explore the respective influence of density and viscosity changes (by decoupling them). Figure 8 depicts the analytical solution (for the case with "high" injection rates) with gas properties taken at initial aquifer pressure (120 bars, see Table 1), and with gas density and viscosity taken alternatively and then together at high pressure (180 bars, μ ൌ 7.04 • 10 ିହ Pa. s, ρ ൌ 790kg. m ିଷ ).

Keeping the same mass injection rate, changes in gas density induce changes in gas injection volumetric rate; then, at a given time and in the drainage zone, the saturation (and notably the leading one) will be closer to the injection well for higher values of density. However, this statement does not hold for the imbibition front: the saturation propagation in the drainage zones being more rapid when density is low, the water advancement front collide with higher gas saturations.

However, the velocity of this front depends on the reached saturation so that it is higher when the saturation is lower. The imbibition front will therefore reach a higher distance for a higher density value (see on Figure 8). According to equation ( 18), mobile gas mass fraction (as a function of the saturation reached or as a function of the dimensionless time) does not depend on density; however, coming back to dimensioned time requires gas density (through the injection rates ratio).

Concretely, the mobile gas mass fraction left after a given water injection period is less significant for higher density values (see on Figure 8).

Gas viscosity plays a role in the gas fractional flow as a gas with a lower viscosity flows more easily within water (i.e. fractional flow is higher). The velocity of the gas advancement front is higher for lower viscosity values (see on Figure 8). For the same quantity of gas injected (changing viscosity has no influence on the volumetric injection rates), the saturation at a given distance of the well is then lower with a lower viscosity (it is true for the whole profile except close to the advancement gas front). The imbibition front will be faster when the gas is less


Julien Bénabès

Emilie Poirson

Fouad Bennis

Yannick Ravaut

ICED11 1 MODULAR OPTIMIZATION STRATEGY FOR LAYOUT PROBLEMS

Keywords: layout problems, design optimization, modular optimization strategy, genetic algorithm

Layout design optimization has a significant impact in the design and use of many engineering products and systems. The search of an optimal layout configuration is a hard and critical task and designers need specific optimization strategies to achieve it. These strategies have to be generic to be adapted to a maximum of layout problems and a maximum of designer's requirements. These strategies have also to be effective in the search of optimal solutions without being too time consuming. Then, this paper proposes a modular optimization strategy based on the genetic algorithm. This strategy is adapted to multi-objective layout problems, when the designer has to make a final decision on optimal solutions. The optimization modules are used in order to make the search of feasible designs easier, by taking into account the specificities of layout problems. The different algorithms, which result from the modular optimization strategy, are tested and compared on a realworld application. This layout problem deals with the search of an optimal spatial arrangement of components inside a shelter.

INTRODUCTION

Layout problem is inherently a multidisciplinary task [START_REF] Giassi | Multidisciplinary design optimization and robust design approaches applied to concurrent design[END_REF]. It covers all the aspects of the product design life cycle from the conceptual to the detailed stage and makes necessary the collaboration between experts of technical and economical disciplines. Layout problems are usually considered as optimization problems. Although there are, in layout design literature, different definitions of layout optimization problems [START_REF] Cagan | A survey of computational approaches to the three-dimensional layout problems[END_REF], [START_REF] Yi | Vehicle configuration design with a packing genetic algorithm[END_REF], the key idea is always the same: given a set of free form components and an available space, a layout problem consists of finding the best arrangement (location and orientation) of components satisfying geometrical and functional constraints and achieving design objectives. This generic definition can be adapted to all real-world applications. For example, Drira et al. [START_REF] Drira | Facility layout problems: A survey[END_REF] and Wäscher et al. [START_REF] Wäscher | An improved typology of cutting and packing problems[END_REF] have adapted the definition of a layout problem to their respective research domain, meaning the facility layout design and the cutting and packing problems. Cagan et al. proposed in [START_REF] Cagan | A survey of computational approaches to the three-dimensional layout problems[END_REF] a schematic representation of the major constituent parts for solving a generic layout optimization problem. One finds in [START_REF] Bénabès | A new approach for specifying and solving layout problems[END_REF] another representation including three main parts: the description of the problem, the formulation of the problem and the optimization strategy. As it is shown in Figure 1, these three parts are connected. In fact, the problem description defines the dimension of the layout problem (one, two or three dimensions) and describes the layout components, meaning the container and the components that have to be placed into the container. Then, this description and all the expert's requirements have to be translated into design variables, constraints and objectives in order to change the layout problem into an optimization problem. During the optimization process, the solving strategy continuously interacts with the formulation problem in order to compute the design constraints and the objectives. The optimization strategy can also include an interactive tool that allows the designer to interact with the process. This tool needs a graphical environment where the layout components are modeled according to the problem dimension.

ICED11

Figure 1. Schematic representation of a layout problem

This paper focus on a particular step of the layout optimization process: the problem solving. In fact, layout problems are generally considered as non-linear and NP-hard optimization problems. Problems are intrinsically harder than those which can be solved by a non-deterministic Turing machine in polynomial time. One finds multiple single or multi-objective approaches to solve layout optimization problems in two or three dimensions [START_REF] Cagan | A survey of computational approaches to the three-dimensional layout problems[END_REF]. Traditional gradient-based approaches can be used for simple layout problems. For more complex real-world applications, some stochastic algorithms are required to avoid local optima. For example, some optimization strategies use genetic algorithms [START_REF] Yi | Vehicle configuration design with a packing genetic algorithm[END_REF], simulated-annealing algorithms [START_REF] Szykman | Constrained three-dimensional component layout using simulated annealing[END_REF], [START_REF] Shen | Three-dimensional design layout optimization for mechatronic devices using the simulated annealing approach[END_REF] or extended pattern search algorithms [START_REF] Su | An extended pattern search algorithm for three dimensional component layout[END_REF]. Most search algorithms are developed for a specific problem and they provide an effective optimization strategy for it. However, they are not generic and can not be adapted to a lot of layout problems. Some of recent studies deal with the search of efficient generic algorithms for solving layout problems. Jacquenot et al. propose in [START_REF] Jacquenot | 2d multi-objective placement algorithm for free-form components[END_REF] a hybrid algorithm based on a genetic algorithm coupled with a separation algorithm. A variant of this approach is also presented in [START_REF] Bénabès | Interactive optimization strategies for layout problems[END_REF]. In general, the development of an engineering object is considered as a single process involving multicriteria identification of the mathematical model followed by multi-criteria optimization of the object design on the basis of this mathematical model. The direct participation of the designer in the construction of the feasible design and non-formal analysis are the essential stages of the search for the optimal design. For solving the design problem, the designer almost always has to correct either the mathematical model, the dimension of the vectors of design variables and criteria, the design variable ranges, and so on. One can find in [START_REF] Michalek | Interactive design optimization of architectural layouts[END_REF] a significant contribution to this concept applied to the design optimization of architectural layouts. Moreover, interactivity with designer can be used to insert qualitative fitness or user perceptions in the design process [START_REF] Poirson | Integration of user perceptions in the design process: application to musical instrument optimization[END_REF]. In layout design, Brintrup et al. have already developed an interactive genetic algorithm based framework for handling qualitative criteria in design optimization [START_REF] Brintrup | An interactive genetic algorithm-based framework for handling qualitative criteria in design optimization[END_REF]. Actually, this paper proposes a modular optimization strategy, based on the genetic algorithm. This optimization approach consists in combining a genetic algorithm with optimization modules, in order to help the genetic algorithm in the search of designs which respect the design constraints. This strategy is tested on the three dimensional layout problem of a shelter. This paper is organized as follows: the section 2 presents the application studied in this paper. It presents the description and the formulation of the layout problem of a shelter. The section 3 explains why the genetic algorithm is used in the optimization approach proposed in this paper, by defining in particular the complexity of a layout problem. This section also presents the different algorithms which results from the modular approach. The section 4 is dedicated to comparison of the performances of the optimization algorithms, which result from the modular approach, on the layout problem of the shelter. The section 5 concludes this study.

LAYOUT PROBLEM OF A SHELTER

The application studied in this paper deals with the layout problem of a shelter. This problem is a three dimensional layout problem. Eight components have to be located in the shelter, including four electrical and energy cabinets, two desks and two electrical boxes. The CAD model of the shelter is presented in Figure 2.

Figure 2. 3D CAD model of the shelter

This problem has been simplified in two dimensions, according to the fact that the cabinets are full height of the shelter and prevent a superposition of components. One of the possible designs of the shelter, modeled in two dimensions, is presented in Figure 3. The description and the formulation of this model are also described in [START_REF] Bénabès | A new approach for specifying and solving layout problems[END_REF] and [START_REF] Bénabès | Interactive optimization strategies for layout problems[END_REF].

Figure 3. Description of the shelter in two dimensions

The layout components are divided in two categories of components: the material and the virtual components. A material component has a mass and can not overlap with another material component. On the other hand, a virtual component has no mass and can overlap with a material or a virtual component, according to the designer's requirements. In fact, in the layout problem of the shelter, the layout components are made up of: • 8 material components: 4 cabinets, 2 desks and 2 electrical boxes • 6 spaces of accessibility (dotted rectangles in Figure 3) • 1 free space (hatched rectangle in Figure 3) located below the air-conditioner, where no cabinet can be placed • 1 free space in front of the door • 1 free corridor located in the middle of the shelter A space of accessibility is linked to a material component. For example, the space of accessibility of the cabinet 1 is defined as the required space, placed in front of the cabinet 1, and used to insert some materials into the cabinet 1. No material component can be placed in the space of accessibility. The free corridor is used to guarantee that all components are accessible from the shelter's entry. This free corridor, the other free spaces and the spaces of accessibility can be considered as virtual components.

The density of this layout problem, defined as the ratio between the space occupied by the components and the available space in the container, is equal to 105%. This density does not take into account that a virtual component can overlap with another virtual component. C4: non-overlap constraints between the cabinets and the free space located below the airconditioner Because of the dimension of this problem and the rectangular shape of the layout components, the overlap is defined as the intersection area between components. The intersection area between the components i and j is defined by:

)] 2 max( ) 2 min( , 0 max[ )] 2 max( ) 2 min( , 0 max[ j j i i j j i i ij L y L y l x l x A - - + × - - + = (1)
where ( i x , i y ) are the coordinates of the center of gravity of the rectangle i. i L and i l are respectively the length and the width of the rectangle i . Two objectives

) 2 max , 1 (min O O
are also defined. The objective 1 O is used to balance the masses inside the shelter. It means to minimize the distance between the center of gravity of the layout components and the geometrical center of the shelter. The objective 2

O is used to maximize the distance between the cabinet 1 (energy network) and the cabinets 2 and 3 and the electrical box 2 (electrical network). All the distances between components, used in the computation of the two objectives, are the distances between the geometrical centers of components defined by the Euclidean norm. More information about the definition of these objectives is available in [START_REF] Bénabès | Interactive optimization strategies for layout problems[END_REF].

MODULAR OPTIMIZATION STRATEGY

The principle of this modular optimization strategy is to use a genetic algorithm and to insert in this algorithm some optimization modules in order to improve the global performances of the optimization strategy. This section presents the different algorithms which result from this modular approach.

Genetic algorithm

This paper proposes to use the genetic algorithm process because of the great complexity of real-world layout problems. The complexity of a layout problem is defined in [START_REF] Bénabès | A new approach for specifying and solving layout problems[END_REF]. This complexity is linked to the geometry of layout components, the layout density and the problem formulation. The problem formulation has a significant impact on the choice of an appropriate optimization algorithm. More numerous the constraints and objectives are, more complex the search of a feasible design is. The design space is parceled and the designer can not use traditional gradient-based optimization approaches to pass to a feasible region to another one. Stochastic algorithms, as the genetic algorithm, have to be used and the calculation time increases so that the problem is more complex. In this paper, the Genetic Algorithm Omni-Optimizer is used [START_REF] Deb | Omni-optimizer: A generic evolutionary algorithm for single and multiobjective optimization[END_REF]. This algorithm is designed to handle single and multi-objective problems. Given a set of initial individuals, randomly generated, the genetic algorithm uses basically three operators in order to create a set of new solutions. These genetic operators are the selection, the crossover and the mutation. Figure 4 illustrates the convergence of the Genetic Algorithm Omni-Optimizer for the layout problem of the shelter. The axis of ordinates represents, for each generation of the genetic algorithm, the minimal sum of constraints obtained by one design. It means that, for a specific generation, if this minimal sum of constraints is equal to 0, the genetic algorithm has found a "feasible" design, meaning a design which respects all the non-overlap constraints. This algorithm has been initialized by 200 designs, randomly created. Figure 4 shows that the genetic algorithm converges after 60 generations but it can not find a feasible solution, which respects all the placement constraints. The layout problem of the shelter is too complex for this genetic algorithm. Consequently, the main objective of the following sub-sections is to describe some optimization modules, which can be linked to the genetic algorithm in order to help it in the search of feasible designs. These optimization modules are specific to layout optimization problems.

Module 1: "Optimization of the direction of the components"

In lots of real-world applications, as the layout problem of the shelter, one can find material and virtual components. This optimization module is specific to these layout problems, which deals with the spatial arrangement of material and virtual components. Let us consider the application studied in this paper. The shelter is modeled in two dimensions and all the layout components have a rectangular shape. Let us consider a material component 1 (for example the cabinet 1) and a material component 2 (for example the electrical box 2). A virtual component (for example the space of accessibility of the cabinet 1), represented by a dotted rectangle in Figure 5, is linked to the component 1. Figure 5 illustrates the two possible directions ) (λ for the component 1, meaning the two possible positions of the virtual component. The objective of the module 1 "Optimization of the direction of the components" is to optimize the placement of this virtual component, by deleting the discrete variable λ used for the direction and minimizing the non-overlap constraints betweens virtual components and material components. For example, Figure 5 shows that the module 1 chooses the case (a), with a direction equal to 1 because this direction minimizes the non-overlap constraint between the space of accessibility of the component 1 and the component 2.

Module 2: "Separation algorithm"

Given a layout configuration that does not satisfy placement constraints, the objective of the separation algorithm is to minimize the non-respect of overlap and protrusion constraints. In fact, the separation algorithm is designed to solve a single objective optimization problem, where all the placement constraints are gathered into one objective function. This separation algorithm has been tested on a simple two dimensional layout problem, which deals with the search of an optimal spatial arrangement of square components inside a square container. Different optimization problems with different densities have been tested and the results are presented in [START_REF] Jacquenot | 2d multi-objective placement algorithm for free-form components[END_REF].

Module 3: "Local disruption"

This module is used to help the separation algorithm to find a feasible design, meaning a design which respects all the placement constraints. In fact, the separation algorithm, proposed in this paper, only modifies the position variables, it means the continued variables. It has not effect on the orientation or direction variables. Consequently, for some layout configurations, the separation algorithm does not find a feasible solution. It finds a local minimum which does not respect all the placement constraints. One or several components have to be rotated. The main idea of this module 3 is to associate to the separation algorithm a local disruption in the layout configuration. This local disruption is used to randomly change the orientation of some components. The components which are rotated are randomly selected. The number of rotated components is fixed by the designer. For example, it can be fixed to 30% of the total number of components of the layout problem. In fact, a local "disruption" already exists in the genetic algorithm process. This disruption is realized by the genetic operator "mutation". However, the mutation is not dedicated to specific optimization variables. The module 3 is dedicated to the variables which define the orientation of components. Consequently, this optimization module is more effective in the search of feasible designs.

Modular optimization strategy

This sub-section describes the different algorithms which results from the modular optimization approach proposed in this paper. Figure 6 illustrates two optimization algorithms, created by combining:

• Figure 6 (a) : the genetic algorithm with the modules 1 and 2 • Figure 6 (b) : the genetic algorithm with the modules 1, 2 and 3 The module 3 can be repeated max j times, where max j is fixed by the designer. Between two applications of the module 3, the optimization process tests if the current design respects all the placement constraints. If so, the optimization process evaluates the design and moves to the next solution according to the genetic operators. Otherwise, the module 3 is run if

max j j < . Moreover, if max j j =
and if the separation algorithm has not found a feasible solution, the optimization process continues by considering the design that minimizes the placement constraints among all the designs generated by the separation algorithm. Moreover, a hybrid algorithm, used for solving layout problems, is described in [START_REF] Su | An extended pattern search algorithm for three dimensional component layout[END_REF]. This algorithm is obtained by combining the genetic algorithm with the optimization module 2. Firstly, Figure 7 shows that all the optimization algorithms have converged because the algorithms are not able to find more variants. For each algorithm, the maximal number of iterations has also been attained on the period of calculation time defined by the graph. Secondly, the results underline the great influence of the module 1 on the performances of the algorithms. In fact, the algorithm A has generated 9 variants whereas the genetic algorithm without optimization modules was not able to find feasible solutions. By the same way, the results also underline the great influence of the module 2 and the combination of the two optimization modules 1 and 2. The optimization algorithm C found more variants than the algorithm B, which found more variants than the algorithm A. Moreover, Figure 7 shows that the optimization modules 1 and 2 allow the algorithm to generate variants more quickly. On the other hand, it can be relevant to evaluate the different Pareto-optimal variants, generated by these three algorithms, according to their objective values. Let us plot a scatter graph with all the Pareto-optimal solutions found by the three algorithms. Figure 8 illustrates this scatter graph. Figure 8 shows that the Pareto-optimal solutions generated by the algorithm C have better performances than the other solutions. In summary, the optimization algorithm C, based on a genetic algorithm combined with the modules 1 and 2 has the best performances. Let us now consider the algorithms C, D, E and F. By the same way as Figure 7, Figure 9 illustrates the results obtained by these four algorithms. Firstly, in this simulation, on the period of calculation time defined by the graph, the maximal number of iterations is only attained for the algorithms C and D. In fact, the algorithms E and F use the optimization module 3, with the parameter max j respectively equal to 10 and 20, and this optimization module is time consuming. When the module 3 is run, the module 1 and 2 are also run. Each generation of the algorithm takes more time to generate designs and the maximal number of generations is not attained on the period of calculation time defined by the graph. Secondly, it is difficult here to define which optimization algorithm has the best performances because it depends on the designer's requirements. The optimization algorithms can be sorted by ascend with the number of generated variants in the order C, D, E and F. It means that the module 3 has a great influence on the search of feasible variants because the local disruption helps the separation algorithm in the search of feasible designs. However, the algorithm C finds variants before the other optimization algorithms. Consequently, the choice of the optimization algorithm depends on the designer's requirements, if he wants to find lots of solutions without caring about the calculation time or if he wants to generate few solutions in a short period of calculation time. Moreover, the Pareto-optimal variants generated by the optimization algorithms have similar performances.

CONCLUSION

Layout problems are generally considered as complex problems and stochastic optimization algorithms have to be used. The complexity of layout problems is usually linked to the problem formulation. The design space is parcelled so that the designer needs specific optimization algorithms as the genetic algorithm in order to find feasible solutions. However, the application studied in this paper shows that the genetic algorithm is not completely adapted to very complex layout problems. The layout problem of the shelter deals with the search of an optimal spatial arrangement of material and virtual components. This problem has a big density and its formulation is complex. Consequently, this paper proposes a modular optimization approach, based on the genetic algorithm, in order to help the algorithm in the search of feasible designs. Three optimization modules are described: the optimization of the direction of the components, the separation algorithm and the local disruption. The designer can decide to insert these modules in the genetic algorithm process, according ICED11 to the specificities of his layout problem. For example, the module 1 is only dedicated to layout problems which deal with the spatial arrangement of material and virtual components. The results described in this paper show that the optimization modules enable the designer to generate more feasible variants, meaning more designs which respect all the placement constraints. The optimization modules are time consuming. The designer has to choose the most appropriate module according to his requirements: to quickly generate variants with average performances or to generate lots of variants with good performances. In the two scenarios, the calculation time is not the same. The development of this modular optimization approach is part a global design approach used to solve complex layout problems. These optimization algorithms are dedicated to multi-objective optimization problems. The number of variants, generated by these algorithms is very big. However, the designer has to make a final choice, by choosing one of these optimal variants. Outlooks are now dedicated of the development of methods and tools which enable the designer to make this final choice. These tools provide interaction between the designer and the optimization process, in order to locally modify a design, modify the problem formulation or insert the personal judgment of the designer in the search of optimal solutions.

Figure 4 .

 4 Figure 4. Convergence of the genetic algorithm

Figure 5 .

 5 Figure 5. Direction of a layout component

Figure 7 .

 7 Figure 7. Comparison of the algorithms A, B and C

Figure 8 .

 8 Figure 8. Pareto-optimal variants for the algorithms A, B and C

Figure 9 .

 9 Figure 9. Comparison of the algorithms C, D, E, and F

ACKNOWLEDGMENTThe authors would like to acknowledge Thales Communications for the application study.

In order to compare the performances of the different algorithms, the number of feasible variants generated by each algorithm is measured for each generation. In fact, the design j is a new variant if it differs from the design i by at least one of the following criteria:

• one of the components of the layout has been displaced from at least Δ mm along one of the axis X and Y ( Δ is fixed to 500mm in the following simulations) • one of the components has been rotated • the minimum difference between the objective values of the two designs is bigger than a limit, fixed to 10 cm This number of variants is plotted according to the calculation time. By associating the different modules with the genetic algorithm, with different parameters, six optimization algorithms are created: • optimization algorithm A: genetic algorithm with the module 1 • optimization algorithm B: genetic algorithm with the module 2 • optimization algorithm C: genetic algorithm with the modules 1 and 2 • optimization algorithm D: genetic algorithm with the modules 1, 2 and 3, with max j =3 • optimization algorithm E: genetic algorithm with the modules 1, 2 and 3, with max j =10 • optimization algorithm F: genetic algorithm with the modules 1, 2 and 3, with max j =20 The parameters of the genetic algorithm Omni-Optimizer are the same for each simulation. The size of the population is fixed to 200 individuals. The maximal number of iterations is fixed to 100. Firstly, let us consider the comparison of the optimization algorithms A, B and C. Figure 7 illustrates the results obtained for these three algorithms.