Bernard Blavoux
email: bernard.blavoux@gmail.com

Patrick Lachassagne
email: patrick.lachassagne@danone.com

Abel Henriot
email: abel.henriot@danone.com

Bernard Ladouche
email: b.ladouche@brgm.fr

Vincent Marc
email: vincent.marc@univ-avignon.fr

Jean-Jacques Beley
email: jean-jacques.beley@danone.com

Gérard Nicoud
email: gnicoud@live.fr

Philippe Olive
email: olive87@wanadoo.fr

A Fifty-year chronicle of tritium data for characterising the functioning of the Evian and Thonon (France) glacial aquifers

Keywords: Evian, lumped model, natural mineral water, inverse modelling, tritium, glacial aquifers

to determine if important hydrologic information is available. The interest and limits of such a modelling, also for other constituents than tritium, along with the future for tritium as a tracer are discussed and it also provides new insights on the structure and functioning of alpine paleo glacial hydrosystems.

Highlights

-interpretation of an exceptionally long tritium chronicle (rainfall, natural mineral water and superficial springs) -accurate characterization of a hydrological transfer function -discussion of the interest and limits of lumped models -new insights on the structure and functioning of glacial aquifers

Abstract 1

Using lumped models and a transfer function model, this paper deals with the interpretation of exceptionally long (up to 50 years) and precise tritium chronicles characterizing the rainfall, recharge (efficient rainfall) and outflow from various types of glacial aquifers from the French Alps (Evian-Thonon area). The efficient rainfall tritium chronicle was computed from tritium measurements performed for 11 years (1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)) in a lysimeter. The evapotranspiration induces a mean 15% drop of the annual tritium signal. The three superficial glacial aquifers (two fluvio-glacial kame terraces and a lateral till) provide similar results: a best fit with an exponential flow model (EM) (playing the major role) combined in parallel with a piston flow model (PFM), and a rather short mean transit time (T: 5 to 7 years). The deepest mineral aquifer (Evian) can only be fitted with the in a series combination of a highly dispersive model (DM; T 68 years, DP=0.5) and a piston flow model (T 2.5 years) or, better, by the in a series combination of an EM (T 8 years) modelling the subsurface aquifer and a DM (T 60 y, DP=0.75) and the same piston flow model (T 2.5 years) modelling the deep mineral aquifer, this latest combination of models providing the following parameters: T 70 years and median transit time 45.5 years. It is also to be noted that a very small part of the efficient rainfall, about 1.3%, avoids both the EM and the DM, and directly enters the PFM (at the Northern limit of the Gavot Plateau). These models are very sensitive regarding the T (±1 y, 0.25 y for the PFM), less so with DP. These results will prompt hydrologists to (re)work historical data

Introduction -History of isotope measurements in the Evian-Thonon area

To evaluate the potential hydrogeological vulnerability of hydrosystems to anthropic pollution and predict any future changes, the pattern of the transit time of water in the aquifers should be defined. Isotopic measurements are presently used in field investigations and research projects to improve qualitative understanding of hydrosystems and identify the processes involved in water transfer [START_REF] Turner | Modeling of isotopes and hydrogeochemical responses in catchment hydrology[END_REF]. Environmental tracers ( 18 O, 2 H and tritium) and specific flow models have been used by numerous authors to estimate the transit time of water in hydrological systems [START_REF] Herrmann | The environmental tracer approach as a tool for hydrological evaluation and regionalization of catchment systems[END_REF][START_REF] Stewart | Modeling baseflow soil water residence times from deuterium concentration[END_REF][START_REF] Maloszewski | Application of flows models in an alpine catchment area using tritium and deuterium data[END_REF][START_REF] Maloszewski | Isotope hydrogeological study of mean transit time in alpine basin (Wimbacht, Germany)[END_REF][START_REF] Maloszewski | Isotope hydrology investigations in large refuse lysimeters[END_REF][START_REF] Maloszewski | Identifying the flow systems in a karstic-fissured-porous aquifer, the Schneealpe, Austria, by modelling of environmental 18O and 3H isotopes[END_REF][START_REF] Maloszewski | Determining the turnover time of groundwater systems with the aid of environmental tracers, I-models and their applicability[END_REF][START_REF] Holko | Stable environmental isotopes of 18 O and 3 H in hydrological research of mountainous catchment[END_REF][START_REF] Amin | A general lumped parameter model for the interpretation of tracer data and transit time calculation in hydrologic systems[END_REF][START_REF] Rodhe | Transit times for water in a small till catchment from a step shift in the oxygen 18 content of the water input[END_REF][START_REF] Vitvar | Estimation of mean water residence times and runoff generation by 18O measurements in a pre-alpine catchment (Rietholzbach, eastern Switzerland)[END_REF][START_REF] Etcheverry | Direct simulation of groundwater transit-time distributions using the reservoir theory[END_REF][START_REF] Uhlenbrook | Hydrograph separations in a mesoscale mountainous basin at event and seasonal timescales[END_REF][START_REF] Viville | Isotope hydrological study of mean transit time in the granitic Strengbach catchment (Vosges massif, France): application of the FlowPC model with modified input function[END_REF][START_REF] Hagedorn | Estimating recharge in fractured aquifers of a temperate humid to semiarid volcanic island (Jeju, Korea) from fluctuations, and Cl, CFC-12 and 3H chemistry[END_REF][START_REF] Molicova | Groundwater age for identification of baseline groundwater quality and impacts of land-use intensification -The National Groundwater Monitoring Programme of New Zealand[END_REF]. One difficulty in estimating the transit time of water is the need of the input function, i.e. the concentration of tritium in the infiltrating water while only tritium in precipitation and precipitation amount are available. [START_REF] Mcguire | A review and evaluation of catchment transit time modeling[END_REF] tried to construct better model to estimate the input function. Moreover, due to the insufficient length of input and output chronicles, hydrosystems with transit times of several decades sometimes cannot be studied.

Tritium ( 3 H) is a natural isotope of hydrogen (half-life 12.32 years; Lucas and Unterwerger, 2000), naturally present in the atmosphere at low concentrations (5 to 10 T.U.), and at higher concentrations (up to about 6000 T.U.) due to above ground nuclear weapons testing beginning in the mid 1940s, then dramatically increasing with the development of the hydrogen bomb (first test was in November 1952) until 1963, the year when the atmospheric atomic test was banned (the few small nuclear tests above ground since 1980 have just been too small to make any impact on the tritium burden). However the beginning of tritium recording started in August 1953 [START_REF] Maloszewski | Lumped parameter models for the interpretation of environmental tracer data, User's guide of Flowpc program[END_REF]. It has been widely used in hydrology to estimate aquifers' mean transit time (or residence time) since it is a constituent of the water molecule and its concentration in rainfall has varied during the past sixty years, giving a well-known input function for hydrological models. Tritium better captures the longer term components of flow and recharge than other typical approaches (chloride or stable isotopes for instance) [START_REF] Stewart | Truncation of stream residence time: How the use of stable isotopes has skewed our concept of streamwater age and origin[END_REF].

The first tritium measurements of groundwater from the Quaternary deposits of the Geneva Lake basin were taken in the fifties, and have been densified since 1964, shortly after monthly analyses of tritium from rainfall began in June 1963 at the Thonon (France, Fig. 1) Centre de Recherches Géodynamiques (CRG, a branch of the Paris VI University), which was a member of the International Atomic Energy Agency (IAEA) isotopic monitoring network. As early as 1963 this new isotopic tool proved useful in characterizing the hydrogeological functioning of the complex glacial aquifers of the area [START_REF] Blavoux | Hydrodynamique et teneur en tritium des eaux du bassin d'Evian[END_REF][START_REF] Fontes | Oxygène-18 et tritium dans le bassin d'Evian[END_REF]. A tritium content of several hundred TU (Tritium Units) was found in 1964 in the Thonon unconfined aquifer while in contrast the water from the Evian confined aquifer located about ten kilometres eastwards was still only a few TU. Sampling was subsequently performed quarterly in this area to monitor the evolution of the tritium signal until recently, providing probably the longest available detailed temporal record of tritium in groundwater.

In that framework, this paper deals with the interpretation of such an exceptionally long tritium chronicle characterising both the recharge and the outflow from various types of glacial aquifers, using lumped models and a transfer function model. It discusses the interest and limits of such modelling and also provides new insights on the structure and functioning of ancient alpine glacial hydrosystems. Such glacio-fluvial hydrosystems, either quite old, or "fossils" (see for instance [START_REF] Blavoux | Les sondages dans le complexe quaternaire du Bas-Chablais et leurs enseignements stratigraphiques. Leur intérêt pour l'hydrogéologie et l'hydrochimie régionales[END_REF][START_REF] Dray | Les terrasses de Thonon : aspects géologiques de la déglaciation würmienne et intérêt hydrogéologique[END_REF][START_REF] Caballero | Hydrological characteristics of slope deposits in high tropical mountains (Cordillera Real, Bolivia)[END_REF][START_REF] Nicoud | Les complexes détritiques de marge glaciaire active dans le Bas Chablais (Bassin lémanique, France). Implications hydrogéologiques[END_REF][START_REF] Henriot | Geological structure and hydrogeological functioning of a deeply incised glaciary margin complex (Gavot Plateau, Evian, France) from a pluridisciplinary approach[END_REF][START_REF] Bayer | Three-dimensional high resolution fluvio-glacial aquifer analog : Part 1 : Field Study[END_REF] or subactual (see for instance McClymont et al. 2011, Rodhe and[START_REF] Rodhe | Groundwater dynamics in a till hillslope: flow directions, gradients and delay[END_REF][START_REF] Clow | Groundwater occurrence and contributions to streamflow in an alpine catchment, Colorado Front Range[END_REF] are of high interest for their water resources in the many regions of the world that were glaciated, mostly during the Quaternary era, or are still at the margins of glaciated areas: the northern part of Eurasia, North and South America, alpine regions, etc. (see for example [START_REF] Roy | Multiple, distinct groundwater flow system of a single moraine-talus feature in an alpine watershed[END_REF]Beckers andFrind, 2000, 2001).

Geological and hydrogeological setting of the Evian-Thonon area

The studied area is about 20 km long by 10 km wide along the southern shore of Lake Geneva (Fig. 1 & 2). As in most of the Lemanic basin [START_REF] Fiore | Würmian deglaciation of western Lake Geneva (Switzerland) based on seismic stratigraphy[END_REF], the land surface is almost exclusively covered by thick (up to 400 m or more in Evian; [START_REF] Triganon | Contrôle de la construction de l'ensemble détritique de la région d'Evian par trois phases glaciaires durant le Würm[END_REF] Quaternary glacial sediments deposited over a major glacial erosion surface truncating the cemented, folded and thrusted Alpine rocks. The lake was dug during Quaternary glaciations by ice originating from the upstream Rhône valley catchment basin. The deep valley of the Dranse River (Fig. 1), the bottom of which exposes much older pre-Quaternary geological formations, divides the study area into two independent zones both from a geological and a hydrological viewpoint.

East of the Dranse River, the Evian area Quaternary geological formations are composed of glacial sediments deposited along the Rhône glacier's southern margin [START_REF] Nicoud | Les complexes détritiques de marge glaciaire active dans le Bas Chablais (Bassin lémanique, France). Implications hydrogéologiques[END_REF], well known from several kilometres of cored boreholes [START_REF] Blavoux | Les sondages dans le complexe quaternaire du Bas-Chablais et leurs enseignements stratigraphiques. Leur intérêt pour l'hydrogéologie et l'hydrochimie régionales[END_REF][START_REF] Blavoux | L'occupation de la cuvette lémanique par le glacier du Rhône au cours du Würm[END_REF][START_REF] Triganon | Contrôle de la construction de l'ensemble détritique de la région d'Evian par trois phases glaciaires durant le Würm[END_REF]

(Figs 1 & 2a):
-the shape of the pre-Quaternary substratum is the result of erosion from the pre-Wurmian pulsations of the Rhône glacier. This paleotopography is irregular, reaches an elevation less than 100 m amsl at Evian, and progressively rises towards the South, up to a few pre-Quaternary outcrops at the highest points of the Gavot Plateau; -this depression is filled by about 100 m of thick sands and gravels (older than the last glacial cycle), corresponding to the northern tip of a paleodelta of the Dranse River, covered by a subglacial till up to 200 m thick from the Würmian maximum (65'000 to 35'000 y BP, [START_REF] Triganon | Contrôle de la construction de l'ensemble détritique de la région d'Evian par trois phases glaciaires durant le Würm[END_REF], during which the glacier reached an elevation of about 1300 m amsl in the Evian area. This till surely prevents any groundwater flow between the pre-Quaternary substratum and the sediments described here below; -glacio-lacustrine sediments (small local streams paleodeltas within lacustrine silts) from the "Inferior Complex", older than 30'000 y BP, sedimented during a major retreat phase of the Rhône glacier, overlie this till up to about 390 m amsl; -the Gavot Plateau glacial margin formation, over 200 m thick, results from the progressive rise of the Rhône glacier during the Würmian "Lemanic stage" (30'000-27'000 y BP, [START_REF] Triganon | Contrôle de la construction de l'ensemble détritique de la région d'Evian par trois phases glaciaires durant le Würm[END_REF][START_REF] Fiore | Würmian deglaciation of western Lake Geneva (Switzerland) based on seismic stratigraphy[END_REF] when the glacier reached the 7 Geneva city area. It is composed of a subglacial till in the Evian area and up to 13 superimposed subglacial and lateral tills with peatbog sediments below the plateau; -the partial glacier retreat down to an elevation of about 480 m amsl, and then to about 430 m amsl at Evian, was responsible for the deposition of the glacio-lacustrine sediments of the "Terminal Complex" (25'000-21'000 y BP, [START_REF] Triganon | Contrôle de la construction de l'ensemble détritique de la région d'Evian par trois phases glaciaires durant le Würm[END_REF] with highly permeable (10 -3 to 10 -4 m/s) sands and gravels in the fans of the local paleostreams and low hydraulic conductivity silts laterally; -finally, the last flood of the Rhône glacier, with the western tip of the glacier again reaching Western Lake Geneva -"le Petit-Lac" -(Copet and Nyon readvances, [START_REF] Fiore | Würmian deglaciation of western Lake Geneva (Switzerland) based on seismic stratigraphy[END_REF], covered all these sediments with a final subglacial till (the "terminal till") of 30 to more than 50 m thick up to an elevation of 650 to 750 m amsl in Evian. This till locally comprises kame terraces, mostly deposited at the mouth of local streams during the final deglaciation, and lateral tills.

The Natural Mineral Water (NMW) emerging from the confined aquifer of the "Terminal Complex", particularly at the Evian-Cachat spring, the highest-yield spring in Evian at 400 m amsl, has been famous since 1789. Since that time, it has been widely used for its medicinal properties (hydrotherapy). The Société Anonyme des Eaux Minérales d'Evian now bottles "évian ® " NMW from these mineral springs. The terminal till prevents any infiltration of surface water into the aquifer and confines it in the downstream part of the aquifer. Recharge occurs on the Gavot Plateau, above the 650-750 m amsl maximum elevation of the last glacier advance [START_REF] Blavoux | Etude du cycle de l'eau au moyen de l'oxygène 18 et du tritium[END_REF][START_REF] Henriot | Geological structure and hydrogeological functioning of a deeply incised glaciary margin complex (Gavot Plateau, Evian, France) from a pluridisciplinary approach[END_REF]. Then the water slowly flows down and northward within the Gavot Plateau glacial margin formation that comprises several superimposed small confined and unconfined aquifers flowing mostly in the lateral tills, separated one from the other by the basal till layers that act as semi-confining units [START_REF] Henriot | Geological structure and hydrogeological functioning of a deeply incised glaciary margin complex (Gavot Plateau, Evian, France) from a pluridisciplinary approach[END_REF]. The mean elevation of the catchment area of the aquifer was estimated at 950-1000 m amsl by comparing the 18 O content of the Evian NMW aquifer with the 18 O content of precipitation at several altitudes and with the 18 O content of low discharge subsurface springs with a well-known catchment area altitude [START_REF] Blavoux | Etude du cycle de l'eau au moyen de l'oxygène 18 et du tritium[END_REF]. Mean annual rainfall on the catchment area rises to 1'200 mm, and effective rainfall to about 750 mm. The specific hydrochemical characteristics of the Evian NMW, calcium (80 mg/l) and magnesium (26 mg/l) (that dominate the cations) on one hand, and silica (15 mg/l) on the other hand, are respectively the result of water-rock interaction with rocks of the aquifer (limestone and dolomitic sands) that were eroded locally, and with rocks (granites, schists, green rocks…) from buried lateral tills that were imported from the Rhône valley catchment by the glacier.

The pH is neutral (7.2) and bicarbonates (360 mg/l) and sulphate (13 mg/l) dominate the anions. Due to the high inertia of the hydrogeological system, the mineral content and temperature (11.6°C) at the spring are perfectly stable (below the analytical thresholds) and its discharge shows less than 5% annual variations.

Small superficial aquifers are locally imbedded within the sands and gravels from the alluvial kame terraces interstratified in or deposited above the terminal till. These aquifers are recharged by the local precipitation. These waters can be easily distinguished from Evian Natural Mineral water as the mineral content is different from that of Evian NMW, for example in the Ca/Mg ratio and silica content. West of the Dranse River, the Thonon slope is mainly covered by the fifteen [START_REF] Gagnebin | Les terrains quaternaires des environs de Thonon (Haute-Savoie)[END_REF]) Dranse River glacio-fluvial kame terraces (sands, gravels, and -limestone and less numerous crystalline -large pebbles) progressively deposited westwards against lateral tills (Figs. 1 & 2b), parallel to the Rhône glacier, during glacier retreat phases. The highest terraces deposited during the Gavot Plateau retreat (Geneva stage: 25'000 to 23'000 y BP) and the lowest during the final deglaciation of the "Petit Lac" stage (21'000 to 18'000 y BP).

They all lie on a subglacial till contemporaneous with the formation of the Gavot Plateau ("Lemanic" stage), and the lowest lie on the till associated with the last flood of the Rhône glacier ("Petit Lac" stage), which there confines the oldest terraces. The terrace aquifer is generally unconfined and flowing northward. It exhibits a rather high hydraulic conductivity (about 5.10 -4 m.s -1 ; [START_REF] Blavoux | Etude du cycle de l'eau au moyen de l'oxygène 18 et du tritium[END_REF]. Numerous springs discharge from the terraces and some of them, "Les Blaves" for instance, are impounded for drinking water by Thonon town waterworks (Figs. 1 & 2b). The "Versoie spring" discharges from one of the lowest terraces at an elevation of 465 m amsl. Known since Roman times, it has been exploited as natural mineral water, known as "Thonon Water", for bottling and hydrotherapy since 1864. The "spring" is tapped with a 15 m deep radial well equipped with a submerged pump. The mean elevation of the Versoie spring catchment area has been estimated at about 600 m amsl from natural 18 O tracing [START_REF] Blavoux | Etude du cycle de l'eau au moyen de l'oxygène 18 et du tritium[END_REF][START_REF] Dray | Les terrasses de Thonon : aspects géologiques de la déglaciation würmienne et intérêt hydrogéologique[END_REF]. It corresponds to permeable kame terraces located between 650 and 465 m amsl, where mean annual precipitation reaches 1000 mm and effective rainfall 550 mm.

Material and methods -Sampling and analyses

Monthly samples of rain water have been collected since July 1963 at the CRG Thonon rain gauge station, at 385 m amsl on the shore of Lake Geneva and 10 km away from Evian (Fig. 1 and Fig. 3). Sampling consists of collecting and protecting from evaporation in a cumulative rain gauge all the rain falling on a 2000 cm 2 plastic rain gauge during each month. The water catchwork of the Evian-Cachat spring, built in 1911, reaches the sandy confined aquifer at the end of a 66 m long horizontal gallery crosscutting the till cover. As the aquifer is confined, the gallery is completely sealed and neither the aquifer nor the water is visible there. The artesian flow is limited by a valve to maintain the pressure. The samples were collected at a sampling tap located in the gallery; Natural Mineral Water is not treated, and of course not sterilised, consequently the samples were not subject to any treatment. Some sporadic measurements were made for comparison on the other boreholes of the same mineral hydrogeological system. The Versoie samples were collected on the pumping well from a tap located at the top of the pumping tube. The Maxima spring (Figs. 1 & 2) emerges from a superficial aquifer. It is located near the Lake's shore at 380 m amsl within a sandy kame terrace of the terminal till on the lowest part of the Evian area; its more than 40 years of tritium data are worth studying. The Larringes-Pugny well is located on the Gavot Plateau at 790 m amsl. This well exploited for drinking water taps a superficial aquifer of the Gavot Plateau glacial margin formations (4 m thick sand and gravel locally covered by 2 m of clay).

All the water samples were analysed for tritium at the Centre de Recherches Géodynamiques (CRG) in Thonon until it was closed in April 2005, and since then have been analysed at the Hydrogeology Laboratory of the University of Avignon. Tritium measurements were obtained by liquid scintillation counting after electrolytic enrichment. The associated error on the data varies with the tritium concentration in the samples; it is 1 TU for recent samples and was about 1 or 2 tens of TU for the samples which contained several hundred TU (see Figs 4A, 4B and 5). The limit of detection of the CRG equipment was 5 to 8 TU during the earliest years and about 1 TU since 1972 (see Fig. 5).

Theory/calculation -Tritium restitution modelling

Lumped parameter models are used in this study to estimate the transit time distribution of the groundwaters discharging from the Evian -Thonon area. In the lumped parameter approach, the output concentration is related to the input concentration using the following convolution integral (Equ. 1) [START_REF] Maloszewski | Determining the turnover time of groundwater systems with the aid of environmental tracers, I-models and their applicability[END_REF]Zuber, 1986(a) and (b)):

ܥ ௨௧ ሺݐሻ =  ܥ ሺݐ -ݐ ᇱ ሻ݁ ିఒ௧ ᇲ ݃ሺ′ݐሻ݀ݐ ᇱ ஶ (1)
where ݐ is the sample date, ′ݐ is the transit time of single tracer particle through the system, ߣ is the decay rate of the tracer, per year, C in _ and C out _ are the input and output tracer contents, and g(t') is the weighting function defining the transit time distribution.

Lumped parameter models assume a steady-state flow and that the selected tracers move with the water molecules, which is why tritium is an ideal tracer. Such a steady-state flow assumption is suitable for modelling the Evian and Thonon area aquifers; in fact, the water volume in these aquifers does not change significantly on a seasonal and yearly time scale as the variability in interannual discharge measurement is less than 5 % of the total annual discharge for each spring. In addition, as two of these springs (Evian-Cachat and Versoie) are bottled NMW, the spring discharges are monitored. In this paper, it is also considered that the aquifer system is mono-porous, thus that the portion of immobile water is negligibly small, which permits to use the mean transit time found from tritium as that of water.

The following lumped parameter models were considered [START_REF] Yurtsever | An overview of conceptual model formulations for evaluation of isotope data in hydrological systems[END_REF]IAEA, 1996): the piston flow model (PFM), exponential flow model (EM), linear model (LM), and dispersive model (DM). According to the tritium restitution curves and the hydrogeological conceptual models (see chapter 2.), some of these models were also combined in parallel, thus constituting Binary Mixing Models (a terminology proposed for instance by [START_REF] Katz | Timescales for nitrate contamination of spring waters, northern Florida, USA[END_REF], which means a simple mixing of two -or more -endmembers occurring at the spring/outlet. Piston flow and exponential model (named after BMMEP) were combined.

Some of these models were also combined in series. All these models take into account the natural radioactive decay of tritium. The choice of model/combination of models was made as a function of the appearance of the observed signal at each spring and the available knowledge on their hydrogeological functioning. Among the different computer codes available in the hydrological literature (for example FLOWPC, Maloszewski andZuber 1994, 1996), we first used LUMPED [START_REF] Ozyurt | LUMPED: a Visual Basic code of lumped-parameter models for mean residence time analyses of groundwater systems[END_REF], available at http://www.sukimyasilab.hacettepe.edu.tr/english/software.shtml). In a second phase, the complex signal measured at the Evian-Cachat spring was interpreted using the TEMPO software developed at BRGM by [START_REF] Pinault | Manuel utilisateur de TEMPO : logiciel de traitement et de modélisation des séries temporelles en hydrogéologie et en hydrogéochimie[END_REF] (see for instance also [START_REF] Dörfliger | Inverse modeling approach to allogenic karst system characterization[END_REF], in which the equations of the models described by Maloszewski andZuber (1994, 1996) have been implemented (the models developped by Maloszewski and Zuber have also been used to validate TEMPO). FLOWPC and TEMPO have the same input requirements and provide similar outputs. The subsurface aquifers of the Versoie, Larringes-Pugny, and Maxima springs were characterised using the exponential model (EM), and finally with a Binary

Mixing Model combining exponential and piston flow models where the exponential part plays the major role. The Evian-Cachat aquifer was characterised with an in a series model combining dispersive and piston flow models (Fig. 10), and secondly with an in a series model combining exponential, dispersive and piston flow models (Fig. 13).

The piston flow model (PFM) assumes that all the flow lines have the same transit time and thus, the response function is the Dirac delta function g(t)=δ(t-Τ), with Τ: mean transit time.

Hydrodynamic dispersion and diffusion are negligible (Maloszewski andZuber, 1982, 1994).

13

Tracer concentration is only affected by radioactive decay. It ideally corresponds to a confined aquifer with a narrow recharge area far from the sampling point.

In the EM, the distribution function of the transit time of water in the aquifer is exponential (g(t)= T -1 exp(-t/T), with T = mean transit time, thus the average age of all parcels of water that entered the system and reached the discharge point). Water travelling along the shortest flow line has a transit time equal to zero while water travelling along the longest flow line has an infinite transit time. An exponential distribution of transit times in an aquifer can corresponds to multiple hydrogeological situations, including those where hydraulic conductivity decreases with depth or the aquifer has uniform permeability and porosity [START_REF] Eriksson | The possible use of tritium for estimating groundwater storage[END_REF][START_REF] Maloszewski | Determining the turnover time of groundwater systems with the aid of environmental tracers, I-models and their applicability[END_REF]. Use of this model assumes that (1) there is no mixing of water among the various flow lines during transfer in the hydrosystem and (2) mixing of the water from the various flow lines occurs at the level of the sampling point.

The dispersive model (DM) allows the effects of transit time dispersion in the flow direction to be taken into account. The distribution function of transit times in this model uses the onedimensional solution of the dispersion equation for a semi-infinite medium (Equ. 2) [START_REF] Maloszewski | Determining the turnover time of groundwater systems with the aid of environmental tracers, I-models and their applicability[END_REF]:

݃ሺݐሻ = ቀ ସగ ௧ ் ቁ ି భ మ ݔ݁ - ఛ ସ ் ቀ1 - ௧ ఛ ቁ ଶ ൨ ݐ ିଵ (2)
with T = mean transit time and DP (without dimension) the dispersion parameter: the dimension of variance of the dispersive distribution of the transit time. This additional calibration parameter DP depends mainly on the distribution of transit times of the water (shape of the flow lines) in the aquifer. It thus includes apparent dispersion of the tracer in the system plus also apparent dispersion of the tracer in the recharge zone (recharge zone is in the lumped parameter approach reduced to one point, while space distribution of the input signal through the surface of the recharge zone, yields additional "dispersion"), and also the variable infiltration conditions over the recharge area. DP is generally greater (by one or more orders of magnitude) than the hydrodynamic dispersion parameter which usually represents only the dispersive properties of the material [START_REF] Maloszewski | Lumped parameter models for the interpretation of environmental tracer data, User's guide of Flowpc program[END_REF].

The distribution functions of transit times providing the best correspondence of the simulated concentrations to the experimental data were obtained by trial and error until a graphically good fit and a "good" Nash coefficient were obtained (at best, the Nash coefficient is 1).

Results

Observation results

The annual tritium levels of the precipitation at Thonon CRG (Table 1) have been computed

from the monthly values weighted by the precipitation rate. For values prior to 1963, a regression (Equ. 3) was established with the annual levels measured in Ottawa from 1964 to 1985, a station at the same latitude as Thonon that has operated since 1953 [START_REF] Brown | A fifty-year chronicle of tritium data for characterising the functioning of the Evian and Thonon (France) glacial aquifers[END_REF].

[ 3 H] Thonon = 0.892[ 3 H] Ottawa + 9, n=22 , R²=0.9 [START_REF] Molicova | Groundwater age for identification of baseline groundwater quality and impacts of land-use intensification -The National Groundwater Monitoring Programme of New Zealand[END_REF] ( between 30 and 20 TU until 1996; this is due to the feedback from the oceans and surface waters that contain higher levels of tritium [START_REF] Fine | Source fonction for tritium transport models in the Pacific[END_REF][START_REF] Fine | The penetration of tritium into the tropical Pacific[END_REF][START_REF] Michel | Residence times in river basins as determined by analysis of long-term tritium records[END_REF] and also the influence of the civilian nuclear industry. Since 2000, it seems that the rain signal at Thonon has entered a permanent regime around a value of 9 TU. The tritium content in rainfall shows a seasonal variation. Monthly levels in Thonon have a maximum in June-July and a minimum in November-December. For the temperate climate zone of the Northern hemisphere, this springtime increase results from direct discharge of tritium in the stratosphere toward the troposphere due to the discontinuities of the tropopause [START_REF] Olive | Contribution à l'étude géodynamique du cycle de l'eau dans l'hémisphère nord par la méthode du tritium[END_REF]. The tritium-rich rain of the summer months is much more affected by evapotranspiration than that of autumn-winter, poor in tritium, therefore the transformation of the rain signal upon infiltration must be investigated [START_REF] Foster | The interpretation of tritium in the chalk unsaturated zone[END_REF]. A 1 m 2 lysimetric box covered with grass and filled with fluvio-glacial earth 1 m deep was installed in Thonon at the site of the rain gauge. The drainage of this box, equivalent to the effective rainfall, was monitored quantitatively and qualitatively on a monthly basis from 1969 to 1979. During these 11 years, the average annual rainfall was 944 mm and its weighted tritium content was 132 TU, while the drainage at the base of the box was 537 mm and showed a weighted annual tritium content of 114 TU. The average inter-annual evapotranspiration of 407 mm thus affects the tritium rain signal by a drop of 15%. This mechanism was repeated each year in similar proportions during the 11 years of monitoring.

Consequently, the annual tritium levels in precipitation used as input data for modelling were systematically adjusted by a drop of 15% to estimate the infiltration signal toward the unconfined aquifers (Table 1). This correction is only really significant for the high tritium values close to the 1963 peak (Fig. 6). It seems illusory to specify the infiltration signal more precisely, in that the lumped models involve a steady-state regime that does not take into account the inter-annual variability of input volumes.

Results of modeling

Given the relatively large variation of the tritium signal at the Versoie spring, at least during the initial years of monitoring [START_REF] Blavoux | Hydrodynamique et teneur en tritium des eaux du bassin d'Evian[END_REF]( -1974) (Figs. 4) (Figs. 4), the calculated annual tritium levels at the spring (Table 1) were averaged from the numerous individual levels measured on samples taken during the initial years of monitoring, every month until 1971, then every quarter until 1979. Finally, far from the 1963 thermonuclear peak in precipitation, the values used correspond to a single annual check. The shape of the observed tritium signal (Fig. 7) can be described using an exponential model (EM), taking a mean transit time (T) of 5 years (Fig. 7).

This value is close to that of 4 years put forward by [START_REF] Molicova | Groundwater age for identification of baseline groundwater quality and impacts of land-use intensification -The National Groundwater Monitoring Programme of New Zealand[END_REF] using a Gamma law on a shorter series, which was later confirmed by [START_REF] Hubert | Modélisation par une loi Gamma de la distribution des temps de séjour de l'eau dans des systèmes hydrologiques en régime permanent[END_REF] and [START_REF] Olive | Estimation pratique de l'âge des eaux souterraines en Europe par le tritium[END_REF]. A PFM only is definitely not appropriate (Fig. 7). Nevertheless, the signal measured at the Versoie spring shows a pronounced peak as of the beginning of monitoring in 1965, which probably indicates the existence of a component of water transfer by piston effect. Use of a Binary Mixing Model combining in parallel an EM with a PFM (BMMEP) thus allows the measured concentrations to be well reproduced (Fig. 7) with a 5.1 years mean transit time, and with the volume and flow parameters for the two reservoirs shown in Table 2. This transfer model is perfectly compatible with the hydrogeological context of the Versoie spring ( §2. and Fig. 2) that involves two aquifers: (1) most of the flow (about 80%) comes from an unconfined aquifer that is attributed to the paleo-Dranse river lower fluvio-glacial kame terraces resting on the Petit Lac Stage impervious morainic substrate, and (2) a lower part of the flow (about 20%) is attributed to the earlier terraces aquifer, which is confined by the Petit Lac Stage impervious basal till and appears to be recharged upstream.

Figure 7 -Modelling the Versoie spring tritium time series (using LUMPED, [START_REF] Ozyurt | LUMPED: a Visual Basic code of lumped-parameter models for mean residence time analyses of groundwater systems[END_REF] Table 2. -Summary of the parameters used to fit the Evian-Thonon area springs. Note that for the exponential-piston flow model, the two models are placed in parallel

The tritium signal measured at the Maxima spring of Amphion (380 m amsl), although including less data, has an appearance very similar to that of the Versoie spring (Figs 4). The same Binary Mixing Model (BMMEP) provides a very good fit (Fig. 8). Neither the PFM nor the EM are able to simulate the observed data. The simulation provides a mean transit time of 6.7 years (Table 2). This result is in good accordance with the hydrogeological structure of the Maxima spring, which outflows from a shallow semi-captive aquifer (local kame terrace) interstratified within the terminal basal till that appears to comprise both a confined and an unconfined component. The water molecule stable isotope content shows that this aquifer is very local (mean elevation of the recharge area of about 450 m amsl).

Figure 8 -Modelling the Maxima spring tritium time series (using LUMPED, [START_REF] Ozyurt | LUMPED: a Visual Basic code of lumped-parameter models for mean residence time analyses of groundwater systems[END_REF] The tritium signal at the Larringes-Pugny well has been simulated with the same BMMEP with a good fit (Fig. 9) insofar as the scarce data available allow input and output data to be compared and a mean transit time of 6.3 years (Table 2). About 80% of the discharge of the spring is explained by the EM, while the last 20% is the result of PFM. The well also taps a shallow subhorizontal semi-confined aquifer flowing into a Gavot Plateau lateral till lying on a basal till which constitutes its impervious substratum. Hydrogeological data show that this aquifer is local, with a mean elevation of the recharge area of about 820 m amsl, similar to that of the well (790 m amsl). As the available data at this well are scarce, they can also be well fitted with an EM with T= 7 years (Fig. 9 and Table 2), which is not incompatible with the hydrogeological context (only locally confined aquifer). The signal is excessively attenuated; the maximum observed tritium levels are approximately 20 times lower. As it propagates over a greater time period, it appears to result from a transfer in a highly dispersive system. Moreover, the response of the hydrogeological system to infiltration of "tritiated" water is bimodal. Even if, due to their low activities, tritium measurements show a high uncertainty, an initial maximum, low in absolute value, is observed around 1966; the second, of significantly greater amplitude and very extended in time, appears much later, toward the 1980s. These two responses of the hydrogeological system to recharge by tritiated effective rainfall are the result of two different functional processes of the reservoir: (i) mainly a very highly dispersive aquifer system with very slow flow dynamics and (ii) a modest contribution of water with a less slow flow dynamics that doesn't flow through the highly dispersive aquifer (Fig. 10). After several trials, a conceptual model with two models in a series (a DM and then a PFM) (Fig. 10) was adopted. In addition, a modest (about 1.3%) contribution of efficient rainfall directly enters the PFM and consequently avoids the DM. The other models tested are not acceptable; in particular, a DM significantly overestimates tritium levels during the decline that follows the initial maximum of 1966. Although not perceptible with the classical hydrochemical approaches, the computed parameters of this first component are consistent with the hydrogeological knowledge:

-the PFM and the perfect mixing before the outflow is consistent with the fact that the aquifer of the "terminal complex" is confined by the terminal till at and upstream from the spring. An EM (unconfined aquifer) would have been totally implausible;

-with the known parameters of this aquifer in the vicinity of the spring (hydraulic conductivity about 2.10 -3 m.s -1 , hydraulic gradient about 1.10 -3 , effective porosity about 5%), a 2.5 years mean transit time for the first component yields a recharge area about 3 km from the spring. This order of magnitude of a few kilometres (northernmost part of the Gavot Plateau) is in agreement with the maximum southern extension of the terminal till that, for a maximum elevation of about 650 to 750 m amsl ( § 2.), respectively, lies at a horizontal distance of about 1.5 to 1.9 km from the spring (Figs. 1 & 2a). Consequently, this first component is very probably a very small part of the total recharge (about 1.3 %) that infiltrates on the northernmost part of the Gavot Plateau and (this second condition is very important) is able to rapidly (in about one year) reach the "terminal complex" high hydraulic conductivity aquifer; in fact, as the Gavot Plateau is composed of several superimposed basal and lateral tills, its mean hydraulic conductivity is quite low and this explains, with the capping terminal till that avoids any infiltration downstream, the long transit time of Evian water (see the following part of this paper). Consequently, the infiltration of this first component 21 occurs through few layers of superimposed basal and lateral tills, which is consistent with the geological structure of the northernmost tip of the Gavot plateau (Fig. 2a).

The effect of this less slow flow component is significant up to the mid-1980s; however, it is less than 5 TU starting from the mid-1970s.

Different types of isotopic transfer models (PFM, EM, BMM combining a PFM with an EM, DM) were tested to simulate the long term second part of the observed signal. Only the DM allows the restitution curve to be satisfactorily reproduced (Fig. 12) with T = 68 (±2; Fig. 12A) years, DP = 0.5 (±0.025; Fig. 12B), and a flow of 98.7% of the total discharge in the DM, the 1.3% remaining corresponding to the less slow flow component presented above.

The sensitivity analysis shows that the DM is sensitive to small variations of both T and DP (Fig. 12A andB). This behaviour is in good agreement with the hydrogeological knowledge.

The mineral aquifer recharge area, the 200 m-thick (Fig. 2a) Gavot Plateau, is composed of up to 13 superimposed Würmian basal and lateral till sequences. These formations respectively exhibit very low and quite low hydraulic conductivity. Hydrogeological and hydrochemical data from deep piezometers drilled on the plateau [START_REF] Henriot | Geological structure and hydrogeological functioning of a deeply incised glaciary margin complex (Gavot Plateau, Evian, France) from a pluridisciplinary approach[END_REF] show the presence of several superimposed quite low hydraulic conductivity aquifers (flowing in the lateral tills)

with hydraulic heads decreasing with depth, which is a strong clue for a vertical component of the groundwater flow. The groundwater is thus able to flow down through the basal till layers These subsurface aquifers constitute the first step in the functioning of the mineral hydrosystem. In this case, the tritium input signal of the deep mineral reservoir is transformed by its stay in the subsurface aquifer and so can in fact be different from the tritium chronicle resulting from effective rainfall. To simulate this new input signal, we thus assumed that the functioning of the various superficial aquifers of the Evian watershed area (Gavot Plateau) can be described using a mixing model (EM) with a transit time on the order of several years (Fig. 13). Various simulations were performed, taking among others mean transit times in the subsurface aquifer increasing from 6 to 10 years (Fig. 14). The monitoring of tritium levels in one of these subsurface aquifers at Larringes-Pugny from 1966 to 1969 is actually compatible with this choice (Fig. 9).The use of these data from the EM as new input data for the Evian-Cachat dispersive model (the two models are in a series) requires the parameters of this latest model to be revised. The tritium levels in the Evian-Cachat spring can only be well simulated if it is taken into account that the mean transit time in the EM is about 8 years (Fig. 14A), with a T of 60 (±1) years (Fig. 14B) in the deep aquifer system, DP = 0.75 (±0.025; Fig. 14C), and still a flow of 98.7% of the total discharge in the EM and DM, the 1.3% remaining corresponding to the less slow flow component presented above that only flows in the PFM (Figure 13). The first result (8 years in the EM) is in good accordance with the fit obtained at the Larringes-Pugny well (7 years) tapping one of the Gavot Plateau superficial aquifers (see section 5. and Fig. 9). Consequently, the mean transit time in the whole Evian hydrosystem computed from the 3 models in a series is about 69. -during the following 14 years the difference between the two models is rather low.

The effect of the EM is nevertheless visible, with, during the first 8 years (up to a total duration of about 12 years), a restitution higher for M2 (EM8 + DM 60) than for M1, mostly showing the influence of the EM and of the lower T of the DM;

-after that, the differences are governed by the DM with logically a higher peak for M2.

Due to the difference of the two curves shapes, even if the T of both M1 and M2 are rather similar (about 70 years), their median are significantly different with a higher one for M1 (48 y) than for M2 (45.5 y). Consequently, the two different models provide significantly different results particularly during some important phases of the output, and particularly during the second part of the restitution (beyond 18 years) with significant absolute differences reaching about -0.12% (M1-M2), and consequently relative differences of about 8% ((M1-M2)/(Average of M1,M2)). 6. Discussion

On the need of anthropogenic tracers long term database in hydrological studies

It appears that, even in quite heterogeneous and complex aquifers such as those of the Evian-Thonon area, the use of quite simple lumped models leads to very accurate results. Such tools allow transfer model types to be obtained, their parameters inferred, and aquifer transit time distributions computed.

Use of the exceptionally long and precise Evian-Thonon rainfall, effective rainfall and spring tritium content time series allows existing knowledge about the structure and functioning of such complex aquifers to be much better characterised and/or assessed. Among other results, the computed median transit time in the Evian-Cachat aquifer, about 45 years, appears to be much longer than previously expected from simpler approaches. Knowledge of this distribution function of transit times, allowing mass transfer of chemical elements within the hydrosystem to be imaged, has been especially useful in establishing a policy to protect the mineral resource (see for instance [START_REF] Buric | Payment for environmental services: first global inventory of schemes provisioning water for cities[END_REF], Devier et al., 2013;[START_REF] Lachassagne | The 20 years technical and socio-economic Evian experience conciliating groundwater quality preservation, collective responsibility for environment protection and local development, and its transposition to other Danone water resources in the world[END_REF].

The long tritium chronicle available at the Evian-Cachat spring even allows a 1.3% contribution to the spring discharge that, due to its low value, was not identifiable by other hydrodynamical or hydrogeochemical means. However, this identification results from a very opportune and exceptional combination of circumstances: the fact that tritium monitoring at the Evian-Cachat spring began only a few years (2 years) after the main period of tritium injection into the atmosphere, and was carried out continuously for several decades. The scientists at the Centre de Recherches Géodynamiques were able to convince stakeholders to carry out such monitoring.

These results should encourage the owners of the old tritium database, even if sparse or incomplete, to make use of them, at least with such simple computational models. It should also encourage young scientists to be on the lookout for emerging tracers, and consequently to acquire the corresponding data on such high-stakes aquifers. A few years ago, noble gases and chlorofluoromethanes were such emerging tracers. These methods are presently being tested on the Evian aquifer and the results will be compared to those of this study.

This research also shows that the standard simple models can easily be improved with some kind of plug-in, for instance, to compute the behaviour of superficial aquifers. The next steps at Evian will be to more precisely compute the role of the high hydraulic conductivity aquifer of the "Terminal Complex" in the framework of a complete process-based modelling of the aquifer.

What is the future for Tritium as a tracer?

The tritium tool was especially effective in the vicinity of the 1963 rainfall peak and up to approximately 1980. But new prospects are emerging. Since 2000, the mean annual tritium signal in rainfall at Thonon has oscillated around a level of 9± 1TU. Comparison with the data available at the Orléans-BRGM station (Parisian Basin) shows the same trend around a level of 6 TU, and good correlation. Thus there seems to be a return to the case of a constant influx of tracer. In this case, the transit time of the water would be accessible starting from one or more tritium measurements using simple empirical models that require only this parameter;

here, the exponential model (EM) and the Piston Flow model (PFM). For unconfined aquifers likely to function according to an exponential model (of course if the vadose zone is not too thick, in order to allow the arrival of water with a theoretically zero residence time), the tritium level could thus be translated into the apparent age of the water, on the express condition of having good knowledge of the local input rainfall signal. Given the low levels and the necessary precision, correlations with a too-distant rain gauge station are now inadequate, and a denser network would be necessary. This projection gives 6.7± 0.7 TU for Versoie in the immediate future. In 2009, its content was measured at 7.1 ± 0.6 TU.

Measurements on such systems already studied near the 1963 peak would be worth repeating today to establish these new prospects. It would be interesting to verify whether slight differences still exist between families previously defined by their tritium level during former regional prospecting, so as to offer new standards for interpretation.

Conclusion

This paper deals with the interpretation, with lumped models and a transfer function model, of With this study the interest and limits of such a modelling and the future for tritium as a tracer can be discussed and it also provides new insights on the structure and functioning of alpine paleo glacial hydrosystems.

Figure captions

1

  Abbreviations: amsl: above mean sea level (elevation) CRG: Centre de Recherches Géodynamiques DM: dispersion model DP : dispersion parameter EM: exponential flow model BMMEP: in parallel combined piston flow and exponential models IAEA: International Atomic Energy Agency LM: linear model LPM: combined linear flow and piston flow model NMW: Natural Mineral Water T: Mean Transit Time PFM: piston flow model TU: Tritium Unit y BP: year Before Present

Figure 1 -

 1 Figure 1 -Location and schematic geological map of the Evian-Thonon area. Location of

Figure 2 -

 2 Figure 2 -Schematic geological cross-sections of the Evian (2a -top) and Thonon (2b -

Figure 3 CRG

 3 Figure 3 CRG Thonon tritium monthly rainfall time series (since june 1963; previous annual

  This chronicle of annual levels constitutes the input signal for the hydrogeological systems of the Evian-Thonon area (Figs.4A and 4B). Since the start of thermonuclear weapons testing in the early 50's, continued more and more intensively until 1963, atmospheric tritium levels increased enormously. The levels in rainwater at Thonon reached several thousand TU in 1963, with a spectacular rise between 1961 and 1963 (Fig.4B). From 1963 to 1966, nearly all atmospheric nuclear tests were stopped and levels dropped sharply with a pseudo-period close to one year, thus defining a very characteristic peak on the multiyear scale. From 1967 to 1980, atmospheric nuclear tests were of very low capacity, with less than 3% in total of all the tritium produced artificially, nevertheless leading to fallout that decreased slowly from 200 to 100 TU. Since 1980, there have been no atmospheric nuclear tests. With a pseudo-period of one year and according to tritium natural radioactive decay (half-life: 12.32 years), the level should have reached the natural level (5 to 10 TU) around 1985. However, it remained

Figure 5 -

 5 Figure 5 -Evian-Cachat spring tritium chronicle (note the increase in measurement

Figure 9 -

 9 Figure 9 -Modelling the Larringes-Pugny well tritium time series (using LUMPED, Ozyurt

Figure 10 -

 10 Figure 10 -Conceptual model of the Evian-Cachat spring (in a series dispersive and piston

Figure 11

 11 Figure 11 Modeling (TEMPO) of the Evian-Cachat spring. Fig. 11A: Sensitivity analysis of

(

  either with leakage flow and/or in places where the basal till layer(s) is (are) missing). The hydrogeochemical data sampled in these deep piezometers ( 13 C and 14 C data, but also and mainly silica content, and Mg/Ca ratio, which are all residence time indicators) also show that most of the NMW mineralisation process occurs within this Gavot Plateau multilayer low hydraulic conductivity system, which is not capped by the thick low hydraulic conductivity superficial terminal till deposited below 650 to 750 m amsl. Consequently, such large dispersion behaviour is the result of the various pathways the future NMW can follow within the Gavot Plateau formation. The water then fast flows within the Terminal Complex aquifer. Consequently, with this set of parameters, the mean transit time in the Evian hydrosystem is about 70 years (T: 70.05 years): mostly (about 68 years) within the subsurface and deep aquifers of the Gavot plateau (DM) and about 2.5 years within the high hydraulic conductivity Terminal Complex (PFM). Statistical processing of the output data shows that the median transit time of the Evian NMW (a less biased estimator than the arithmetic mean for such non Gaussian distributions) is about 48 years (48.0 years).

Figure 12 -

 12 Figure 12 -Modelling (TEMPO) and sensitivity analysis of the slow flow component of the

  5 years: 8 years within the subsurface aquifers of the Gavot plateau (EM), 60 years within the deep aquifers of the Gavot plateau (DM) and 2.5 years within the Terminal Complex (PFM). The median transit time of the Evian NMW is then 45.5 years. More in detail, the comparison of the results obtained with both in a series modelling approaches ((Model 1 -M1) DM+PFM and (Model 2 -M2) EM+DM+PFM) for a perfect tracer (without the tritium radioactive decay) shows interesting results (Fig. 15): -the 4 first years are governed by the PFM and show no difference between M1 & M2;

Figure 13 -

 13 Figure 13 -Conceptual model of the Evian-Cachat spring (in a series exponential, dispersive

  exceptionally long (up to 50 years) and precise tritium chronicles characterizing the rainfall, the recharge (efficient rainfall) and the outflow from various types of glacial aquifers from the French Alps (Evian-Thonon area).The evapotranspiration induces a mean 15% drop of the tritium signal. The three superficial glacial aquifers (two fluvio-glacial kame terraces and a lateral till) provide similar results: a best fit with an exponential flow model (EM) (playing the major role) combined in parallel with a piston flow model (PFM), and rather short mean transit time (T: 5 to 7 years). The deepest mineral aquifer (Evian) can only be fitted with the in a series combination of a highly dispersive model (DM; T 68 years, DP=0.5) and a piston flow model (2.5 years) or, better, by the in a series combination of an EM (T 8 years) modelling the subsurface aquifer and a DM (T 60 y, DP=0.75) and the same piston flow model (2.5 years) modelling the deep mineral aquifer (note that a very small part of the efficient rainfall, about 1.3%, avoids the DM (and the EM), and directly enters the PFM). This latest in a series combination of models is rarely used in hydrology even if it is conceptually highly credible in such a hydrogeological context.It provides the following parameters: T 69.5 years and median transit time 45.5 years. These models are very sensitive regarding the T (±1 y, 0.25 y for the PFM), less so with DP.

Figure 1 -

 1 Figure 1 -Location and schematic geological map of the Evian-Thonon area. Location of

Figure 2 -

 2 Figure 2 -Schematic geological cross-sections of the Evian (2a -top) and Thonon (2b -

Figure 3 -

 3 Figure 3 -CRG Thonon tritium monthly rainfall time series (since june 1963; previous annual

Figure 4 -

 4 Figure 4 -Rainfall and springs tritium mean annual content at the Thonon CRG rain gauge

Figure 7 -Figure 8 -Figure 9 -

 789 Figure 7 -Modelling the Versoie spring tritium time series (using LUMPED, Ozyurt and

Figure 12 -

 12 Figure 12 -Modelling (TEMPO) and sensitivity analysis of the slow flow component of the

Table 1

 1

	-(i) Rainfall tritium mean annual content at the Thonon CRG rain gauge, (ii)
	computed tritium efficient rainfall signal, and (iii) mean annual tritium content at the Evian-
	Thonon area springs
	Figure 4 -Rainfall and springs tritium mean annual content at the Thonon CRG rain gauge
	and at the Evian-Thonon area springs (top: logarithmic scale; bottom: arithmetic scale)

Acknowledgments

The authors are grateful to BRGM and Evian Volvic World Sources, which co-funded this research project, but also to Paris VI University and colleagues from CRG Thonon who greatly contributed to the development of such long time series. The constructive and detailed comments of the reviewers greatly helped to improve the quality and the strength of this paper. We are grateful to E. Le Lann for revising the English text and translating some early parts written in French.


Ph Bolon

E Trouvé

I Pétillot

G Vasile

M Gay

L Bombrun

J.-M Nicolas
email: nicolas@tsi.enst.fr

F Tupin

T Landes
email: tania.landes@insa-strasbourg.fr

M Koehl

P Grussenmeyer

Monitoring Alpine Glaciers with ALOS SAR and Optical data

Keywords: SAR differential interferometry (D-InSAR), SAR polarimetry, temperate glacier

The monitoring of temperate glacier activity is one of the applications which require the combined use of optical and SAR data. It should become feasible thanks to the ALOS PRISM and PALSAR sensors. In this paper, we present the interest of using L-band polarimetric interferometric SAR data and high resolution panchromatic stereoscopic data to monitor moving temperate glaciers made of ice, snow and rocks. The scientific issues are described and a test-site located in the Mont-Blanc region in the Alps is proposed for specific experiment and validation. The Mont-Blanc test-site includes two well-known glaciers (Mer-de-Glace and Argentière glaciers) where a large data set has been collected: ERS, RADARSAT and ENVISAT data, airborne E-SAR data (X, C, L and P bands), airborne photographs, GPS and GPR in-situ measurements…

INTRODUCTION

The monitoring of temperate glacier fast evolution is an important issue for economical and security reasons and as an indicator of the local effects of global climate change. Compared to sparse terrestrial ground measurements, remote sensing is expected to allow regular observations of glacier activity and to provide dense measurements of physical parameters which are necessary to detect significant changes and to constrain glacier flow models. Optical data can be used to derive by correlation techniques high accuracy digital terrain model (DTM) and surface displacement fields, but only during the summer season and below the equilibrium lines where snow disappears for a few months [START_REF] Berthier | Mountain glaciers surface motion derived from satellite optical imagery[END_REF]. Synthetic Aperture Radar data is a complementary source of information since SAR Interferometry (InSAR) or alternative techniques can be applied to measure glacier topography and displacement fields which can reach several decimetres per day in the Alps. Spaceborne data from ERS-1/2 tandem mission have been successfully used to derive velocity fields, mainly during the cold season because of the strong temporal decorrelation in summer [2]. To investigate the potential and processing issues of high resolution or polarimetric SAR data in different bands (X, C, L and P), two airborne SAR campaigns took place over the proposed test-site in October 2006 and February 2007. E-SAR repeat pass interferometric, polarimetric and multi-band data have been acquired thanks to the collaboration between the DLR Microwaves and Radar Institute and the MEGATOR group (http://www.gipsa-lab.inpg.fr/megator). The first results confirmed the high potential of SAR polarimetry in the analysis of the glacier complex surface and sub-surface made of rocks, snow, firm and ice [3]. The repeat pass space-borne data expected from the new generation satellites (ALOS, TerraSAR-X, RadarSAT-2, COSMO-Sky-med) should allow us to investigate different issues such as coherence preservation, atmospheric perturbations, uncertainty assessment…

MONT-BLANC TEST-SITE

The "Mont-Blanc" test site is located in the Alps, near the borders between France, Italy and Switzerland. It goes from "Les Bossons" glacier (45°50'N/6°51'E) in the South, up to the "Trient" glacier (46°59'N/7°01'E) in the North. It covers several glaciers near the Mont-Blanc and the Chamonix valley, including the Mer-de-Glace (45°55'15''N/6°55'45''E, 1600-3000 m) and Argentière (45°56'15'''N/7°00'30''E, 1500-3200 m) glaciers, which are instrumented and observed by many scientific teams. This test site is very interesting for the processing of remote sensed data and the validation of glacier monitoring feasibility, for several reasons:

-the geophysical specificity of the glacier behaviour (continuous motion and temporal changes), -the large datasets acquired in this specific area from different sensors in the last 15 years and expected in the three coming years, -the accessibility of the studied glaciers which makes ground measurements easier for the interpretation and the validation of the extracted information. The high velocity of these Alpine glaciers (in average 270 mm a day in summer on the Mer-de-Glace and Argentière) and their strong topography (in the surroundings of Mont-Blanc) are challenging features for SAR differential interferometry. Moreover, the velocity field is far from being uniform on the whole glacier surface and also varies with the seasons. One important feature of the temperate glacier is that, along the different layers within the glacier cross-section, the temperature profile lies constantly at 0°C. Different ice and snow mixtures, as well as rocks of different sizes and orientations are spread all over the surface of the glacier (cf. Fig. 1). This diversity, observed on the proposed glaciers, leads to different backscattering mechanisms. Their effect will be observed either in PALSAR dual polarization mode and their influence on the coherence level of associated interferograms will be investigated, or in fully polorimetric mode making the use of PolInSAR techniques such as coherence optimization even more interesting. This will allow us to observe the penetration of the L-band in the different parts of the glaciers which should be partly snow covered at the time of the expected fully polarimetric acquisitions (end of spring). These results will be very interesting to correlate with the bottom displacement measured in cavities beneath the Argentière glacier.

TEAM AND PROJECTS

The ALOS ADEN proposal on this test-site has been initiated by a group of 4 laboratories (cf. author affiliations) working in a French national project (MEGATOR, 2004(MEGATOR, -2007) ) dedicated to the development of processing methods for temperate glacier monitoring by optical and SAR remote sensing. This research group has been enlarged with new partners (IETR -Université de Rennes, LGIT -Grenoble/Chambery and Laboratoire de Géologie, ENS Paris) in a new project (EFIDIR, 2008(EFIDIR, -2011)).

The goals of these projects are:

-to develop and combine optical and SAR data processing techniques, -to investigate new sensor potential: high resolution, SAR polarimetry... -to gather experts in information processing and geophysics to improve SAR displacement measurements for different applications: glaciers, earthquake, volcanoes... The partners of MEGATOR project have been developing an experimental data base on the proposed test-site (http://www.demorecherche.univ-savoie.fr/megatorrecherche). It includes space-borne SAR images from ERS1/2 satellites used to derive displacement fields by differential interferometry (cf. Fig. 2), airborne photographs used to derive high resolution DTM by photogrametry (cf. Fig. 3) and airborne SAR data acquired during two E-SAR campaigns together with ground measurements such as GPS velocity profiles (cf. Fig. 4), snow stratigraphic profiles… The E-SAR data include X and C band intereferometric data and L and P band polarimetric interferometric (Pol-InSAR) data. A first PolSAR analysis of an L-band image acquired over Argentière glacier is illustrated in Fig. 5. The PRISM stereoscopic panchro images will allow the computation of accurate up to date DTM's when the snow disappears between May and September. The associated high resolution ortho-images will provide information on surface cover and allow the detection of different features such as rocks, crevasses and glacier limits which are necessary for risk management.

Measured every year, such information will allow comparison and change detection, providing important information on the global evolution with much higher resolution than ground measurements performed only in a few points. These results will also be useful for the processing of the PALSAR data, allowing, for instance, to subtract the topographic fringes in interferograms or to analyze backscattering mechanisms with a priori information on the glacier surface. The PALSAR data are expected i investigate different issues at different seasons:

-single polarisation interferometric images in winter --To n strategy, season to analyse the coherence preservation and, if the level is sufficient, to derive 46-day displacement fields, dual polarisation images in spring and summer seasons for target detection and tracking by bivariate correlation, Quad polarisation interferometric images in winter season to apply coherent/non-coherent decompositions for surface/subsurface feature classification. comply with the ALOS systematic observatio the sought-after PALSAR data have to be acquired in ascending passes. Because of the high topography and the main orientation of the glacier valleys in this area, the glacier visibility is reduced in ascending passes with low incidence angle (cf. Tab. 1). According to the simulations performed from a DTM (cf. Fig. 7), the Quad-Pol mode operated at about 20° would image only a small part of the studied glaciers. Moreover, the Quad-pol data acquired in this area show that the Mont-Blanc test-site is not imaged by the regular tracks (it is between tracks 643 and 644). Accordingly, despite of the lower resolution, specific acquisitions, with about 30° incidence angle, are necessary to image the studied glaciers with ALOS Quad-pol data.

CONCLUSIONS AND PERSPECTIVES

The m ical:

- -understanding of snow/ice L-ba -assessment of L-band interferometric potential and limits on heterogeneous fast moving temperate glaciers. e Mont-blanc area has been proposed as a test-site for "new generation" spaceborne SAR data. Results obtained with ALOS data will be compared with results obtained in X-band (TerraSAR-X and COSMO-SkyMed) and C-band (ENVISAT and Radarsat 2), and benefit from the in-situ measurements and experiments performed in this test-site. For instance, a corner reflector (CR) provided by the DLR and a continuous GPS recording station have been fixed on the ice on the upper part of the Argentière glacier during the February 2007 E-SAR campaign (cf. Fig. 8). Two other non-moving continuous GPS stations have been installed, one near the glacier at 2700 meter height and one down in the Chamonix valley. Despite the difficulties to maintain this equipment on the moving and changing glacier surface, it should allow us to detect the CR with a high accuracy in SAR images whilst its position is known by the GPS data. The CR displacement will be used either as ground truth for performance assessment or as complementary data for solving D-InSAR phase unwrapping ambiguity. The 3 continuous GPS stations will be used to study the electromagnetic path delay variations and to correct the atmospheric perturbations in the SAR based glacier motion measurements. Finally, the analys different sensors, performed with the available prior information on glacier activity, will allow experts to investigate many issues in glacier modeling, evolution tracking and associated risk management.

n Argentière glacier; ENVISAT VV image, February 2007

Figure 1 .

 1 Figure 1. Mer-de-Glace glacier surface, May 2004

Figure 2 .

 2 Figure 2. Mer-de-Glace and Argentière glacier 1-day displacement measured from ERS tandem data, March 96.

Figure 3 .

 3 Figure 3. Mer-de-Glace and Tacul glaciers, high resolution DTM from airborne photographs, August 95.

Figure 4 .

 4 Figure 4. Argentière glacier, 3-day displacement vectors observed by differential GPS, October 2006

Figure 7 .

 7 Figure 7. Glacier visibility for ALOS PALSAR data, La .

  Figure 8. Corner reflector and continuous GPS installed o

Table 1 .

 1 Glacier visibility in PALSAR images

	Incidence angle	20°	30°	45°
	Ascending pass	36 % 79 % 99 %
	Descending pass	76 % 96 % 94 %

Acknowledgement

This research will be co ALOS ADEN proposal # 3638 by researchers of the MEGATOR (http://www.gipsa-lab.inpg.fr/megator) and EFIDIR (http://efidir.fr) projects which are supported by the French National Research Agency.


Jean-Denis Durou
email: durou@irit.fr

Jean-Marc Schlenker
email: schlenker@math.ups-tlse.fr

Shape-from-shading for Surfaces Applicable to Planes

Under the classical assumptions of shape-from-shading, we show that the image of any "applicable surface" (surface applicable to a plane) is also the image of a 1dimensional manifold of applicable surfaces, provided the image contains no singular point. Moreover, we show that the knowledge of a normal in the image suffices to reconstruct the whole shape of the surface, since the problem can be reformulated as an ordinary differential equation w.r.t. the normal, in this case. The usefulness of this theoretical result to document image analysis is straightforward.

Introduction

Among various computer vision techniques for 3Dreconstruction, shape-from-shading (SFS) is sometimes considered as a "toy problem", since its results are unsatisfactory on real images [25,[START_REF] Durou | Numerical Methods for Shape-from-shading: A New Survey with Benchmarks[END_REF]. As for photometric stereo, a main problem of SFS is to get the reflectance properties of the materials which constitute the scene. However, even if this knowledge is available, a second problem with SFS is ill-posedness. In the definition by Hadamard, a problem is well-posed if it has a unique solution and if the solution is stable w.r.t. a small perturbation in the data. Many inverse problems are ill-posed because of instability.

In this paper, we focus on the question of uniqueness of the solution. More precisely, we show that for a class of surfaces, namely the applicable surfaces which contain no singular point, there exist a 1-dimensional manifold of solutions to the SFS problem. Moreover, the knowledge of the normal to the surface at one point is enough to compute the whole 3D-shape. This result is of great interest in the framework of document image analysis.

In Section 2, we recall the classical assumptions of SFS as well as its basic model, and we tackle the question of uniqueness of the solution. In Section 3, we show that the image of an applicable surface without singular point is also that of a 1-dimensional manifold of applicable surfaces. This is the theoretical contribution of the paper. In Section 4, a straightforward application of this result con-sists in reconstructing the 3D-shape of a curved document from its shading. Finally, Section 5 concludes our study and states several perspectives.

Shape-from-shading and Well-posedness

Shape-from-shading Model

We attach to the camera a three-dimensional coordinate system Oxyz, such that Oxy coincides with the image plane H. Let us denote x = (x, y, z) the current point of a scene surface. Under the assumption of orthographic projection on H, the visible part of the scene is a graph z = u(x, y). The SFS problem consists in computing the function u for (x, y) ∈ Ω ⊂ R 2 , from one image only. As is well known [START_REF] Horn | Obtaining Shape from Shading Information[END_REF], it can be modeled by the "image irradiance equation":

R(N(x)) = I(x, y), (1)
where I(x, y) is the greylevel at a point (x, y) in the image (in fact, I(x, y) is the irradiance, but both quantities are usually proportional) and R(N(x)) is the reflectance function, giving the value of the light re-emitted by the surface as a function of its orientation i.e., of the unit normal N(x) to the surface. This normal can easily be expressed as:

N(x) = 1 1 + p(x, y) 2 + q(x, y) 2 (-p(x, y), -q(x, y), 1),
(2) where p = ∂ x u and q = ∂ y u. Assume that there is a unique light source at infinity whose direction is indicated by the unit vector ω = (ω 1 , ω 2 , ω 3 ) ∈ R 3 . Recalling that, for a Lambertian surface of uniform albedo equal to 1, R(N(x)) = φ ω, N(x) , where φ is the "density" of the illumination, Eq. ( 1) can be written, using (2):

I(x, y) φ 1 + ∇u(x, y) 2 + (ω 1 , ω 2 ), ∇u(x, y) -ω 3 = 0, (3)
which is a first order non-linear partial differential equation (PDE) of the Hamilton-Jacobi type. Points (x, y) ∈ Ω such that I(x, y) = φ correspond to the particular situation where ω and N(x) point in the same direction. These points are usually called the "singular points" [START_REF] Horn | Obtaining Shape from Shading Information[END_REF].

Let us consider the equation which appears in most of the papers and which corresponds to vertical light source at infinity i.e., ω = (0, 0, 1). Then (3) becomes the "eikonal equation" [5]:

p(x, y) 2 + q(x, y) 2 = f (x, y), (4)
with f (x, y) = φ 2 /I(x, y) 2 -1. In this particular case, the singular points are those where where p(x, y) = 0 and q(x, y) = 0.

Let us mention that the eikonal equation is not only the SFS model for Lambertian materials, since it also holds for any radially-symmetric reflectance [START_REF] Oliensis | Uniqueness in Shape from Shading[END_REF]. Finally, more sophisticated models have appeared in the last few years, but we concentrate in this paper on the simple model of the eikonal equation.

Uniqueness of the Solution

The number of solutions to the SFS problem has been addressed in various situations. In [START_REF] Belhumeur | The Bas-Relief Ambiguity[END_REF], Belhumeur et al. determine the set of images of an object under all possible lighting conditions, a problem currently referred to as the "bas-relief ambiguity". The existence of impossible shaded images is proved in [START_REF] Oliensis | Uniqueness in Shape from Shading[END_REF]4,[START_REF] Horn | Impossible Shaded Images[END_REF]. However, our concern in this paper is uniqueness.

First of all, it is obvious that a translation of the surface in the direction of Oz does not modify the first member of Eq. (4). Moreover, Eq. ( 4) is also invariant by symmetry w.r.t. the image plane H. This is known as the "concave/convex ambiguity". In his pionneering work [START_REF] Horn | Obtaining Shape from Shading Information[END_REF], Horn shows that the presence of a singular point in the image much reduces the number of possible ambiguities but, whereas N is completely known at a singular point, only its vertical component N z = N, e z is known at other points. For this reason, SFS is traditionally considered as an ill-posed problem.

Bruss uses a Taylor series expansion at a singular point in order to characterize the solutions [5], and states the first proof of uniqueness, which requires the presence of a closed "occluding contour" i.e., of a contour on which I(x, y) = 0 or, equivalently, p(x, y) 2 + q(x, y) 2 = +∞. On such a contour, as well as at a singular point, the normal is computable [START_REF] Ikeuchi | Numerical Shape from Shading and Occluding Boundaries[END_REF], which is of great importance to constrain the problem. Blake et al. state another uniqueness result [2] in the case where u is known on a closed contour (Dirichlet boundary condition).

Oliensis concentrates on the occluding contours. He claims that, in many situations, regularization techniques (see [START_REF] Horn | The Variational Approach to Shape From Shading[END_REF]) should be avoided in SFS, since the problem is well-posed and regularization may distort the solution. He states a uniqueness theorem which is more generic than Bruss' one [START_REF] Oliensis | Uniqueness in Shape from Shading[END_REF], and partially extends it to the case of any direction of illumination [START_REF] Oliensis | Shape from Shading as a Partially Well-Constrained Problem[END_REF], a situation which is also explored by Saxberg in [23].

Brooks et al. exhibit a counterexample to Bruss' theorem in [3] and address the problem of radially-symmetric eikonal equations. Finally, most of these contributions, which generally use the "characteristic strips method" (see [START_REF] Horn | Obtaining Shape from Shading Information[END_REF] and Section 3.3.1), are updated by Kozera in [START_REF] Kozera | Uniqueness in Shape from Shading Revisited[END_REF].

The problem of uniqueness in SFS has also been addressed in the framework of "viscosity solutions", which are solutions in a weak sense [START_REF] Lions | Shape-from-Shading, viscosity solutions and edges[END_REF]. Of course, SFS should be less-constrained if u is a viscosity solution than if it is differentiable. Indeed, uniqueness of the viscosity solution requires, most of the time, Dirichlet boundary conditions as well as the knowledge of u at each singular point [START_REF] Lions | Shape-from-Shading, viscosity solutions and edges[END_REF], since precisely a singular point improves the number of viscosity solutions! Therefore, the concept of "maximal viscosity solution" [START_REF] Falcone | An algorithm for the global solution of the Shape-from-Shading model[END_REF]6] eliminates the need for a priori additional information. More recently, Prados et al. prove an existence and uniqueness result for the perspective SFS problem coupled to state constraint boundary condition, in the case where the light source is located at the center of projection [22,[START_REF] Prados | A generic and provably convergent Shape-From-Shading Method for Orthographic and Pinhole Cameras[END_REF].

Invisible Continuous Deformations

A somewhat different problem is the search for invisible continuous deformations i.e., manifolds of surfaces which are solutions to the same eikonal equation, and which can therefore be continuously distorted without any change in the image. More precisely, an invisible continuous deformation is a family {u c } c∈R of functions Ω → R, depending on a real parameter c, which all are solutions to Eq. (4). Differentiating (4) w.r.t. c, it comes:

∂ c ∇u c (x, y) 2 = 0, (5)
or:

∇u c (x, y), ∇∂ c u c (x, y) = 0. (6)
Nevertheless, this equation does not prove that invisible continuous deformations exist.

Starting from a solution u 0 , a construction of an invisible deformation could proceed as follows. Consider the family of curves of equations u c (•) = h, for all h ∈ R, and the family of curves of equations ∂ c u c (•) = h, for all h ∈ R. Assume that these two families are orthogonal, as required by Eq. ( 6). This ensures that each u c is a solution to Eq. (4). Of course, in order to make a proof from this heuristics, one would still have to show that the differential equation (6) could indeed be integrated, especially at a singular point where all the lines ∂ c u c (•) = h, h ∈ R, converge. Furthermore, all the solutions which are built with this procedure could represent surfaces without real deformation. Consider for instance f (x, y) = 1 and u 0 (x, y) = x. From (6), one can get u c (x, y) = x cos c -y sin c, but u c is the composition of u 0 and of a rotation.

Nevertheless, although this might seem to be almost impossible, the existence of an invisible continuous deformation in presence of a singular point is proved in [START_REF] Durou | Ambiguous Shape from Shading with Critical Points[END_REF]. In the next section, we will focus on the existence of invisible continuous deformations of applicable surfaces.

Invisible Continuous Deformations of Applicable Surfaces

Basic Equations

We deal hereafter with some classical tools of differential geometry. We refer the interested reader, e.g., to [START_REF] Brooks | Impossible and Ambiguous Shading Patterns[END_REF]7].

Let S ⊂ R 3 be an applicable surface. It is well known that S is ruled. As a simplification we suppose that through each point x ∈ S goes exactly one line D x contained in S (there might also be flat pieces in S). The lines D x project on H to lines ∆ x which foliate Ω. The surface S can be locally parameterized as follows:

x(s, t) = t v(s) + a(s), (7)
where v(s) is a unit vector (along D x ). We can suppose additionally that a (s) is orthogonal to v(s). Let w(s) be the orthogonal projection of v(s) on H.

Under the simple SFS model presented in Section 2.1, the vertical component N z of N follows from the shading. Let II be the second fundamental form of S.

Remark At all x ∈ S, the kernel of II is the set of vectors parallel to D x . It follows that N is constant along the lines D x . Thus, N is a function of s only, and ∂ s N is orthogonal to v at each point.

To sum up, at each point:

• N (s) is orthogonal to N(s) and to v(s), so that it can be expressed in terms of N z (s):

N (s) = N z (s) N(s) × v(s), e z N(s) × v(s). (8)
Note that the denominator in the second member of Eq. ( 8) vanishes if N is parallel to e z i.e., at a singular point.

• v(s) is orthogonal to N(s) and its projection w(s) to H is known, so that it can be recovered from N(s) and w(s):

v(s) = iw(s) × N(s) iw(s) × N(s) . (9)
Here iw(s) means that we consider the complex structure on H identified with R 2 and with C.

Deformations with Constant Shading

Our main claim is that many applicable surfaces can be deformed without changing the shading.

Theorem 1 Suppose that S is an applicable surface which is a graph over Ω, that Ω is bounded, and that N is nowhere vertical. Then S has a 1-parameter deformation as an applicable surface such that N z , considered as a function over Ω, is constant.

The proof rests on a converse of the remarks in section 3.1.

Remark Suppose given w, N z and let v, N be functions of s satisfying equations ( 8) and [START_REF] Durou | Ambiguous Shape from Shading with Critical Points[END_REF]. Then there exists a (unique) applicable surface S, going through a given point m ∈ R 3 over Ω, such that the projection of v is w and that the z-component of N is N z .

The condition that the surface goes through a given point m is necessary to choose one solution among a 1-parameter family invariant under vertical translation.

Proof of Theorem 1

The starting point is a classical result from differential geometry: given two vector fields u, v in R 3 , linearly independent at each point, there exists a foliation of R 3 by surfaces to which u and v are everywhere tangent if and only if [u, v] ∈ vect(u, v) at each point, where

[•, •] is the Lie bracket. Let v ⊥ = N × v.
Then v ⊥ is constant along the integral lines of v i.e.:

D v v ⊥ = 0, ( 10
)
where D is the flat connection of R 3 (directional derivative). It is possible to suppose (in the parametrization of S described above) that a (s) is orthogonal to v. It follows that a (s) is parallel to v ⊥ . Consider v, N as two vector fields defined on R 3 , invariant by vertical translation. Note now that N (s) is orthogonal to v(s) by Eq. ( 8). Since N(s), v(s) = 0:

N (s), v(s) + N(s), v (s) = 0, (11)
so that v (s) is orthogonal to N(s). Thus, the Lie bracket:

[v, v ⊥ ] = D v v ⊥ -D v ⊥ v = -D v ⊥ v (12)
is orthogonal to N. This means that the distribution of 2planes spanned by v and v ⊥ is integrable. Integrating it yields a foliation of the union of vertical lines intersecting Ω by surfaces. This foliation is invariant under vertical translation.

Moreover the definition shows that the vector field w lifts to those surfaces as vector fields v so that the integral curves of v are straight lines, and the vector field N considered in the construction is orthogonal to those surfaces. This implies that the second fundamental form of those surfaces has non-trivial kernel -the kernel is generated by v -and so the surfaces are applicable, by the Gauss equation.

Examples

Let us consider the following eikonal equation:

p(x, y) 2 + q(x, y) 2 = 1 + x 2 . ( 13
)
The lines ∆ x are parallel to Oy. Let us solve Eq. ( 13) using a classical method of resolution of the first order non-linear PDEs.

Characteristic Strips Method

Since p 2 + q 2 > 0 everywhere, there is no singular point in the image (cf. Fig. 1) associated to Eq. ( 13). The five ordinary differential equations of the so-called "characteristic strips method" [START_REF] Horn | Obtaining Shape from Shading Information[END_REF] are, in our case:

ẋ = 2 p, (14a)
ẏ = 2 q, (14b)
ż = 2 (1 + x 2 ), (14c)
ṗ = 2 x, (14d)
q = 0, (14e)
where the dot denotes differentiation w.r.t. a real parameter s which could be, for example, an arc-length. In fact, the five unknowns (x, y, z, p, q) depend on another parameter than s, but Eqs. (14a)-(14e) suppose that this second parameter t is constant along each characteristic strip. Integrating this system, we can reconstruct the whole shape from an "initial curve" on which the 5-tuple (x, y, z, p, q) is known. Moreover, at any "starting point" on the initial curve, the vector ( ẋ, ẏ), which by definition is tangent to the current characteristic strip, must be non-tangent to the initial curve. Finally, this last vector must be everywhere nonnull, which here is guaranteed by the form (13) of p 2 + q 2 . Let us use, for instance, the axis Oy as initial curve. We can thus define t as the y-component of the starting point. For a given starting point (x, y) = (0, t), the 2-tuple (p, q) cannot be uniquely determined from Eq. ( 13). Let us choose (p(0, t), q(0, t)) = (f (t), g(t)), with f (t) 2 + g(t) 2 = 1 and f (t) = 0, so that (p, q) is not parallel to Oy. On the one hand, Eqs. ( 14a) and (14d) give:

∂ 2 s x(s, t) -4 x(s, t) = 0, ( 15
)
with the initial conditions x(0, t) = 0 and ∂ s x(0, t) = 2 f (t). Then:

x(s, t) = f (t) sinh(2 s). ( 16
)
On the other hand, Eqs. ( 14b) and (14e) give:

∂ 2 s y(s, t) = 0, (17)
with the initial conditions y(0, t) = t and ∂ s y(0, t) = 2 g(t). Then: y(s, t) = 2 s g(t) + t.

Finally, Eqs. (14c) and ( 16) give:

∂ s z(s, t) = 2 -f (t) 2 + f (t) 2 cosh(4 s), (19)
with the following initial condition:

z(0, t) = t 0 g(u) du + z(0, 0). ( 20
)
Because z has to be fixed at some point, we set z(0, 0) = 0. Therefore, from ( 19) and (20):

z(s, t) = 2 -f (t) 2 s + f (t) 2 sinh(4 s) 4 + t 0 g(u) du. ( 21
)

Applicable Solutions

For any applicable solution to Eq. ( 13), it has been proved in Section 3.1 that the normal is uniform along Oy, and so are f and g as well. Knowing that f (t) 2 + g(t) 2 = 1, there exists a real θ = 0 such that f (t) = cos θ and g(t) = sin θ.

Eqs. ( 16), ( 18) and ( 21) can thus be rewritten:

x(s, t) = cos θ sinh(2 s), (22a)
y(s, t) = 2 s sin θ + t, (22b)
z(s, t) = (2 -cos 2 θ) s + cos 2 θ sinh(4 s) 4 + t sin θ. ( 22c
)
This is a parametric representation of the solution S θ , which is parameterized by θ. Its explicit equation is straightforward, since s and t can be expressed from (22a) and (22b):

s = 1 2 argsh x(s, t) cos θ , ( 23a
) t = y(s, t) -sin θ argsh x(s, t) cos θ . ( 23b
)
Replacing s and t by these expressions in (22c), the explicit equation of S θ is:

z = cos 2 θ 2 argsh x cos θ + x 2
x 2 + cos 2 θ+y sin θ. [START_REF] Brooks | Impossible and Ambiguous Shading Patterns[END_REF] from which it arises that all surfaces S θ are cylinders, whose generatrices are parallel to the plane Oyz. Therefore, a continuous modification of θ in ] -π/2, π/2[ gives rise to an invisible continuous deformation of the surface S θ . In Fig. 2, three cylinders among this 1-dimensional manifold are shown, for s ∈ [-0.9, 0.9] and t ∈ [-5, 5].

Other Solutions

Of course, other solutions to Eq. ( 13) exist, that are not applicable. For instance, for a given real T = 0, let us choose f (t) = cos(t/T ) and g(t) = sin(t/T ), with t/T ∈ ] -π/2, π/2[. Eqs. ( 16), ( 18) and ( 21) can be rewritten:

x(s, t) = cos t T sinh(2 s), (25a)
y(s, t) = 2 s sin t T + t, (25b) z(s
, t) = (2 -cos 2 t T )s + cos 2 t T sinh(4 s) 4 + T (1 -cos t T ). (25c)
This is a parametric definition of non-applicable solutions which we call T T . In Fig. 3, three of these non-applicable solutions are shown, for s ∈ [-0.9, 0.9] and t/T ∈ [-π/4, π/4].

On the one hand, the parametric definition (25a)-(25b)-(25c) is not as easy to transforme in an explicit equation as for the applicable solutions. On the other hand, if the surfaces T T form a 1-dimensional manifold of solutions, parameterized by T , they do not constitute all the nonapplicable solutions to Eq. (13).

Presence of an Occluding Contour

Now let us consider the following eikonal equation:

p(x, y) 2 + q(x, y) 2 = 1 x 2 , for x = 0. ( 26
)
This case is interesting because p 2 + q 2 → +∞ when

x → 0, which means that each solution to Eq. ( 26) has Ox as occluding contour. Therefore, N → ±e x when x → 0. Nevertheless, this knowledge on N does not constrain the applicable solutions. As a proof, here is the explicit equation of a 1-dimensional manifold of applicable solutions to Eq. ( 26), which are parameterized by a parameter c:

z = 1 -c 2 x 2 -argch 1 c x + c y. (27)
In Eq. ( 27), however, one has to impose 1/(c x) ≥ 1, since

argch is defined on [1, +∞[.

Application to Document Image Analysis

Theorem 1 can be straightforwardly applied to document image analysis and, for example, to the simulation of document flattening.

Simulation of Document Flattening

The digitization of documents currently knows an increasing popularity, because of the expansion of Internet browsing. The traditional process, which uses a flatbed scanner, is satisfactory for flat documents, but is unsuitable for curved documents like for example a thick book, since some defects will appear in the digitized image. Several specific systems have been designed, but such systems are sometimes intrusive with regard to the documents and, before all, they cannot be referred to as consumer equipments. An alternative consists in simulating the flattening of curved documents i.e., in correcting the defects of images provided by a flatbed scanner or a digital camera. In order to successfully simulate the document flattening, it is necessary to compute its surface shape.

Among many papers dealing with this problem, two only deal with the most general case of applicable surfaces [20,[START_REF] Gumerov | 3D structure recovery and unwarping of surfaces applicable to planes[END_REF] and, moreover, none of both relies on SFS. Our result could potentially lead us to a new algorithm of flattening simulation.

A New Method of 3D-reconstruction

As a concrete application of Theorem 1, we can design a new method of 3D-reconstruction from one shaded image of a document. Given the normal along an "initial curve", which should be one of the straight lines ∆ x , 3Dreconstruction could of course be carried out through using the characteristic strips method (cf. Section 3.3.1), but the implementation would be tedious: new characteristic strips would have to be created when they separate too much, or deleted when they approach each other too closely [START_REF] Horn | Obtaining Shape from Shading Information[END_REF].

We propose a new method of 3D-reconstruction of a document shape that is much easier to implement, because the choice of an integration path is free. Once the normal is fixed at some point in the image, and knowing that the normal is uniform on each line ∆ x and that these lines foliate Ω, the integration of Eqs. ( 8) and ( 9) can be carried out along any path in Ω, provided this path crosses all the ∆ x . After that, the integration of the computed normal field is very simple, since Theorem 1 tells us that the normal field of any applicable surface is integrable.

Example

We print a regular black grid on a square sheet (cf. Fig. 4), and we simulate the 3D-warping of this sheet on the three cylinders of Fig. 2, as shown in Fig. 5. It could seem that the three shapes of Fig. 5 are not those of Fig. 2. It happens that, among the three surfaces viewed in Fig. 2, only that on the left would be rectangular after flattening. The others would not, and this is why our appreciation of these shapes in erroneous. The images of the three shapes of Fig. 5, seen from above under the classical assumptions of SFS, are shown in Fig. 6. In the images of Fig. 6, the sheet contour, as well as the deformation of the grid, tell us much about the shape of the sheet, but this information may lead us to an erroneous interpretation of the shape, as in the case of Fig. 2! Using only the shading information in the images of Fig. 6, we cannot avoid the ambiguity expressed in Theorem 1. Moreover, the non-uniform albedo of these scenes requires a pretreatment in which the inked points, whose albedo is equal to 0, have to be separated from the non-inked points, whose albedo is equal to 1. This could be done, e.g., through adaptive thresholding but, since this is out of scope, we use in fact the non-textured image of Fig. 1 instead of the images of Fig. 6. The simplest integration path is Ox. Three profiles obtained from three possible normals at (0, 0) are shown in Fig. 7. This simple example shows, if necessary, that SFS remains ill-posed in many situations since, even in the framework of the applicable surfaces, one cannot decide between the three shapes of Fig. 7 and many others.

Conclusion and Perspectives

In this paper, we show that SFS is ill-posed, under the simple model of the eikonal equation, even if the scene surface is supposed to be applicable to a plane. If no singular point is visible, the applicable solutions form a 1dimensional manifold. The knowledge of one normal thus suffices to render the problem well-posed. A straightforward application of this theoretical result is document image analysis, and especially the simulation of document flattening.

Many prospects remain open. First of all, we should examine how much a singular point constrains the problem. We must extend Theorem 1 to applicable surfaces with flat pieces. It is also necessary to address more realistic SFS models: perspective projection, as well as various directions of illumination, should be explored. It could happen that another model better constrains the problem (see [22]).

A remaining problem of our new method of 3Dreconstruction is that the knowledge of one normal is required, at some point that can be neither a singular point nor a point lying on an occluding contour. In order to render the 3D-reconstruction well-posed, we now address the information given by the contour, as mentioned in Section 4.3. Moreover, we already proved the following interesting result (which is not detailed in this paper, due to lack of space): knowing the image of a straight line on the flattened surface (which is a "geodesic segment" on S) can be used to recover S without ambiguity from its shading information.

Figure 1 .

 1 Figure 1. Image associated to the eikonal equation (13).

Figure 2 .

 2 Figure 2. Three cylinders giving the same image as Fig. 1: S0 (left), S π/6 (middle); S π/3 (right).

Figure 3 .

 3 Figure 3. Three non-applicable surfaces giving the same image as Fig. 1: T5 (left), T7.5 (middle); T10 (right).

Figure 4 .

 4 Figure 4. Regular black grid printed on a square sheet.

Figure 5 .

 5 Figure 5. Warping of the sheet of Fig. 4 on the three cylinders of Fig. 2: S0 (left), S π/6 (middle); S π/3 (right).

Figure 6 .

 6 Figure 6. Images of the three scenes of Fig. 5, seen from above under the classical assumptions of SFS.

Figure 7 .

 7 Figure7. Three profiles reconstructed from the image of Fig.1, using the normals at (0, 0) of: S0 (left), S π/6 (middle); S π/3 (right).