Charlotte Thierion
email: charlotte.thierion@mines-paristech.fr

Laurent Longuevergne

Florence Habets

Emmanuel Ledoux

Philippe Ackerer

Samer Majdalani

Etienne Leblois

Simon Lecluse

Eric Martin

Solen Queguiner

Solen Queguiner

Pascal Viennot

Assessing the water balance of the Upper Rhine Graben hydrosystem

Keywords: groundwater-surface water relations, groundwater recharge/water budget, Upper Rhine alluvial aquifer, hydrogeological model parameters sensitivity

come   L'archive ouverte pluridisciplinaire

Introduction

Alluvial hydrosystems have been recognised as important water resources but also very vulnerable systems to pollution or change in water availability [START_REF] Allen | Groundwater and climate change: a sensitivity analysis for the Grand Forks aquifer, southern British Columbia, Canada[END_REF].

Knowledge of the different components of the water balance is therefore of significant importance to optimize water management and/or to estimate the impacts of anthropogenic activities or climate changes. These water balance components are difficult to measure, as the observations of water fluxes provide only local information [START_REF] Sanford | Recharge and groundwater models: an overview[END_REF][START_REF] Kalbus | Measuring methods for groundwater-surface water interactions: a review[END_REF]. Thus there is a need to assess these fluxes by modelling, as they are linked to the boundary conditions, the sink/source terms and the hydrodynamic parameters, which are mostly unknown and estimated through model calibration. Model calibration leads to non-unique solutions and therefore, the evaluation of some components of the water balance is uncertain [START_REF] Konikow | Ground-water models cannot be validated[END_REF]Bredehoeft 1992, Beven 2006).

The components of the water balance for unconfined groundwater are recharge due to precipitations, exchanges with surface water, in/out-flows through the boundaries (either prescribed heads or prescribed fluxes), storage in the aquifer, and water uptake by pumping wells. Recharge and surface-groundwater interactions constitute key processes in the behaviour and the evolution of these systems [START_REF] Sanford | Recharge and groundwater models: an overview[END_REF][START_REF] Woessner | Stream and Fluvial Plain Groundwater Interactions : Rescaling Hydrogeologic Thought[END_REF][START_REF] Sophocleous | Interactions between groundwater and surface water : the state of the science[END_REF]. These fluxes can be estimated by fully coupled surface-subsurface models including the unsaturated zone such as ParFlow [START_REF] Kollet | Integrated surface-groundwater flow modeling: A freesurface overland flow boundary condition in a parallel ground water flow model[END_REF], Hydrogeosphere [START_REF] Therrien | HydroGeoSphere: A Three-dimensional Numerical Model Describing Fully-integrated Subsurface and Surface Flow and Solute Transport[END_REF], and CATHY (Camporese et al., 2008) among others. However, applications of these fully coupled models are generally limited to small areas. For larger systems, coupled surface-subsurface models are usually simplified [START_REF] Hu | Simulated groundwater interaction with rivers and springs in the Heihe river basin[END_REF][START_REF] Rushton | Representation in regional models of saturated river-aquifer interaction for gaining/losing rivers[END_REF]) and the recharge is estimated separately using simplified hydrological models based on 1D Richards equation (e.g. MIKE SHE, Reefsgard 1997) or more conceptual schemes (e.g. MARTHE, [START_REF] Noyer | Modelling aquifer salinity in the Potash Basin (Alsace)[END_REF].

This work aims at estimating the different components of the water balance and its associated uncertainties for the Upper Rhine aquifer. The studied hydrosystem is one of the most important water resources in Western Europe, and several regional scale models were previously developed, as for example in the LIFE project in 1996 [START_REF] Lfu | Action de démonstration portant sur la protection et la gestion des réserves en eau souterraine dans la partie franco-germano-suisse de la vallée du Rhin Supérieur -Rapport final (Demonstration action on protection and management of groundwater resources in the french-german-swiss part of the Upper Rhine valley -Final report[END_REF][START_REF] Lfu | MONIT: Développement des outils de prévision (MONIT : developement of prevision tools[END_REF], and in the INTERREG project Monit [START_REF] Lubw | Projet INTERREG III Monit : Modelisation hydrodynamique et transport des nitrates (INTERREG III project Monit : Hydrodynamic modelling and nitrates transport)[END_REF] focusing on nitrate contamination of groundwater. These two studies led to different findings about aquifer recharge: river infiltration was found to be the main recharge component in LIFE, whereas in Monit it was effective rainfall. However, these discrepancies may be linked to the fact that the models were limited to the alluvial plain, and used simple methods to assess the subsurface flows from the Vosges and the Black Forest surrounding mountains.

In the present work a modelling of the whole Rhine basin between Basel and Lauterbourg is performed (figure 1). This enables to better constrain the lateral water inputs to the aquifer, and thus, to have a homogenous representation of the whole Upper Rhine graben basin.

Moreover it allows assessing the impacts of climate change on the hydrosystem, which could lead to altered conditions over the mountainous catchments.

The Upper Rhine alluvial aquifer groundwater dynamics is modelled using different assumptions and calibration methods based on trial and errors or inverse approach. The use of several likely parameter values for a given model leads to different likely water balances and provides an estimate of the associated uncertainty. Since several possible calibrated models are obtained, we associated classic statistical criteria (average error, Nash criteria) with a more advanced statistical method, the Karhunen Loève transform (KLT) to evaluate the model performance [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences[END_REF][START_REF] Longuevergne | Extracting coherent regionel information from local measurements with Karhunen-Loève transform : case study of an alluvial aquifer (Rhine Valley, France and Germany)[END_REF]. The KLT allows the decomposition of piezometric time series into several independent temporal vectors, corresponding to the main influential processes. This gives a different insight into the performances of the simulations, focused on the representation of the main influential processes over the piezometric head temporal evolution, and helps in discriminating and quantifying the model ability to represent these processes. Using the results of this analysis, the aquifer water budget and recharge components are assessed, as well as the uncertainty in their computation. This article is organised as follows: first, the characteristics of the basin are presented. Then, the models used in this work are described (section 3), and the calibration of the surface water budget is presented (section 4). In section 5, the different parameter sets are presented, and the sensitivity of model results is analysed using the classic statistical criteria and the Karhunen Loève transform. Eventually, section 6 presents the estimation of the aquifer water budget and stock variations and the uncertainty on the recharge/discharge fluxes.

Basin characteristics

Geographical and geological setting

The Upper Rhine graben hydrosystem, situated at the French-German border (Figure 1), is composed of three units: i) a tectonic graben containing Secondary to Tertiary sediments, several kilometres thick (Figure 2) and forming the substratum of the Rhine alluvial plain [START_REF] Illies | The Rhine graben rift system-plate tectonics and transform faulting[END_REF], ii) the graben shoulders with much older materials such as crystalline, metamorphic and ancient sedimentary rocks constituting the Black Forest mountains in the East and the Vosges mountains in the West, and iii) the Sundgau hills in the South, formed by marly Oligocene substratum of the Quaternary sediments which was relatively less subsided than in the rest of the graben during the continental rifting.

The graben contains Tertiary marls and clays, covered by Quaternary alluvium deposited by the Rhine River and forming the studied alluvial aquifer (cf. figure 2) [START_REF] Duprat | Sciences Geologiques -La nappe phreatique de la plaine du Rhin en Alsace (The phreatic water table of the Rhine Plain in Alsace)[END_REF][START_REF] Bauer | Isotopic and hydrochemical studies of groundwater flow and salinity in the Southern Upper Rhine Graben[END_REF]. The hydrographic network is very dense in the plain, due to the presence of many groundwater-fed streams, allowing significant exchange of water between rivers and the aquifer [START_REF] Schmitt | Typologie hydro-geomorpholoqique fonctionnelle de cours d'eau : recherche methodologique appliquee aux systemes fluviaux d'Alsace (Hydrogeomorphologic functionnal typology of streams : methodologic research applied to the fluvial systems in Alsace)[END_REF]. The main Rhine's tributary in this part of its basin, the Ill River (figure 1), has its source in the Sundgau, while some small rivers having their sources in the same area fully infiltrate towards the aquifer. The other most important rivers have their source in the Vosges and Black Forest mountains West and East of the plain.

The material of the Rhine alluvial aquifer is mainly coarse quaternary gravels and sands with good hydrogeological properties. The hydraulic conductivities are in the order of 10 -4 -10 - 3 m/s [START_REF] Duprat | Sciences Geologiques -La nappe phreatique de la plaine du Rhin en Alsace (The phreatic water table of the Rhine Plain in Alsace)[END_REF][START_REF] Lfu | Action de démonstration portant sur la protection et la gestion des réserves en eau souterraine dans la partie franco-germano-suisse de la vallée du Rhin Supérieur -Rapport final (Demonstration action on protection and management of groundwater resources in the french-german-swiss part of the Upper Rhine valley -Final report[END_REF][START_REF] Lubw | Projet INTERREG III Monit : Structure hydrogéologique et caractéristiques hydrauliques (INTERREG III project Monit : Hydrogeological structure and hydraulic characteristics[END_REF]. On the other hand the rocks forming the mountains are rather impervious material containing only local aquifer formations of small extension. The alluvial aquifer thickness reaches more than 200 m in its centre, East of Colmar [START_REF] Lubw | Projet INTERREG III Monit : Modelisation hydrodynamique et transport des nitrates (INTERREG III project Monit : Hydrodynamic modelling and nitrates transport)[END_REF]. Mineral composition and thus structure and texture of the sediments deposited by the Rhine River are marked by their Alpine origin, whereas towards the edges of the graben, Vosgian and Black Forest rivers have deposited alluvial material having different hydrogeological properties. These elements contribute to the heterogeneity of the aquifer material, which in turn impacts the geometry of the river-aquifer interactions.

The aquifer substratum consists of the Oligocene marls considered impervious. The groundwater flows from South to North and is closer to the surface in its northern part. Water table depths range from 0 to 20 m (Hardt area), and groundwater fed wetlands are present in the middle and the northern parts of the plain.

Climatic conditions

Precipitations show large contrasts over the basin, with annual values ranging from 550 mm/year in the plain, to more than 2000 mm/year in the Vosges and Black Forest mountains (Figure 1). Snowfalls account for around 3 % of total precipitations in the plain, and up to 37 % at the mountains tops. Thus snow accumulation and melting are important processes for the dynamics of river flows in these catchments. Moreover potential evapotranspiration, computed with Penman formula on the 1983-2006 time period, is weaker on the mountains tops (around 600 mm/year) than in the plain (around 800 mm/year). As a consequence the rivers flowing from the mountainous catchments carry an important proportion of the water involved in the hydrosystem budget and play a crucial role in the recharge of the Rhine alluvial aquifer.

Daily precipitations and potential evapotranspiration (PET) for the modelling are provided by the SAFRAN analysis performed by Météo France, based on the definition of climatologically homogeneous zones and on altitudinal variation of the meteorological variables (Quintana-Seguì et al 2008, Vidal et al., 2010). They are used here at a daily time step from 1985 to 2003, which is about the period of time used in the previous Monit project [START_REF] Lubw | Projet INTERREG III Monit : Modelisation hydrodynamique et transport des nitrates (INTERREG III project Monit : Hydrodynamic modelling and nitrates transport)[END_REF]. For the present study the SAFRAN analysis that is usually provided on a regular grid of 8x8 km 2 was projected over a specific grid with cells varying from 1 to 8 km size in order to better account for the topographic gradient of precipitation (Figure 1).

Available data

Daily river discharges are observed in about 20 to 30 river gauges within the studied area, depending on the time period. These gauging stations are managed mainly by the DREAL-Alsace (Direction Régionale de l'Evironnement, de l'Aménagement et du Logement in Alsace) on the French side, and by the LUBW-Baden-Württemberg (Landesanstalt für Umwelt, Messungen und Naturschutz in Baden-Württemberg) on the German side.

Piezometric levels data consist in numerous observation wells (more than 200 points) sampled weekly, only a few of them having daily data. These observation wells are managed mainly by the APRONA (Association pour la PROtection de la Nappe d'Alsace) on the French side and by the LUBW on the German side.

Seasonal piezometric levels fluctuations are generally limited to one meter or less, except in the South and near the Eastern and Western borders where the water table is deeper. No long term tendency appears for the piezometric levels over the simulated period, from 1985 to 2003. Nevertheless piezometric levels are rather low between 1990 and 1995, and rather high between 2000 and 2003. The same holds true for the mountain rivers discharges (cf. figure 7).

Land use is characterized by large agricultural areas in the plain, with extensive cultures of cereals (some of them being irrigated), vineyards in the piedmont, whereas the mountains are largely covered by coniferous and deciduous forests (Figure 3). Agricultural lands cover approximately 40 %, and forests 45 % of the basin. Artificial areas represent almost 10 % of the area, which is densely populated with more than 200 inhabitants/km². The main cities are Basel, Colmar, Mulhouse and Freiburg in the South, Strasbourg in the North. Important water withdrawals, around 3.6.10 8 m 3 /yr, are located close to these cities for drinking and industrial purpose. Irrigated agricultural areas represent around 6 % of the region, mostly in the Southern central part of the plain.

Models presentation

Two different models were used: one classical hydrogeological model (MODCOU, [START_REF] Ledoux | Spatially distributed modeling: conceptual approach, coupling surface water and groundwater[END_REF], and one hydrogeological model including an inversion of the model parameters capacity (HPP-INV, [START_REF] Majdalani | Identification of Groundwater Parameters Using an Adaptative Multiscale Method. Application to the Upper Rhine aquifer[END_REF].

In this section, the models are presented briefly, with a special focus on MODCOU which is mainly used in this study.

MODCOU

MODCOU is a spatially distributed coupled hydrogeological modelling tool that was already used in several large basins in France [START_REF] Golaz-Cavazzi | Comparison of two hydrological simulations of the Rhone basin[END_REF][START_REF] Ledoux | Agriculture and groundwater nitrate contamination in the Seine basin. The STICS-MODCOU modelling chain[END_REF][START_REF] Korkmaz | Application of the coupled model to the Somme river basin[END_REF].

It is built up on several interconnected modules. The first one is dedicated to the estimation of the surface water budget computed with a simple reservoir scheme, the production functions (Figure 4). Inputs to this conceptual model are the daily precipitation and potential evapotranspiration (PET) data on a spatial grid. Outputs are actual evapotranspiration (AET), infiltration and runoff on each grid cell and at a daily time step. More details on the production functions are given in the next section where the calibration of the surface water budget is presented.

The direct infiltration of water from effective rainfall through the unsaturated zone is simulated with a Nash cascade [START_REF] Philippe | Improvement of the solute transfer in a conceptual unsaturated zone scheme: a case study of the Seine River Basin[END_REF].

The aquifer representation is based on a square grid discretization, with a finite differences resolution of the diffusivity equation. The rivers flows are simulated according to a Muskingum scheme [START_REF] David | Rapid applied to the SIM-France model[END_REF].

There are two kinds of interactions between the surface and the groundwater: 1) the aquifer can be drained where piezometric level reaches the soil level, and 2) drainage or infiltration of water towards the aquifer can occur on river cells. Owing to the major role of the river-aquifer interactions on the Upper Rhine Graben aquifer, the description of these processes in MODCOU is further detailed here.

On river cells, the river-aquifer exchange flow is proportional to the head difference between the river and the aquifer, to permeability and width of the riverbed, and inversely proportional to its thickness. As these characteristics of the river bed are mostly unknown, they are integrated in two parameters called the transfer coefficient (T p ), for the case when the aquifer and the river are connected, and the maximum river infiltration flow (Q lim ) for the disconnected case [START_REF] Rushton | Representation in regional models of saturated river-aquifer interaction for gaining/losing rivers[END_REF]. The exchange flow is therefore computed as follows:

    dis riv gw P exch Q Q H H T Q     ; ; max lim , (1)
With Q exch the exchange flow between the river and the aquifer which is negative when water flows from the river to the aquifer, H gw the piezometric head in the aquifer, H riv the water head in the river, and Q dis the available flow in the river cell. Q lim is a maximum infiltration flow from the river to the aquifer, corresponding to the case when the aquifer head is below the river bed [START_REF] Rushton | Representation in regional models of saturated river-aquifer interaction for gaining/losing rivers[END_REF]. A steady river head H riv is usually considered in MODCOU, however, due to the impact of river-aquifer interactions in the Rhine hydrosystem, time varying heads based on observations were included for the Rhine River, as well as for the Hardt Canal. This waterway is used for irrigation purposes, and preliminary results have shown it was important to explicitly simulate its losses to the aquifer.

MODCOU uses a daily time step to allow a better estimation of the river-aquifer exchange dynamics.

MODCOU can also be coupled with the soil-vegetation-atmosphere transfer scheme ISBA (Noilhan and Planton 1989) to form the SIM model [START_REF] Habets | The SAFRAN-ISBA-MODCOU hydrometeorological model applied over France[END_REF]. SIM is used here to have another estimation of the effective rainfall infiltration over the aquifer.

HPP-INV

The HPP-INV model (Chardigny 1999, Majdalani and[START_REF] Majdalani | Identification of Groundwater Parameters Using an Adaptative Multiscale Method. Application to the Upper Rhine aquifer[END_REF]) is a modelling tool using an inversion algorithm in order to assess aquifer hydrodynamic parameters. It is based on a classical hydrogeological model, and includes an adaptative multiscale triangulation for the inversion algorithm [START_REF] Majdalani | Identification of Groundwater Parameters Using an Adaptative Multiscale Method. Application to the Upper Rhine aquifer[END_REF].

Piezometric levels simulation is based on a finite elements scheme and the river flows are simulated with a Muskingum scheme [START_REF] Majdalani | Identification of Groundwater Parameters Using an Adaptative Multiscale Method. Application to the Upper Rhine aquifer[END_REF]. Rivers-aquifer interactions are computed from the difference between river and aquifer heads, with transfer coefficients values calibrated during the inversion process [START_REF] Lubw | Projet INTERREG III Monit : Modelisation hydrodynamique et transport des nitrates (INTERREG III project Monit : Hydrodynamic modelling and nitrates transport)[END_REF]. In the present application, this model is run at a monthly time step.

There is no computation of a surface water budget, and groundwater recharge by rainfall is provided by MODCOU. Moreover, the model is limited to the aquifer in the plain, and the lateral inflows from bordering catchments are calibrated.

The inversion method used to derive the spatial distribution of permeability is based on minimization of the error between simulated and observed piezometric heads and river flow data. Several distributions can be obtained according to different inversion conditions.

Model set up

Discretization

The area modelled with MODCOU corresponds to the Rhine basin between Basel in the South and Lauterbourg in the North. The basin total area is approximately 13900 km², and the area of the alluvial aquifer is 4655 km² (Figure 3). Two different boundary conditions are imposed: the Rhine river discharge at Basel, and fixed piezometric heads at the northern and southern boundaries.

The basin was discretized with square cells ranging in size from 200 m for the river cells to 1600 m. There are 99715 surface cells, and 34180 aquifer cells (Figure 3). The alluvial aquifer is represented with a single layer.

Subsurface flow at the western and eastern borders of the aquifer is known to be a significant part of the recharge [START_REF] Lfu | Action de démonstration portant sur la protection et la gestion des réserves en eau souterraine dans la partie franco-germano-suisse de la vallée du Rhin Supérieur -Rapport final (Demonstration action on protection and management of groundwater resources in the french-german-swiss part of the Upper Rhine valley -Final report[END_REF][START_REF] Lubw | Projet INTERREG III Monit : Modelisation hydrodynamique et transport des nitrates (INTERREG III project Monit : Hydrodynamic modelling and nitrates transport)[END_REF]. It is represented by specific boundary conditions in the model: the surface runoff reaching the aquifer surface on the bordering cells infiltrates to the water table and thus contributes to the aquifer recharge. This process is called hereafter lateral subsurface flow. Uncertainty on the components of the water budget is evaluated by exploring several equally probable hydrodynamic parameter sets (permeability, transfer coefficient, maximum river infiltration flows) as well as two estimations of the water budget from the SIM land surface scheme, and the MODCOU production functions. Consequently, only the best parameter set for the production functions, obtained through calibration, was considered to run MODCOU.

The calibration method is presented in the following section.

Calibration of the surface water budget model

The parameter calibration process can be problematic, because of non-uniqueness (see for example [START_REF] Beven | Prophecy, reality and uncertainty in distributed hydrological modelling[END_REF][START_REF] Gupta | Toward improved calibration of hydrologic models : Multiple and noncommensurable measures of information[END_REF][START_REF] Refsgaard | Parameterisation, calibration and validation of distributed hydrological models[END_REF]. To avoid this kind of issues during the calibration process, the number of parameters to be adjusted is usually reduced. A manner to do this is to assess values for parameters based on their physical significance in the hydrosystem. However the spatial distribution of these physical parameters is often not well known, and the point measurement values may not be representative of grid cell scale parameters. Moreover, model parameters are not always easily linked to physical properties of the system. This is the case for the production function parameters of the MODCOU surface model (Figure 4), which allow partitioning the incoming rainfall between actual evapotranspiration and effective rainfall, and of effective rainfall between infiltration and surface run-off. The values of the parameters are not connected to observed variables and are most often calibrated. Calibration of the water budget is based on the reduction of errors between observed and simulated rivers flows. River gauges out of the Rhine alluvial plain were used to calibrate the production functions, considering that the contribution of aquifer-rivers interactions at these gauging stations is negligible.

This calibration was carried out in a two steps process, with first calibration of the water volumes flowing at the gauging stations, and then of river flow dynamics.

Ten gauging stations in the Vosges Mountains and in the Sundgau area were selected, based on the repartition of the main production functions on their drainage areas. Twenty one production functions were defined on this area, crossing information on land use (Corine Land Cover) and geological soil types (from geological maps). The parameters were calibrated for the fourteen production functions representing more than 20 km² on the gauging stations watersheds, as the influence of the other production functions was not sufficient over these areas. The parameters from the ten remaining production functions were set according to similar soil and vegetation types in previous applications (Golaz et al.2001[START_REF] Gomez | Un outil de modélisation intégrée du transfert des nitrates sur un système hydrologique : application au bassin de la Seine (An Integrated Modelling Tool for Nitrates Transport in a Hydrological System: Application to the River Seine Basin)[END_REF][START_REF] Korkmaz | Application of the coupled model to the Somme river basin[END_REF].

Calibration was performed in a semi automatic way, by testing several predefined values for each parameter. The calibration period was from August 2000 to December 2003, as this period presents contrasted climatic conditions, and has almost complete daily river flow measurements for the selected gauging stations.

For each of the fourteen production functions, five reservoir parameters were adjusted. First the parameters determining the total volume of water flowing at a gauging station (CRT and DCRT) were calibrated based on the relative difference between simulated and observed water volumes over the calibration period following equation 2.

stations n i i obs i sim i obs n V V V stations 1 1 , , ,     , ( 2
)
n stations being the number of gauging stations used for calibration, V sim, i the simulated, and V obs, i the observed volumes of water flowing at a station i.

Then the parameters affecting only the dynamics of river flow (FN, CQI and CQR) were calibrated by optimizing the Nash efficiency (see § 5.1).

Only the best set of parameters was kept, as the sensitivity of the model to groundwater recharge is assessed by using the water budget from the SIM model, as explained below.

Assessment of the hydrogeological simulations

In order to address the uncertainty in the aquifer water budgets and specifically recharge components, several simulations using different equiprobable parameter sets were performed.

Some well constrained parameters were fixed (e.g. calibrated production functions parameters, specific storage), while unconstrained hydrodynamic parameters that influence 1/ the transfer dynamic, 2/ the computation of rivers-aquifer interactions, and 3/ the aquifer recharge were investigated.

1/ Transmissivity: Several distributions were obtained with HPP INV, by using various initial conditions or criteria. These distributions are in good agreement with experimental data from pumping tests [START_REF] Majdalani | Identification of Groundwater Parameters Using an Adaptative Multiscale Method. Application to the Upper Rhine aquifer[END_REF]. Four contrasted distributions were tested in the present study to assess the sensitivity of the model to this parameter (Figure 5 and Table 1).

Two sets have similar mean values, but contrasted standard deviations, and the two others have higher and lower mean values.

2/ To investigate the uncertainty on the river aquifer interactions, the sensitivity of the model to three parameters was tested:

a) The maximum river infiltration rate to the aquifer, for river cells of 200 m: Q lim = 0 (i.e. no infiltration), 25, 50 and 100 L/s (see Equation 1). This parameter influences surface / groundwater interactions, exclusively in the direction of river infiltration towards the water table.

b) The transfer coefficients between the surface and the groundwater parts of the model: the initial value of 0.05 m²/s was multiplied by 2 and by 10. This parameter has an influence over the surface-aquifer interactions in both directions.

c) The Rhine River and the Hardt canal heads: the impact of the variation with time of these river heads on the aquifer and streamflows behaviours was assessed. Although the Hardt canal has far smaller discharges than the Rhine River, preliminary results and KLT analysis have

shown that it was important to represent this irrigation waterway. For this canal it was also tested not taking into account variable water flows.

3/ Uncertainty in the surface water balance is taken into account by using two estimates of the infiltration rates: one from MODCOU (conceptual model), and the other from SIM (ISBA land surface scheme). The latter leads to an infiltration reduced by about 23 %. The specific storage distribution used is the same as in [START_REF] Majdalani | Identification of Groundwater Parameters Using an Adaptative Multiscale Method. Application to the Upper Rhine aquifer[END_REF], with values of 0.05 everywhere, except along a band a few kilometres wide around the southern part of the Rhine River, where the value is 0.12.

A reference simulation, based on mean range parameters (Table 1), is considered in order to assess the impact of different parameters variation independently: eleven other simulations were run by modifying a single parameter values at a time.

Standard analysis

Results obtained with the different parameterizations are quantitatively evaluated in terms of river flows using Nash criterion, and of piezometric levels using biases.

Comparison of the riverflows

Nash criterion [START_REF] Nash | River flow forecasting through conceptual models, 1 : A discussion of principles[END_REF]) is a measure of fit between observed and computed values (Equation 3). Its optimal and maximal value is 1 and its minimum value is -∞. A zero value indicates that the model does not perform better than the mean value of the observations.

            n i obs i obs n i i sim i obs Q Q Q Q Nash 1 2 , 1 2 , , 1 , (3)
with n the number of observed values, Q obs, i the observed river flows, Q sim, i the simulated river flows, with i denoting the time step, and Q obs the arithmetic mean of the observed river flows.

The distribution of Nash criteria for the twelve simulations is shown at each gauging station (Figure 6). General agreement between simulated and observed river flows is found on most mountainous catchments, although some stations have low criterion values. Main discrepancies are located in the plain, which can be explained by the difficulty to reproduce the complex interactions between rivers and the aquifer in this area, and also by an important regulation of discharges in some rivers. The poor results on the German side can be linked to the poorer reconstitution of the atmospheric forcing, as less local climatic data were used in the SAFRAN analysis in Germany. There is no variation between the different simulations in the mountainous catchments as the different parameters tested have no influence on these catchments. The criterion presents some variability in the plain, especially on the Ill River. However, the parameters have a more pronounced impact on low flows. Indeed, the 5-year return period low flow (QMNA5) can almost double from a simulation to another, especially in the southern part of the plain (not shown).

Piezometric levels comparison

Piezometric levels biases are computed as the average difference between simulated and observed piezometric heads.

    n i i obs i sim H H n Bias 1 , , ) ( 1 , ( 4
)
Here n is the number of observed values, H sim, i the simulated piezometric levels, H obs, i the observed piezometric levels, i denoting the time step. Figure 7 shows the spatial distribution of the piezometric levels biases for the reference simulation as well as the mean bias and its standard deviation for the 12 simulations in meter.

The piezometric levels are underestimated in the South Western part and rather well estimated in the centre of the plain. The sensitivity of the parameters tested is higher in the Southern part of the aquifer, with large standard deviations. Indeed in this area piezometric levels are strongly influenced by river infiltration. Some simulations better represent the observed piezometric levels in the South, as illustrated with the reference simulation results. Compared to the reference simulation, a maximum infiltration rate set to Q lim =100 L/s leads to an overestimation of the piezometric levels, with more observation wells having bias values ranging from 1 to 5 m, while the reduced value of this maximum infiltration rate (Q lim = 25 L/s) leads to slightly less values close to 0. The modification of the transfer coefficient for surface/groundwater exchange (T p ) leads to a greater number of negative biases values, but mainly in the range between -2.5 and -1 m.

It is rather difficult to distinguish the results obtained by the four transmissivity distributions, and by the two infiltration rates.

Not taking into account the variations of the Rhine water levels leads to more overestimated piezometric levels, whereas steady values for the Hardt Canal water levels and flow lead to more underestimated levels.

Similar results were obtained for the root mean square error (not shown) for which the maximum dispersion is reached for the Q lim = 0 simulation.

From these comparisons, it is clear that classic statistical results are not sufficient to discriminate the various sets of parameter and that all the simulations except the one with no river infiltration give rather comparable results. In order to go further in the analysis of the simulations, we used a more advanced tool, the Karhunen Loève Transform.

Analysis using the Karhunen-Loève transform method

Principle of the piezometric data analysis

Longuevergne et al. ( 2007) have performed a statistical analysis (Karhunen-Loève transform, hereafter KLT) to separate the main influent processes on piezometric variability from observed time series. KLT, similar to Principal Components Analysis [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences[END_REF], involves summing up a set of possibly correlated time series into a few uncorrelated variables called eigenvectors k Y , related to the original variables by an orthogonal transformation maximizing explained variance. Piezometric head time series, expressed as the variation around the mean level, may then be simply reconstructed as a sum of the temporal behaviours (eigenvectors

) (t Y k ) with associated spatial weight obs k a :    k k obs k obs t Y a t H ) ( ) (
This method was already applied on the Upper Rhine aquifer piezometric observations [START_REF] Longuevergne | Extracting coherent regionel information from local measurements with Karhunen-Loève transform : case study of an alluvial aquifer (Rhine Valley, France and Germany)[END_REF]. Several eigenvectors describing the global temporal behaviour of the aquifer were extracted and interpreted as different contributions to groundwater variations, namely mountain rivers interactions, effective rainfall and Rhine river interactions for the three more important ones. The analysis of the weights a k obs at 190 observation wells allows describing the spatial distribution of these different influences. The results showed that nearly 70% of the aquifer behaviour is determined by interactions with rivers, with a distinction between the Rhine River and the other ones.

In the present study, the Karhunen Loève method is also applied to each simulation in order to assess how well it reproduces the processes involved in the observed piezometric time series.

For this purpose, the model results were projected on the eigenvectors k Y obtained from observations to compute spatial information sim k a for each model and compare them with spatial information from observations, for each separate process described by the eigenvectors, allowing a multi-criterion analysis of the model quality.

Decomposition of the simulated piezometric time series

The five main eigenvectors k Y (interpreted as specific processes) determined by KLT of the piezometric data were able to explain more than 80% of the aquifer behaviour [START_REF] Longuevergne | Extracting coherent regionel information from local measurements with Karhunen-Loève transform : case study of an alluvial aquifer (Rhine Valley, France and Germany)[END_REF].

         5 1 2 2 2 ) ( ) ( ) ( k k sim k obs k sim obs t Y a a t H t H
, where  is a truncation error (only the five first eigenvectors are considered here). Minimization of the root mean square error may be carried out on the spatial information, for each eigenvector associated with highlighted processes.

The first eigenvector (Figure 9), related to the contribution of mountain rivers, explains nearly 50% of the aquifer behaviour. Specifically, variations in annual river discharge at the entrance of the plain where the aquifer lies are closely related to long-term signature of the first eigenvector (Figure 9 a). Spatially, the contribution of the first eigenvector to observations is more important on the mountain piedmont than along the Rhine River (Figure 9 b). The same map can be drawn for the reference simulation, considering the same eigenvector (Figure 9 c), and the difference between the two maps highlights the model error (Figure 9 d). Globally, the agreement is quite good and differences are localized in specific areas. For example, the impact of mountainous rivers seems to be underestimated on the mountain piedmont north of Colmar only. On contrary, the impact is overestimated along the Rhine River and in the areas which are known to be highly anthropized (canals north of Basel, Mulhouse area affected by potash mines and artificial lakes south of Strasbourg).

The method allows pointing out areas where the contribution of each process to the head variability is not properly reproduced by the model. This spatial information may also be interpreted globally as statistical information to assess the ability of the different simulations to capture the global behaviour of the aquifer system (Figure 10). Ideally, the distribution of

obs k sim k a a 
for the 190 observation wells and for a given simulation should be centred on zero with reduced dispersion, for each eigenvector.

For example, considering no river infiltration drives to significant underestimation of the impact of the first eigenvector, corresponding to the contribution of mountainous rivers.

Figure 10 shows that the maximum river infiltration flow of 100 L/s, gives larger errors, and thus, that this parameter set is less probable than others.

The largest contribution to RMSE for all simulations is concentrated on the second eigenvector (Figure 10) previously interpreted as the contribution of effective rainfall. The examination of spatial errors highlights strong differences in the Southern and Eastern parts of the domain. However, the error is reduced when taking into account even partially impact of irrigation in the Hardt area [START_REF] Noyer | Modelling aquifer salinity in the Potash Basin (Alsace)[END_REF], since the error is larger for the simulation without varying water flows and levels in the Hardt irrigation canal (Figure 10) Analysis of the differences on the third eigenvector, related to the impact of the Rhine River, clearly indicates that both infiltration and Rhine river level variations are important processes to describe the interactions between the Rhine River and the aquifer. It also underlines that varying the transfer coefficient values (T P ) can be interesting for a better representation of the Rhine river influence over the aquifer. The four transmissivity distributions tested perform equally well to describe the impact of the mountainous rivers. Hpp103 shows slightly larger underestimation North of Colmar, but also smaller overestimation in Strasbourg area. In the meantime this distribution, which has the lowest mean value, generates less error on the third eigenvector, which has more energy along the Rhine River, and more error on the second eigenvector, representing the rainfall influence.

The statistical analysis of the simulations has shown that one of the twelve simulations, "Q lim =0", failed to reproduce the behaviour of the basin, and that two others, "Steady H&Q Hardt" and "Q lim = 100 L/s", are less probable although they cannot be completely rejected.

The nine remaining simulations cannot really be distinguished statistically, and, although there are still some errors in the representation of the hydrosystem, these nine simulations represent equally well the basin functioning.

Analysis of the water fluxes

Surface-aquifer exchange flux

Surface-aquifer exchange is known to be one key component of the aquifer water budget, although only few quantified estimations based on observations are available.

Figure 11 shows the average and the standard deviation of the simulated mean exchange flows on the period for the nine acceptable simulations. This map highlights the main drainage and river infiltration zones of the hydrosystem. Simulated river infiltration flows (blue cells on Figure 11) are high in the Southern part of the aquifer, mainly along the Rhine River. They are also important in Strasbourg area. These correspond to zones where the infiltration of rivers is known to be important, especially in the Southern part of the plain where several studies have given quantitative estimations of the recharge from the rivers to the aquifer [START_REF] Esteves | Etude et modelisation des relations aquifere-riviere dans le Ried de Colmar (Haut-Rhin, France) (Study and modelling of the aquifer-rivers relations in the Colmar Ried)[END_REF][START_REF] George | Carte hydrogéologique Suisse a 1/100000 -Feuille Basel/Bale (The hydrogeologic map of Switzerland -Basel leaf)[END_REF].

In most part of the domain, the surface-aquifer exchange present low variability (standard deviation values close to 0) except in some parts near the centre of the aquifer (Figure 11). Water table drainage zones (green to red cells on the left and middle maps of Figure 11) occur in the Southern part of the aquifer, near its borders and along the Rhine River. Drainage zones extensions are quite important in the centre of the plain between Colmar in the South and Strasbourg in the North. This corresponds to the Ried area, an important wetland with shallow water table, leading to an important aquifer overflow and the presence of groundwater-fed streams.

In the Northern part of the aquifer, important drainage flows from the aquifer occur, due to the proximity of the water table to the topographic surface.

Comparisons were made with values found in the literature [START_REF] Piren-Eau-Alsace | Recherche methodologique sur les hydrosystemes pour optimiser la gestion des ressources en eau dans la region Alsace (Methodologic research on hydrosystems to optimize water resources management in the Alsace region[END_REF][START_REF] Esteves | Etude et modelisation des relations aquifere-riviere dans le Ried de Colmar (Haut-Rhin, France) (Study and modelling of the aquifer-rivers relations in the Colmar Ried)[END_REF][START_REF] George | Carte hydrogéologique Suisse a 1/100000 -Feuille Basel/Bale (The hydrogeologic map of Switzerland -Basel leaf)[END_REF][START_REF] Schmitt | Typologie hydro-geomorpholoqique fonctionnelle de cours d'eau : recherche methodologique appliquee aux systemes fluviaux d'Alsace (Hydrogeomorphologic functionnal typology of streams : methodologic research applied to the fluvial systems in Alsace)[END_REF]. These values often correspond to limited space and time extents, but it is interesting to know whether the orders of magnitude of the simulated exchange are plausible. [START_REF] Esteves | Etude et modelisation des relations aquifere-riviere dans le Ried de Colmar (Haut-Rhin, France) (Study and modelling of the aquifer-rivers relations in the Colmar Ried)[END_REF] indicates an average infiltration flow of 200 to 300 L/s on the Fecht River between its entry into the plain and its confluence with the Weiss River. On this reach an infiltration flow of 1000 L/s is simulated on average. Downstream the drainage flow varied from 100 to 300 L/s on average [START_REF] Esteves | Etude et modelisation des relations aquifere-riviere dans le Ried de Colmar (Haut-Rhin, France) (Study and modelling of the aquifer-rivers relations in the Colmar Ried)[END_REF], whereas MODCOU results give an average of 1635 L/s. Thus for this river the intensity of the exchange is certainly overestimated by the model.

Along the Ill River, [START_REF] Esteves | Etude et modelisation des relations aquifere-riviere dans le Ried de Colmar (Haut-Rhin, France) (Study and modelling of the aquifer-rivers relations in the Colmar Ried)[END_REF] gives a value of 200 L/s/km infiltrating from the river to the aquifer between Ensisheim and Colmar. The [START_REF] Piren-Eau-Alsace | Recherche methodologique sur les hydrosystemes pour optimiser la gestion des ressources en eau dans la region Alsace (Methodologic research on hydrosystems to optimize water resources management in the Alsace region[END_REF] gives a value of 130 L/s/km infiltrating between Ensisheim and Oberhergheim, and 37 L/s/km between Oberhergheim and Horbourg-Wihr, for the period 1971-1974[START_REF] George | Carte hydrogéologique Suisse a 1/100000 -Feuille Basel/Bale (The hydrogeologic map of Switzerland -Basel leaf)[END_REF] estimate the infiltration between Mulhouse and Ensisheim amounts to 65 L/s/km in average.

The simulated exchange flows give an average infiltration of 124 L/s/km between Mulhouse and Ensisheim, 93 L/s/km between Ensisheim and Oberhergheim, 87 L/s/km between Oberhergheim and Horbourg-Wihr and 79 L/s/km between Ensisheim and Colmar.

Thus simulated infiltration of the Ill River between Mulhouse and Colmar is less variable than estimated by observations, and it is overestimated upstream of Ensisheim and downstream of Oberhergheim, while it seems underestimated between Ensisheim and Oberhergheim, i.e. on the section just downstream of the confluence with the Thur River. However most of these estimates are made by punctual gauging, a method lacking precision, and giving results with time limited representativeness. Similar results are found on the Thur and Doller rivers, at the South western end of the basin [START_REF] George | Carte hydrogéologique Suisse a 1/100000 -Feuille Basel/Bale (The hydrogeologic map of Switzerland -Basel leaf)[END_REF][START_REF] Schmitt | Typologie hydro-geomorpholoqique fonctionnelle de cours d'eau : recherche methodologique appliquee aux systemes fluviaux d'Alsace (Hydrogeomorphologic functionnal typology of streams : methodologic research applied to the fluvial systems in Alsace)[END_REF].

It comes to the conclusion that the MODCOU model with the selected parameter sets is able to represent the main spatial pattern of the aquifer-surface interactions. However, it has some difficulties capturing the spatial variability of the exchange even by using several parameter sets. On average this flux is overestimated in some sections of the hydrographic network. This is probably due to the use of spatially homogeneous distribution for the parameters T p and Q lim directly influencing the river-aquifer exchange.

Water Budgets

In order to address uncertainties in the functioning of the Upper Rhine aquifer, an analysis of the water budgets computed by the different simulations is performed. +/-37 m 3 /s. Therefore, the KLT analysis allowed reducing the uncertainty from 35% to 14%.

These results are consistent with those obtained in the framework of the LIFE project. A value of 85 % was found with a steady state simulation for a low water situation [START_REF] Lfu | Action de démonstration portant sur la protection et la gestion des réserves en eau souterraine dans la partie franco-germano-suisse de la vallée du Rhin Supérieur -Rapport final (Demonstration action on protection and management of groundwater resources in the french-german-swiss part of the Upper Rhine valley -Final report[END_REF]. On the contrary in the Monit project the main contribution to the recharge of the aquifer was found to be the effective rainfall infiltration over its surface [START_REF] Lubw | Projet INTERREG III Monit : Modelisation hydrodynamique et transport des nitrates (INTERREG III project Monit : Hydrodynamic modelling and nitrates transport)[END_REF].

The lateral subsurface flow is about two times lower than the effective rainfall infiltration, and represents approximately 7% of the water budget.

The main output process in all cases is the drainage of the aquifer by the rivers and by overflow.

Another interesting variable is the aquifer storage variations. To analyse the uncertainty on this variable according to the parameter sets, the evolution of the aquifer storage were converted into water height variations over the entire surface of the aquifer. Thus, a variation of 1 mm over the whole aquifer represents an evolution of the aquifer volume of 4,655,000 m 3 .

These variations are computed by reference to the aquifer stock at a given date, namely the 30 th July 1986. This date was chosen in order not to take into account variations during the first year of simulations, which could be due to inaccurate estimates of initial piezometric levels. The stock variations are not very sensitive to the transmissivity distributions although the differences between the four distributions can reach 10 mm at some periods (Figure 13).

The different transfer coefficient values (T p ) lead to quite similar stock variations, even for the higher value of T p =0.5 m 2 /s, which was shown to be the case having the highest river infiltration flow. In this case, the stock is increased by only a few millimetres. Indeed there is compensation between the river infiltration and drainage rates leading to a similar evolution of the aquifer storage.

The reduced rainfall infiltration computed by SIM also has a limited impact on the aquifer storage variations. Indeed the maximum impact reaches 10 mm at most between the reference simulation and the reduced infiltration simulation.

To summarize, the aquifer stock variations are rather sensitive to the maximum river infiltration rate, whereas they are less sensitive to rainfall infiltration and to transfer coefficients. This can be explained by the fact that a modification of the infiltration over the aquifer surface is balanced by a modification of the exchange with rivers flowing from the mountains. Increasing the transfer coefficients leads to larger exchange between the aquifer and the surface in both directions. This phenomenon can explain why the aquifer stock is not really affected. On the other hand when the maximum river infiltration rate is modified, it only affects the exchange from the rivers to the aquifer, leading to increased aquifer stock variations.

Conclusion

The objective of the study was to increase understanding of the Upper Rhine Graben aquifer Although the developed model is not suited to represent fine scale processes and fully coupled interactions between surface water and groundwater, it represents the main processes that impact the aquifer water budget, and thus, allows an analysis of their spatial distribution.

Twelve model runs were performed using several parameter sets for the transmissivity, the surface/groundwater transfer coefficient, the maximum river infiltration flow, Hardt Canal water levels and flows and Rhine river water levels, and two estimations of the surface water budget in the plain.

These twelve runs were assessed using classic statistical criteria: the daily Nash efficiency for the river flows, and biases for the piezometric levels. These analyses show that for each run, the model had some problems reproducing the daily discharges on the German side, probably due to some weaknesses of the atmospheric analysis using fewer meteorological stations in Germany than in France. The biases on the simulated piezometric heads are rather weak except on the southern part of the basin, where the discrepancy between the twelve simulations is large. These classical analyses were able to highlight that one set of parameters was not coherent with the functioning of the aquifer system, and thus that it is required to take into account the possibility for rivers to infiltrate.

In order to go further in the assessment of the simulations results, a mathematical analysis method based on a Karhunen Loève transform (KLT) of piezometric time series was used to analyse the simulated piezometric levels temporal variations, as well as the observed ones.

This method allows separating the different processes influencing the water table evolution, in the form of eigenvectors representing river infiltration and the effective rainfall influences.

The spatial distribution of the projections stemming from the decomposition of piezometric chronicles as linear combinations of the eigenvectors allows having an idea of the area of dominating influence of the three main processes. Comparison of the projections obtained with simulated and observed piezometric time series also highlights the main areas where difficulties arise to reproduce the water table behaviour. KLT showed that all model runs failed to represent the annual component of the second eigenvector, related to the impact of effective rainfall, mainly due to large errors in the Hardt area in the south of the domain, and in the Eastern part of the plain. This error is however reduced when taking into account the irrigation practises in the Hardt area, with irrigation water flowing through channels with high bed permeability and thus large river-aquifer exchanges. The KLT method underlines the need to better take into account these processes in the model. This method also proved an interesting and original technique to assess the model performances in reproducing the spatial and temporal variability of water table levels, as well as the recharge/discharge processes inducing these variations. It might be considered using this method in order to test spatial patterns of some parameters taken homogeneous for now, such as the transfer coefficient T p and the maximum river infiltration rate Q lim .

Comparison of the aquifer water budgets obtained with different parameter sets shows that there is a balance between the magnitude of river to aquifer infiltration and drainage of the aquifer so that the aquifer stock is rather stable for each run.

The intensity of the simulated river-aquifer exchange was compared to available observations, mainly on the Ill River. Estimated values are in the same order as the observed ones. However the model presents smoother spatial variations of the simulated fluxes than observed, which could be due to the use of spatially homogeneous river bed parameters.

The KLT method has shown that nine over twelve simulations were able to have a similarly satisfying representation of the aquifer behaviour. This allows reducing the uncertainty of aquifer recharge by river infiltration from 35 % to 14 % with an average flux of 115 m 3 /s.

This represents about 80% of the aquifer inflow.

This proportion compares well to the LIFE project estimate for a low water situation, but is more important than the one from the Monit project for a transient state simulation.

This study brings increased comprehension of the functioning of an alluvial aquifer, which is strongly influenced by groundwater /surface water interactions. It shows that modelling, together with inversion and advanced statistical methods are powerful tools to quantify the water budget of an alluvial aquifer.

Figure 1 :

 1 Figure 1 : a) Situation of the basin and topography with main rivers and gauging stations used for model calibration; b) Mean annual precipitation in mm/yr computed on the period 1983-2006 as estimated by the SAFRAN analysis (cf section 2.2).

Figure 2 :

 2 Figure 2 : Geological cross section of the Rhine Graben, showing the sediment succession in the plain and relative thickness (modified from Illies (1972))

Figure 3 :

 3 Figure 3 : a) Discretization of the surface (light red) and the aquifer parts (green) of the model. The border cells where the surface runoff is infiltrating are plotted in dark brown ; b) Land use data extracted from Corine Land Cover 2000 data set.

Figure 4 :

 4 Figure 4 : schematization of a production function in MODCOU for surface water budget computation[START_REF] Ledoux | Spatially distributed modeling: conceptual approach, coupling surface water and groundwater[END_REF]

Figure 5 :

 5 Figure 5 : Two of the four transmissivity distributions tested, obtained through inversion with the HPP INV model.

Figure 6 :

 6 Figure 6 : Spatial distribution of the daily Nash criterion values. At each gauging station a pie is plotted representing the distribution of the Nash criterion values for the 12 simulations.

Figure 7 :

 7 Figure 7 : Comparison between observed and simulated piezometric levels at the 190 observation wells over the 1986-2002 period. Left: bias [m] computed for the reference simulation, middle: mean bias [m] for the 12 simulations, right: standard deviation of the bias over the 12 simulations [m].

Figure 8 :

 8 Figure 8 : Bias histogram computed on the 190 observation wells, for the twelve simulations

  are compared to those obtained by projection of the model results on the eigenvectors k Y , sim k a , for each simulation. As eigenvectors are orthogonal, the root mean square error (RMSE) between model and observations may be simply written as

Figure 9 :

 9 Figure 9: Spatial and temporal analysis with the KLT method for the first eigenvector, associated with the influence of mountain rivers over groundwater variability. a) First eigenvector time series; b) and c) Maps of the weights associated with the first eigenvector for the observations and reference simulation respectively; d) Map of the differences between simulated and observed weights.

Figure 10 :

 10 Figure 10: Errors in the reproduction of the eigenvectors weights at the 190 observation wells for the twelve simulations a) Mean total errors for the five eigenvectors and standard deviations ; b) Contribution to the total error of eigenvectors 1, 2 and 3.

Figure 11 :

 11 Figure 11: Mean exchange fluxes between the surface and the aquifer for the reference simulation (left) and on average on the 9 simulations (centre), and its relative standard deviation (right)

Figure 12 :

 12 Figure 12: Mean annual partition of the aquifer inflows between lateral subsurface flow, effective rainfall and river infiltration (right bars), and of the aquifer outflows between pumping and rivers drainage (left bars).

Figure 12

 12 Figure12shows the water budget for the aquifer recharge and discharge for the nine plausible

Figure 13 :

 13 Figure 13: Monthly aquifer stock evolution from 1986 to 2003 as estimated by the 9 acceptable simulations

  functioning. Modelling of the Upper Rhine graben aquifer was performed with the MODCOU coupled hydrogeological model computing the surface water budget, groundwater flow and interactions between the aquifer and the surface. Distributions of aquifer transmissivities were provided by the inversion algorithm of the HPP-INV hydrogeological model. In the present study, the sensitivity of the MODCOU model results to some model parameters governing the surface / groundwater interactions was assessed, specifically in terms of aquifer behaviour.

Future

  Figure1

Table 1

 1

	summarizes the different

Table 1 : Recapitulation of the different model parameterizations tested

 1

Acknowledgements

This work was supported by the program VMC-2007 of the French Agence Nationale de la Recherche (ANR) under the project VULNAR.

We would like to thank the different persons and institutions who provided us with data for the modelling, especially: Michel Wingering at the LUBW for German hydrometric data; Dr.

Werner Weinzierl and Dr Frank Waldmann at the LGRB for the German soil type data; Philippe Elsass at the BRGM for its help in the understanding of the Upper Rhine hydrosystem; Jean Lett and Claude Husser at the Navigation Service of Strasbourg for the Rhine River levels data; Joëlle Sauter at the ARAA for the French soil data.


Jean-Charles Billaut
email: jean-charles.billaut@univ-tours.fr

Emmanuel Hebrard
email: hebrard@laas.fr

Pierre Lopez
email: lopez@laas.fr

Complete Characterization of Near-Optimal Sequences for the Two-Machine Flow Shop Scheduling Problem

In a two-machine flow shop scheduling problem, the set ofapproximate sequences (i.e., solutions within a factor 1+ of the optimal) can be mapped to the vertices of a permutation lattice. We introduce two approaches, based on properties derived from the analysis of permutation lattices, for characterizing large sets of near-optimal solutions. In the first approach, we look for a sequence of minimum level in the lattice, since this solution is likely to cover many optimal or nearoptimal solutions. In the second approach, we look for all sequences of minimal level, thus covering all -approximate sequences. Integer linear programming and constraint programming models are first proposed to solve the former problem. For the latter problem, a direct exploration of the lattice, traversing it by a simple tree search procedure, is proposed. Computational experiments are given to evaluate these methods and to illustrate the interest and the limits of such approaches.

Introduction

Computing the set of near-optimal solutions of a combinatorial problem has many applications. It is for instance the case when there is some uncertainty in the problem definition. For example, production systems are subject to numerous disturbances deriving mainly from data uncertainty and unexpected events. These disturbances result in lack of raw materials, tardiness, tool failure, etc., and may make the initial planning infeasible. In particular, this motivates research works on flexibility and robustness in scheduling [START_REF] Billaut | Scheduling with Flexibility and Robustness[END_REF]. This is why dynamic or reactive approaches have been developed [START_REF] Billaut | A new method for workshop real time scheduling[END_REF][START_REF] Sabuncuoglu | Analysis of reactive scheduling problems in a job shop environment[END_REF][START_REF] Sabuncuoglu | Reactive scheduling in a dynamic and stochastic fms environment[END_REF]. Among the possible dynamic approaches, some of them are based on the characterization of sets of solutions. In this context, the set of solutions can be used in an interactive and/or dynamic environment to guide the decision making from an optimal solution to another one when an unexpected event occurs, or on a user request.

Another important reason to obtain the characteristics of optimal solutions is to solve multiobjective scheduling problems [START_REF] Gupta | Two-machine flowshop scheduling with a secondary criterion[END_REF]. These characteristics could be incorporated into branch-and-bound procedures, in order to prune nodes more efficiently, hopefully.

In [START_REF] Artigues | Maximization of solution flexibility for robust shop scheduling[END_REF][START_REF] Artigues | Characterization of a set of schedules in a resource-constrained multi-project scheduling problem with multiple modes[END_REF][START_REF] Billaut | A new method for workshop real time scheduling[END_REF], the set of solutions or schedules is given by a sequence of groups of jobs on each machine, where the sequence of jobs inside each group is not fixed (the authors talk about 'groups of permutable jobs'). All the characterized solutions ensure a performance guarantee on a given quality measure. This methodology has been implemented in a software and has been used in practice by several companies [START_REF] Roubellat | Ordonnancement d'ateliers : d'Orabaid à Ordo (in French)[END_REF]. In [START_REF] Aloulou | Worst-case evaluation of flexible solutions in disjunctive scheduling problems[END_REF], a set of semi-active schedules is characterized by a partial order of jobs on each machine. The authors propose a method for computing the best case and the worst case performances. In [START_REF] Briand | A robust approach for the single machine scheduling problem[END_REF], the authors propose a method for characterizing a large set of optimal solutions, based on the analysis of interval structures and on a theorem established in [START_REF] Erschler | A new dominance concept in scheduling n jobs on a single machine with ready times and due dates[END_REF]. The method is applied to a single-machine problem, where jobs have release dates and due dates. All these methods allow the characterization of a subset of the set of optimal or approximate solutions. Another way for this characterization is to provide constraint propagation techniques or dominance properties, in order to maintain the tightest set of remaining consistent decisions [START_REF] Briand | Constraint-based approaches for robust scheduling[END_REF].

It is well known that there potentially exist several optimal solutions to a given scheduling problem. They may even have a huge number of different optimal solutions [START_REF] Billaut | Enumeration of all optimal sequences in the two-machine flowshop[END_REF][START_REF] Sevastyanov | Efficient enumeration of optimal and approximate solutions of the two-machine flow-shop problem[END_REF], for example several hundreds of thousands for mid-size academic instances. Moreover, it is well known as well that it is often not easy to find an optimal solution to a scheduling problem because of its NP-completeness in the general case. However, some scheduling problems can also be solved in polynomial time using specific methods, which generally consist in sorting the jobs according to a simple priority rule.

In this paper, we investigate the possibility to characterize the whole set of optimal solutions of the two-machine flow shop scheduling problem, i.e., to give the analytical characteristics of these solutions. Obviously, the aim is not to enumerate these solutions, but to describe their characteristics. In our approach, this is achieved by the knowledge of dominance rules for the problem under consideration and through the study of the various optimal solutions associated with the vertices of a lattice. Since a realistic goal (for an enterprise for example) is not necessarily to focus on optimal solutions, note that the same methods apply for characterizing the whole set of -approximate solutions, i.e., the set of solutions in the lattice with a performance not worse than a given distance function of from the value of the optimal solution.

The rest of the paper is organized as follows. Section 2 provides the necessary notations, definitions, and properties concerning the main mathematical object used in this work: the "lattice of permutations". Then in Section 3, we recall how all optimal solutions of a class of scheduling problems can be characterized by a subset of vertices of minimal level in the permutohedron. Then we present the problem addressed in this paper, namely: finding a minimum vertex, and finding all minimal vertices. In Section 4 we propose an Integer Linear Programming (ILP) approach as well as a Constraint Programming (CP) approach for the former problem, whereas Section 5 describes an algorithm to solve the latter problem. Section 6 presents the results obtained from computational experiments and Section 7 gives a conclusion and some future research directions.

The Lattice of Permutations: Definitions and Properties

We consider the set {1, 2, . . . , n} of integers and S n the group of all permutations on {1, 2, . . . , n}. We represent the members of S n by strings of integers. As an example, consider n = 4; σ = 4213 denotes a permutation σ where σ(1) = 4, σ(2) = 2, σ(3) = 1, and σ(4) = 3. Using the same notations as in [START_REF] Markowsky | Permutation lattices revised[END_REF], we denote by index(σ, i) the position of integer i in permutation σ. For the previous example, we have index(σ, 1) = 3.

With the elements of S n we define a directed graph where the nodes are the elements of S n . In this digraph, there exists an edge between nodes σ and σ if and only if σ = αijβ with α and β two partial orders, i, j ∈ {1, 2, . . . , n} with index(σ, j) = index(σ, i) + 1, σ = αjiβ, and i < j. In other words, there is an edge between σ and σ if these permutations are the same, except that there exist i and j, two consecutive jobs with i < j, that are in the reverse order in σ .

This digraph is a lattice, called the lattice of permutations or permutohedron [START_REF] Bowman | Permutation polyhedra[END_REF]. Figure 1 gives the lattices of permutations for n = 3 and n = 4.

To each permutation in the lattice can be associated a level. There are at most n(n -1) 2 + 1 levels. By convention we say that permutation (n, n-1, n-2, . . . , 1)

is at level 0 and that permutation (1, 2, 3, . . . , n) is at level n(n-1)/2. We denote by κ(σ) the level of permutation σ. For a given permutation σ, we denote by Γ (σ) the set of couples defined as follows:

Γ (σ) = {(i, j) ∈ {1, 2, . . . , n} 2 | i < j and index(σ, i) < index(σ, j)}
For example, in permutation σ = 4132, we have Γ (σ) = {(1, 3), (1, 2)}.

We now report some properties associated to the lattice of permutations previously defined. Property 1. [START_REF] Bennett | Two families of Newman lattices[END_REF]: For any permutation σ, the level of σ is exactly its number of inversions from permutation at level 0, i.e., the number of times we have i < j and index(σ, i) < index(σ, j): If we consider the elements of Γ (σ) as a set of constraints associated to σ, we can say that all the predecessors of σ have to satisfy at least the same constraints as σ. We then claim that Γ (σ) gives the characteristics of all the predecessors of permutation σ in the digraph.

κ(σ) = |Γ (σ)| Property 2.

Characteristic Solutions

In many scheduling problems, a set of dominant solutions (e.g., left-shifted schedules) can be mapped with the set of permutations of jobs and consequently to the vertices of a permutohedron. Moreover, in a previous paper [START_REF] Billaut | Characterization of all rho-approximated sequences for some scheduling problems[END_REF], we showed that, for a number of classes of such problems, by carefully choosing the labeling of the jobs, one can design a permutohedron with a very desirable property: none of its children can have a strictly better objective value than its own. Thus, the set of optimal solutions with minimal level in the permutohedron is sufficient to characterize all optimal solutions. Indeed, a solution is optimal if and only if it belongs to this set, or is an ancestor of an element in this set.

The problems for which this method holds are those where a simple rule (checkable in O(1)) exists guaranteeing that given two consecutive jobs in a sequence, it is possible to say which job it is always better to put first. Such a rule allows the building of optimal algorithms in O(n log n) time complexity, after sorting the jobs according to this rule. They are therefore all easy to solve.

However, finding all solutions, or a large number of solutions to such problems may still be difficult.

In terms of scheduling, a permutation is a sequence of jobs. Let f (σ) denote the objective value (to minimize) of a sequence σ. We say that an ordering relation R is an improving relation iff for any sequence αjiβ (with α and β two arbitrary subsequences and i and j two consecutive jobs) we have:

iRj ⇒ f (αijβ) ≤ f (αjiβ)
This property is often referred to as "pairwise interchange" mechanism, and it is frequently used for proving properties of scheduling algorithms. Indeed it is easy to see that the sequence σ built by sorting all jobs according to an improving order relation R is optimal. Starting from an arbitrary sequence, one can obtain σ by iteratively swapping consecutive elements while never degrading the objective value.

Assume that sequence (1, 2, . . . , n) is an optimal sequence given by a sorting algorithm based on an improving relation R and consider the permutohedron rooted in (1, 2, . . . , n). By definition, given any sequence σ and one of its successor σ , we have f (σ) ≤ f (σ ). Therefore, if a sequence σ is optimal, all the predecessors of σ are optimal sequences. Moreover, we say that a sequence σ is minimal with respect to a permutohedron, iff it is optimal and none of its children are, i.e., any further swap of two consecutive jobs leads to a sub-optimal solution. Notice that the method that we introduce in this paper applies to any problem for which an improving relation exists (this is true for many scheduling problems).

In this paper, we are concerned with finding all minimal sequences in such permutohedra and, in particular, the minimal sequence with minimum level in the lattice (so-called minimum sequence). Intuitively, the optimal sequence with minimum level covers many optimal solutions. It can therefore be seen as robust in the sense that, when for some reason this solution is no longer valid, one can easily generate optimal solutions by swapping jobs in order to climb up the lattice, with the hope that at least one of these solutions will still be valid. Another possibility is to compute all minimal sequences. Indeed, this is a relatively concise way of storing all the optimal solutions of the problem.

We shall introduce several methods to address the two following problems:

-PB1, finding an optimal sequence with minimum level in the lattice; -PB2, finding all minimal sequences in the lattice.

Example 1. For instance, suppose that the sequence (1234) in Figure 1 is given by a sorting algorithm with respect to an improving relation R and that the sequences (3142) and (2431) are the only minimal sequences (they are underlined in Figure 1). Then we know that the set of optimal solutions is exactly the set composed of themselves and of their ancestors: {(1234), (2134), (1324), ( 1243), (2314), ( 2143), (3124), (1342), (2341), ( 2413), (3142), (2431)} (marked in italic).

Observe that the minimum sequence is (2431) and covers 8 optimal solutions (including itself), whereas (3142) covers only 5.

We consider a set of n jobs to schedule in a two-machine flow shop, where the two machines are denoted by M 1 and M 2 and we denote by p j,1 and p j,2 , the processing time of job j on machine M 1 and on machine M 2 , respectively, 1 ≤ j ≤ n. C j denotes the completion time of job j on machine M 2 . The criterion under study is the makespan C max = max 1≤j≤n C j . The problem that we consider is denoted by F 2||C max and is solved to optimality by ordering the jobs following Johnson's dominance condition [START_REF] Johnson | Optimal two-and three-stage production with setup times included[END_REF].

Considering the makespan minimization (min C max in our case), we define an -approximate sequence ( ∈ R + ) as a sequence σ such that:

C * max ≤ C max (σ) ≤ (1 + ) × C * max
In the following expressions, it is sufficient to replace C * max by (1 + )C * max for characterizing the -approximate solutions.

Finding the Minimum Sequence

In order to find a schedule of minimum level in the permutohedron we solve the scheduling instance with some slight modifications. We assume that the root sequence σ * of the lattice was obtained using Johnson's algorithm, and such that jobs are renumbered with respect to this first optimal sequence. Then we re-solve the initial instance with the following changes:

-The objective function C max to minimize is replaced by the constraint stating that the expression C max should be less than (1 + )C max (σ * ). -We use a new objective function: minimize the level in the lattice.

The complexity of the modified problem is not known. However, there are examples of tractable problems that become NP-hard when adding a similar objective function. For instance monotone-2SAT (the satisfiability of a 2-CNF formula with only positive literals) is NP-hard if the number of atoms set to true is to be minimized [START_REF] Bailleux | DISTANCE-SAT: Complexity and Algorithms[END_REF].

Integer linear programming approaches

We first propose an integer linear programming model with the variables defined as follows: for all 1 ≤ i < j ≤ n, y i,j is a binary variable (0-1) equal to 0 if job i precedes job j in the sequence and 1 otherwise. We also introduce continuous variables: t j,1 and t j,2 are the start times of job j on machine M 1 and on machine M 2 , respectively. HV stands for High Value, and can be set for example to n j=1 2 k=1 p j,k . In addition, it is assumed that the binary variables verify a kind of triangle inequality: y i,k ≤ y i,j + y j,k , ∀i, j, k ∈ {1, . . . , n}, i = j = k

We have:

t j,k ≥ t i,k + p i,k -HV y i,j , (4.1)
t i,k ≥ t j,k + p j,k -HV (1 -y i,j ), (4.2)
(4.1) and (4.2), ∀k ∈ {1, 2}, ∀i, j ∈ {1, . . . , n}, i = j t j,2 ≥ t j,1 + p j,1 , ∀j ∈ {1, . . . , n}

To take account of the objective function, let C * max denote the optimal value -supposed to be known -of the makespan.

We post the following constraints:

t j,2 + p j,2 ≤ C * max , ∀j ∈ {1, 2, . . . , n} (4.3)
Another way to express Property 1 is to say that the level in the lattice is equal to the cardinality of the set of permutations preserved from the root sequence σ * . Notice that the characteristic function of the set Γ (σ) is given by the set of variables y i,j (with the correspondence y i,j = 0 ⇔ (i, j) ∈ Γ (σ)). Finding the sequence with minimum level can therefore be expressed by the following objective function:

MAX n i=1 n j=i+1 y i,j

Constraint programming approach

We propose a constraint programming model similar to the integer linear programming model described in Section 4.1, however solved using a different approach.

As in the ILP model, we introduce a binary variable y i,j for each pair of jobs, taking the value 0 if job i precedes job j and 1 otherwise. We also introduce integer variables t j,1 and t j,2 for the start times of job j on machine M 1 and on machine M 2 , respectively.

We post precedence constraints between the two activities of each job as follows:

y i,j ⇒ t j,k ≥ t i,k + p i,k , (4.4
)

¬y i,j ⇒ t i,k ≥ t j,k + p j,k , (4.5)
(4.4) and (4.5), ∀k ∈ {1, 2}, ∀i, j ∈ {1, . . . , n}, i = j t j,2 ≥ t j,1 + p j,1 , ∀j ∈ {1, . . . , n} Note that we do not post constraints on triplets of binary variables to exclude cycles. Bound consistency is enforced on these constraints, i.e., when the value of the binary variable y i,j is set, the bounds of the integer variables t i,k and t j,k are updated with respect to the constraint t i,k ≥ t j,k + p j,k , or t j,k ≥ t i,k + p i,k according to y i,j 's value. Conversely, if one of the two precedence constraints becomes entailed or disentailed because of the current domain of t i,k and t j,k , then the value of y i,j is set accordingly.

Exactly as in the ILP model, we set a maximum value to the objective function using the same constraints. Moreover, finding the sequence with minimum level is expressed by the same objective function:

MAX n i=1 n j=i+1 y i,j
The method that we used to solve this constraint programming model is essentially the same used for several variants of the job-shop scheduling problem in [START_REF] Grimes | Job Shop Scheduling with Setup Times and Maximal Time-Lags: A Simple Constraint Programming Approach[END_REF][START_REF] Grimes | Closing the Open Shop: Contradicting Conventional Wisdom[END_REF].

Search method: The problem is solved using a dichotomic search on the objective, yielding a satisfaction problem at each step. Each of these steps is bounded by the number of nodes of the search tree that can be expended during search. If the dichotomic search is not conclusive, a branch-and-bound procedure is used starting from the bounds computed during the dichotomic phase.

Variable selection heuristic. We use a slightly modified version of the domain over weighted-degree heuristic [START_REF] Boussemart | Boosting Systematic Search by Weighting Constraints[END_REF] to select the next variable to branch on. Observe that we branch only on the binary variables (y i,j ). Let w(t i ) be the number of times search failed while propagating any constraint involving job i, and let min(t i,k ) and max(t i,k ) be, respectively, the minimum and maximum starting time of t i,k at any point during search. The next variable y i,j to branch on is the one minimizing the value of:

(max(t i,k ) + max(t j,k ) -min(t i,k ) -min(t j,k ) + 2)/(w(t i,k ) + w(t j,k ))

Value selection heuristic. When branching on the variable y i,j we try first the value assigned to this variable in the best feasible solution found so far. If no solution has been found, the value 0 is tried first. This idea is a simplified version of the solution guided approach (SGMPCS) proposed by Beck for job-shop scheduling problems [START_REF] Beck | Solution-Guided Multi-Point Constructive Search for Job Shop Scheduling[END_REF].

Restarts with Nogoods. We use a geometric restarting strategy [START_REF] Walsh | Search in a Small World[END_REF]. When a given number of nodes have been explored, we stop the current exploration and start again from scratch. The limit in number of nodes grows geometrically: it is of the form s, sr, sr 2 , sr 3 , . . . where s is the base and r is the multiplicative factor. In our experiments the base was 256 failures and the multiplicative factor was 1.3. Moreover, after each restart, the dead ends of the previous explorations are stored as clausal nogoods [START_REF] Lecoutre | Nogood Recording from Restarts[END_REF].

Finding all Minimal Sequences

Observe that it is possible to solve this problem using iteratively the method described in the previous section, and avoiding rediscovery of previous solutions with nogoods.

Suppose that a first sequence σ 0 of minimum level has been computed, and suppose that it corresponds to the set of precedences index(σ 0 , a 1 ) < index(σ 0 , b 1 ) and index(σ 0 , a 2 ) < index(σ 0 , b 2 ) and ... and index(σ 0 , a ν0 ) < index(σ 0 , b ν0 ). The clause index(σ 0 , b 1 ) < index(σ 0 , a 1 ) or index(σ 0 , b 2 ) < index(σ 0 , a 2 ) or ... or index(σ 0 , b ν0 ) < index(σ 0 , a ν0 ) can be added to the model in order to avoid finding σ 0 again.

For the CP formulation we can add:

y a1,b1 ∨ y a2,b2 ∨ . . . ∨ y aν 0 ,bν 0 (5.1)
And for the ILP formulation:

y a1,b1 + y a2,b2 + . . . + y aν 0 ,bν 0 ≥ 1 (5.2)
At each iteration, only one additional constraint of this type is produced and no additional variable is generated.

However, we shall see that such an approach is not efficient, hence we propose an algorithm that directly explores the lattice in depth first order and finds all minimal sequences. The only difficulty is to avoid exploring twice the nodes of the lattice, given that we cannot store it explicitly. As in the previously discussed models, the method starts from a schedule σ * given by rule R, then the jobs are renumbered with respect to this first optimal solution. However, instead of solving the scheduling problem, we explore the lattice depth first in Algorithm 1. In other words, a move occurs on the lattice by swapping elements of the optimal order only if the move corresponds to a downward edge in the lattice, and if it does not degrade the objective value below (1 + ) × C max (σ * ).

Let σ be a permutation on {1, . . . , n}. The only operations that we use to move on the lattice is to swap two consecutive elements. We denote by swap(σ,(a, b)) the permutation σ equal to σ except that a and b are swapped. For instance, swap((41325),(1, 3)) = (43125). We denote by opt(σ) the fact that the objective value of the schedule σ is within the tolerated interval. The current permutation is kept in the variable σ, whenever we reach a local optimum, i.e., a permutation σ such that any swap increases the objective value above the acceptable threshold, we print it (Line 4).

Algorithm 1 Explore-Lattice Data: σ, Γ (σ), explored minimal ← True; 1 foreach k ∈ {1, . . . , n -1} do a ← σ(k); b ← σ(k + 1);
We keep track of the current value of Γ (σ) using a set. Last, we use another set denoted explored to store the elements of Γ (σ) that we already explored in previous branches, to avoid visiting twice the same vertex of the lattice.

The first loop (in Line 1) goes over all the possible swaps in the sequence σ. The swaps that satisfy the conditions in Line 2 are actually edges in the lattice leading to an -approximate solution. If such an edge exists, the current node is not a local optimum. Then, the condition in Line 3 ensures that the successor has not yet been explored. Example 2. To illustrate the exploration of the lattice, take again the example depicted in Figure 1 for n = 4. Figure 2 now displays the corresponding sets Γ (σ). Algorithm 1 explores all the solid edges of this digraph. For each dotted edge adjacent to a solid edge, we mark the reason why it is not explored with the following notation:

: The child has an objective value strictly greater than (1 + ) × f (σ * ). Γ : The successor edge corresponds to swapping a for b, however (a, b) ∈ explored.

Theorem 1. Algorithm 1 is correct and explores exactly one node for eachapproximate solution in the instance.

Proof. The correctness of the algorithm is entailed by the properties of the permutation lattice described in Section 3. We merely show here that branches of the search tree pruned because of the set explored do not lead to unseen minimal solutions. A swap (a, b) is added to explored iff it has been previously tried in the same or in an earlier recursive call. In other words, let the current sequence be σ with characteristic set Γ (σ). If (a, b) ∈ explored, then a node π with characteristic set Γ (π) ⊇ Γ (σ) \ {(a, b)} and all its successors have been explored. Now suppose that there exists a successor ω or σ that is not a successor of π and such that (a, b) / ∈ Γ (ω). This implies that Γ (ω) ⊂ Γ (σ) and Γ (ω) ⊂ Γ (π). However, since Γ (π) ⊇ Γ (σ) \ {(a, b)}, we necessarily have (a, b) ∈ Γ (ω), contradicting our hypothesis.

Algorithm 1 never explores a sub-optimal node, that is, a node with objective value strictly greater than (1 + ) × f (σ * ). Indeed all predecessors of an -approximate node are -approximate, and a sub-optimal successor is not visited. Moreover, no node is explored twice. Indeed, suppose that a node ω has already been explored, and that the current node is a direct predecessor σ of ω. Let consider the recursive call at which ω was first explored. This call was at the same level in the recursive tree, and when exiting the branch starting with ω a swap (a, b) ∈ Γ (ω) is added to explored. Therefore, the branch leading to ω from σ will not be explored, since (a, b) ∈ Γ (ω) and (a, b) ∈ explored. Therefore every -approximate node is explored at most once (in fact exactly once since there is no other pruning condition).

Computational Experiments

We generated 30 random instances of two-machine flowshop problems (F 2||C max ) for 3 sizes: 8, 10 and 12 jobs. Then for each of them we computed an optimal sequence by sorting the jobs according to Johnson's rule. We give results for each size class, and for 5 values of , average across the 30 instances. The experiments were run on a MacBook pro dual core 2.3 GHz (no parallelization was used) using Mistral solver, except for the ILP model that was run on a PC clocked at 2.26 GHz with 3.45 GB RAM and solved using IBM ILOG CPLEX.

PB1: Finding a sequence of minimum level

In Table 1, we compare the results of the ILP and CP approaches to solve PB1, i.e., find a solution of minimum level. The first column indicates the number of jobs n, and the value of . A time limit of 10 minutes was imposed on both approaches. For each method, we report the average minimum level of the sequence, the ratio of proven results within the time cutoff, the average CPU time in seconds, and the number of nodes explored. We observe that despite the modest size of the instances, the problem is very hard to solve for the model we used. We believe that the objective function makes the problem extremely difficult to solve to optimality. These models can be efficient on much larger job-shop problems, however with more standard objective functions (such as C max or L max ).

We were unfortunately unable to run the two approaches on identical hardware, however, the two machines were of the same generation and the data we report here is sufficient to conclude that:

the CP approach is slightly faster than the ILP model on smaller instances.

however, when the size grows, and not all instances are solved to optimality within the time cutoff, CPLEX is able to prove optimality more often than Mistral on the CP model.

Another interesting observation is that despite the fact that the proposed approaches cannot prove it for all instances, they both find a sequence of minimum level. The fact that they both find the exact same objective value (level) is already a strong argument, and in fact we were able to verify it in the following set of experiments.

PB2: Finding all minimal sequences

In Table 2, we compare our lattice exploration algorithm with a simple enumeration of all the solutions of the problem using the CP model above. The first column still refers to the number of jobs n, and the value of . The number

Fig. 1 :

 1 Fig. 1: Lattices of permutations for n = 3 (left) and n = 4 (right)

2

  if a < b and (a, b) ∈ Γ (σ) and opt(swap(σ, (a, b))) then minimal ← False; 3 if (a, b) ∈ explored then Explore-Lattice(swap(σ, (a, b)), Γ (σ) \ {(a, b)}, explored); explored ← explored ∪ {(a, b)}; Algorithm 1 is initially called with σ = (123..n), Γ (σ) = {(a, b)/1 ≤ a < b ≤ n} and explored = ∅. It explores the optimal part of the lattice depth first.

Fig. 2 :

 2 Fig. 2: Elements of Γ (σ) in the lattice of permutations (n = 4)

Table 1 :

 1 Minimum level: ILP vs. CP

	Instance	CP Level Opt. Time (s) Nodes	ILP Level Opt. Time (s) Nodes
		= 0	14.40 1.00	0.05	5085	14.40 1.00	0.38	2048
		= 0.05	9.80 1.00	0.04	3945	9.80 1.00	0.32	1340
	n = 8	= 0.1	7.06 1.00	0.02	2597	7.06 1.00	0.24	823
		= 0.15	4.96 1.00	0.01	1510	4.96 1.00	0.30	392
		= 0.2	3.36 1.00	0.01	896	3.36 1.00	0.11	228
		= 0	23.36 1.00	3.27 249195	23.36 1.00	7.57	43679
		= 0.05	16.00 1.00	2.05 170201	16.00 1.00	4.65	28641
	n = 10	= 0.1	12.16 1.00	1.15	91903	12.16 1.00	3.28	19410
		= 0.15	9.06 1.00	0.59	46386	9.06 1.00	2.34	9300
		= 0.2	6.46 1.00	0.30	23350	6.46 1.00	0.92	4086
		= 0	30.26 0.66	72.64	3.83M	30.26 0.70 329.02 942822
		= 0.05	19.63 0.90	57.50	3.42M	19.63 0.93 141.71 506002
	n = 12	= 0.1	14.36 0.96	35.34	2.07M	14.36 1.00	56.85 281196
		= 0.15	10.46 1.00	18.10	1.06M	10.46 1.00	20.78 100744
		= 0.2	7.23 1.00	6.50	0.37M	7.23 1.00	6.53	27201

{(1, 2)(1, 3)} {(3, 4)} {(2, 3)} {(1, 3)} {}

of characteristic solutions (total number of minimal sequences) is reported in column 3 while the CPU time in seconds to compute them all using the lattice exploration algorithm is in column 2. Then, in columns 5 and 4 we give the total number of solutions (total number of -approximate solutions) and the time it takes to list them all using the CP model. Last, we give the ratio of the total number of solutions over the number of characteristic solutions in column 6.

We first observe that the ratio between the size of the whole set of -approximate solutions and the number of minimal solutions required to characterize them seems to grow geometrically both with n and with . It suggests that the benefit of this approach grows with the size of the instances. For n = 12 and > 0.05, the enumeration took longer than the one hour and a half time cutoff that we imposed in every case.

The CPU time needed to compute characteristic solutions with the lattice exploration algorithm is of course much smaller than that of enumerating all solutions. However the factor is relatively constant when the size augments. This is not surprising since the complexity of the lattice exploration algorithm is very much linked to the total number of optimal solutions. Moreover, we observe that finding all minimal solutions using Algorithm 1 is often faster than finding the solution of minimum level with either the CP or ILP model. Clearly, since the minimum solution is minimal, it suggests that the CP and ILP models are sub-optimal. However, as grows, the runtime required for both the CP and the ILP models for solving PB1 becomes lower than that of Algorithm 1 for finding all minimal solutions. Indeed, the complexity of solving these models does not depend directly on . In fact, we observe empirically that both CPLEX and Mistral are faster when grows. Therefore, the two approaches are complementary.

Conclusions and Further Research Directions

In this paper, we propose to characterize the whole set of -approximate schedules for the two-machine flow shop scheduling problem. The main concept to set up our reasoning is the lattice of permutations (total orders of the jobs). To treat the issue, two optimization problems are addressed: (1) finding an optimal sequence with minimum level (maximum depth) in the lattice; (2) finding all minimal sequences. The complexity of these problems remains open.

We propose three different approaches operating the concept of lattice of permutations, namely an integer linear programming formulation, a constraint programming approach, and a direct exploration of the lattice based on depthfirst search.

Computational experiments show that the CP approach with Mistral is faster than the ILP with CPLEX for small instances, but CPLEX can prove optimality more often than Mistral.

This approach can be used for characterizing the set of -approximate schedules for other scheduling problems, if an ordering relation exists between two consecutive jobs. It is the case for some single-machine problems (1||L max , 1|| C j , ...) and for other two-machine flow shop scheduling problems [START_REF] Bouquard | Application of an optimization problem in max-plus algebra to scheduling problems[END_REF].

A future research direction is to go further in the complexity study of the problems and to succeed to prove that they are actually NP-hard for the F 2||C max .

It is unlikely, however, that this method will be helpful for counting the number of optimal solutions. Indeed, counting the number of predecessors of a sequence in the lattice is equivalent to counting the linear extensions of a partial order, which is #P-complete [START_REF] Brightwell | Counting linear extensions[END_REF]. Furthermore, one needs to take into account the intersection of sets of predecessors, which makes the problem even more difficult.