Adel Hamdi
email: adel.hamdi@insa-rouen.fr

Abderrahim Jardani
email: abderrahim.jardani@univ-rouen.fr

Reconstruction of unknown storativity and transmissivity functions in 2D groundwater equations

The paper deals with the identification of unknown storativity and transmissivity distributions within a 2D confined aquifer using pumping tests. We introduce a change of variables that transforms the groundwater equation into a diffusionreaction one, where the diffusion term is the fraction transmissivity/storativity whereas the reaction term yields the right hand side of a second order nonlinear partial differential equation satisfied by the unknown storativity function. Using records of the drawdown at some measuring wells within the monitored aquifer, we establish identifiability results on the introduced diffusion and reaction terms as well as on the storativity values at the employed wells. We develop an identification approach that starts by determining the auxiliary diffusion and reaction variables. Afterwards, this approach uses an assumption related to the incompressibility of water to develop a local determination procedure of the unknown storativity function. Besides, based on the interpolation of its values at the employed wells, a global determination procedure of this function is also developed. The unknown transmissivity is then determined by the product of the identified storativity and fraction transmissivity/storativity functions. Some numerical experiments are presented.

Introduction

Managing effectively groundwater resources requires the knowledge of some hydraulic properties defining the nature of the involved aquifers. For instance, among the main required properties we quote the following, see [START_REF] Meier | A numerical study on the relationship between transmissivity and specific capacity in heterogeneous aquifers[END_REF][START_REF] Todd | Groundwater Hydrology. 2nd Edn[END_REF] Storativity: That is the volume of water released from storage with respect to the change in head (water level) and to the aquifer's surface area. Transmissivity: It represents the aquifer's ability to transmit groundwater throughout its entire saturated thickness. Since those properties affect significantly groundwater movement and storage as well as solute transport, in the literature numerous studies have been devoted to estimate such properties. Those studies are mainly based on the so-called pumping test analysis [START_REF] Bulter | Pumping tests for aquifer evaluation, time for change?[END_REF][START_REF] Hoeksema | An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling[END_REF][START_REF] Meier | A numerical study on the relationship between transmissivity and specific capacity in heterogeneous aquifers[END_REF][START_REF] Yeh | Time to change the way we collect and analyze data for aquifer characterization[END_REF] which consists of analysing data, measured in some surrounding well tests, that represents the aquifer's response to a hydraulic forcing term introduced through one or multi-pumping wells. In practice, accurate estimations of the hydraulic properties lead to employ more appropriate actions in large spectrum of applications that go from groundwater exploration to wastedisposal evaluation as well as to the determination of efficiency and productivity of a well for groundwater extraction. For example, estimating transmissivity could help well field managers to design more energy-efficient pumping schemes since Sterrett reported in [START_REF] Sterrett | Groundwater and wells, 3rd Eddition[END_REF] that the energy requirement for pumping is directly proportional to the hydraulic lift. In addition, an accurate estimation of storativity is important for quantifying groundwater availability in order to satisfy drinking water demand, for instance, municipal wells in Colorado extract over 100 million cubic meters of groundwater per year.

In the literature, the identification of aquifers hydraulic properties has been initiated by Theis in [START_REF] Theis | The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage[END_REF] who used the so-called type-curve matching method to estimate aquifers parameters. Then, Cooper and Jacob employed in [START_REF] Cooper | A generalized graphical method for evaluating formation constants and summarizing well-field history[END_REF] the straight-line method to determine hydraulic properties in groundwater equations. Those techniques, called also graphical approach, are based on matching graphically the recorded data at the pumping tests to the simulations of several analytical models depending on the type of aquifers and the hydraulic conditions. In the last few decades when the use of computer became widely available, we have seen emerging several new approaches extending the graphical approach to solve applications with big data fitting and to estimate wider range of hydraulic properties as well as to explore larger class of aquifers. Those approches employ different techniques to identify the underlined hydraulic properties, for instance, stochastic techniques have been employed in [START_REF] Yeh | Stochastic Modelling of Groundwater Flow and Solute Transport in Aquifers[END_REF] and genetic algorithms in [START_REF] Abdel-Aziz | Parameter estimation of pumping test data using genetic algorithm[END_REF] whereas geostatistical inversion of data by using the Bayesian approach have been used in [START_REF] Beckie | What does a slug test measure: an investigation of instrument response and the effects of heterogeneity[END_REF][START_REF] Castagna | A Bayesianapproach for inversion of hydraulic tomographic data[END_REF][START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF][START_REF] Hoeksema | An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling[END_REF] and probabilistic estiamtions using time series model in [START_REF] Shapoori | Estimating aquifer properties using groundwater hydrgraph modelling[END_REF]. Besides, there is a direct identification approach developed in [START_REF] Nelson | In-place measurement of permeability in heterogeneous media: Experimental and computational considerations[END_REF][START_REF] Richter | An inverse problem for the steady state diffusion equation[END_REF] that consists of considering the transient groundwater equation along the streamlines associated to the gradient of the drawdown as a first order ordinary differential equation of the unknown transmissivity whereas the storativity and the drawdown are both supposed to be known. Provided an initial value of the transmissivity is given for each streamline, this approach transforms the identification of the transmissivity into solving a Cauchy problem. However, the knowledge of the streamlines and of an initial transmissivity value for each streamline are difficult to achieve in the implementation of a real case. In [START_REF] Rogelio | The differential system method for the identification of transmissivity and storativity[END_REF], the authors introduced the so-called Differential System (DS) method that, given the drawdown for three different flows, solves the Cauchy problem in a way that doesn't require anymore the knowledge of streamlines, an initial transmissivity value is needed in only one point and no a priori knowledge of the storativity is required. Nevertheless, some remarks have been made by Beckie in [START_REF] Beckie | What does a slug test measure: an investigation of instrument response and the effects of heterogeneity[END_REF] that the estimation of storativity is a very sensitive problem since it is influenced by the estimated transmissivity. Meier in [START_REF] Meier | A numerical study on the relationship between transmissivity and specific capacity in heterogeneous aquifers[END_REF] reported that the estimation of those properties in heterogeneous media depends on the measurement locations.

In the present study, we focus on identifying unknown storativity and transmissivity distributions within a 2D confined aquifer using pumping tests. Based on the analysis and optimisation of the groundwater partial differential equation governing the drawdown in the considered aquifer, we develop an identification approach that uncouples the determination of the unknown storativity function from the identification of the unknown trans-missivity. Moreover, the developed approach establishes conditions on the used pumping source as well as on the number and the locations of the employed measuring wells to ensure uniqueness of the involved unknown functions. The paper is organized as follows: Section 2 is devoted to introduce the problem statement and to establish some technical results for later use. In section 3, we study the identifiability of the occurring unknown functions. Section 4, is reserved to develop the identification approach that leads to determine the unknown storativity and transmissivity functions. Some numerical experiments on a variant of groundwater equations are presented in section 5.

Problem statement and technical results

Let T > 0 be a finite final monitoring time and Ω be a bounded and connected open subset of IR 2 with Lipschitz boundary ∂Ω. According to [START_REF] Anderson | Applied Groundwater Modeling[END_REF][START_REF] Bear | Hydraulics of Groundwater[END_REF][START_REF] Cooper | A generalized graphical method for evaluating formation constants and summarizing well-field history[END_REF][START_REF] Todd | Groundwater Hydrology. 2nd Edn[END_REF][START_REF] Theis | The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage[END_REF], it follows that due to the similarity between groundwater flow and heat conduction, the hydraulic head (water level or also called drawdown), denoted here by u, in a confined aquifer Ω subject to an external hydraulic pumping source f is governed by the following system:

         L[u](x 1 , x 2 , t) = f (x 1 , x 2 , t) in Ω × (0, T ) u(•, 0) = 0 in Ω P ∇u • ν = 0 on ∂Ω × (0, T ) (1)
where ν is the unit outward vector normal to ∂Ω and L is the second order linear partial differential operator defined by

L[u](x 1 , x 2 , t) := S(x 1 , x 2 )∂ t u(x 1 , x 2 , t) -div P (x 1 , x 2 )∇u(x 1 , x 2 , t) (2)
In (2), S and P designate the storativity and the transmissivity functions defining the nature of the understudy aquifer Ω. In the remainder of this paper, we consider that S and P are two differentiable functions that belong to the following admissible set:

A := 0 < P ∈ W 1,∞ (Ω), 0 < S ∈ W 2,∞ (Ω) and P ∇S • ν = 0 on ∂Ω (3)
We employ a single pumping well that reaches the aquifer at the point a ∈ Ω through which a hydraulic time-dependent forcing function ∈ L 2 (0, T ) is pumped. Therefore, the external time-dependent pumping source f involved in (1) is defined by

f (x 1 , x 2 , t) = (t)δ a (x 1 , x 2 ), for all (x 1 , x 2 , t) ∈ Ω × (0, T ) (4)
where δ a denotes the Dirac mass at the pumping position a. Thus, given S and P elements of the set A together with a ∈ Ω and ∈ L 2 (0, T ) defining f in (4), the forward problem (1)-( 4) admits a unique solution u that belongs to the functional space, see [START_REF] Lions | Pointwise control for Distributed Systems in Control and Estimation in distributed ParametersSystems[END_REF][START_REF] Schwartz | Théorie des distributions[END_REF]:

L 2 0, T ; L 2 (Ω) ∩ C 0 0, T ; H -1 (Ω) (5)
Moreover, the state u is sufficiently regular in Ω \ {a}. Therefore, given the positions of some measuring wells b i=1,...,I ∈ Ω \ {a}, we define the observation operator as follows:

M [S, P ] := u(b i , t) for all t ∈ (0, T ), for i = 1, . . . , I

Later on in this paper, the number I ∈ IN * of measuring wells and their positions b i=1,...,I with respect to the pumping location a ∈ Ω will be further discussed.

The nonlinear inverse problem with which we are concerned here consists of: Given time records d i (t), ∀t ∈ (0, T ) of the state u taken at the measuring wells b i=1,...,I , determine the two unknown functions S and P of A involved in the problem (1)-( 4) that yield M [S, P ] = d i (t) for all t ∈ (0, T ), for i = 1, . . . , I

For later use, to each two differentiable functions S and P of the admissible set A, we associate the intermediate variables ψ, Ψ and ρ defined as follows:

ψ = P S , Ψ = 1 S ∇S and ρ = 1 4 ψ Ψ 2 2 + 1 2 div ψΨ (8)
where • 2 denotes the euclidean norm. Notice that in the case when the two functions ψ and ρ are both known, it follows from (8) that the unknown storativity function S solves in Ω the following second order nonlinear partial differential equation:

ψ 4 Ψ 2 2 + 1 2 ∇ψ • Ψ + ψ 2 div Ψ = ρ ⇔ 1 S ∆S + 1 ψ ∇ψ - 1 2S ∇S • ∇S = 2ρ ψ (9)
Besides, using the auxiliary variables ψ and ρ introduced in (8), we consider the eigenvalue problem: For all n ∈ IN , find µ n and ξ n that solve the system:

   -div ψ∇ξ n + ρξ n = µ n ξ n in Ω ψ∇ξ n • ν = 0 on ∂Ω (10)
For the simplicity of our notations, in the remainder we denote by {ξ n } the set of normalized eigenfunctions solutions of [START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF]. Then, according to [START_REF] Schwartz | Théorie des distributions[END_REF], we have:

Theorem 2.1 (See [START_REF] Schwartz | Théorie des distributions[END_REF]) Let Ω be a bounded open subset of IR 2 with Lipschitz boundary ∂Ω. The normalized eigenfunctions {ξ n } solutions of the system (10) form a complete orthonormal family of L 2 (Ω) and their associated eigenvalues (µ n ) form an increasing sequence of real numbers that tends to infinity.

In the spectral Neumann decomposition Theorem 2.1, the Lipschitz condition on the boundary ∂Ω is required for the compactness of H 1 (Ω) imbedding in L 2 (Ω). Moreover, as far as the first eigenpair of the eigenvalue problem [START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF] is concerned, we quote the following properties: For more details, see for instance [START_REF] Berestycki | Elliptic eigenvalue problems with large drift and applications to nonlinear propagation phenomena[END_REF][START_REF] Chen | Principal eigenvalue and eigenfunction of an elliptic operator with large advection and its application to a competition model[END_REF].

Remark 2.2

The principal eigenvalue µ 0 of the problem (10) is unique i.e., µ 0 < µ n , ∀n ∈ IN * . In addition, for the case when ρ is a real number, the first eigenpair solution of the eigenvalue problem [START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF] is defined by: µ 0 = ρ and ξ 0 (x 1 , x 2 ) = 1/ S(Ω) for all (x 1 , x 2 ) ∈ Ω, where S(Ω) denotes the surface area of the domain Ω.

Furthermore, we introduce what we will refer to in the remainder as strategic position.

Definition 2.3 Let {ζ n } be a complete orthogonal family of continuous functions in L 2 (Ω). We say that

(x 1 , x2 ) ∈ Ω is strategic with respect to {ζ n } if ζ n (x 1 , x2 ) = 0, ∀n.
The notion of strategic position in the sense of Definition 2.3 is well known in the literature. Indeed, this notion has been introduced by El Jai and Pritchard in [START_REF] Jai | Capteurs et Actionneurs dans l'Analyse des systèmes distribués[END_REF] and used by many other authors, for example, in [START_REF] Badia | On an inverse source problem for the heat equation: Application to a pollution detection problem[END_REF][START_REF] Hamdi | Detection and identification of multiple unknown time-dependent point sources occurring in 1D evolutionj transport equations[END_REF]. For later use, to each measuring position b i ∈ Ω \ {a} we associate an impulse response G b i that solves:

   -div ψ∇G b i + ρG b i = δ b i (x 1 , x 2 ) in Ω ψ∇G b i • ν = 0 on ∂Ω (11)
Let G a be the solution of the system (11) with a instead of b i . From multiplying the first equation in [START_REF] Cooper | A generalized graphical method for evaluating formation constants and summarizing well-field history[END_REF] by G a and integrating by parts over Ω using Green's formula, we obtain

G a (b i ) = Ω ψ∇G a ∇G b i + Ω ρG a G b i = -div ψ∇G a + ρG a , G b i = G b i (a) (12)
where , represents the product in the distribution sense. The result in (12) yields a symmetric property of the impulse response solution of the system [START_REF] Cooper | A generalized graphical method for evaluating formation constants and summarizing well-field history[END_REF]. Moreover, using the complete orthonormal family {ξ n }, the solution of ( 11) is given by

G b i (x 1 , x 2 ) = n≥0 ξ n (b i ) µ n ξ n (x 1 , x 2 ), ∀(x 1 , x 2 ) ∈ Ω (13)
Besides, we employ the change of variables:

U (x 1 , x 2 , t) = S(x 1 , x 2 )u(x 1 , x 2 , t) in Ω × (0, T ). That leads to ∇u = S 1 2
1 S ∇U -U 2S 2 ∇S . Afterwards, using the intermediate variables ψ, Ψ and ρ introduced in (8), we get div P ∇u = div S 1 2 ψ∇U -

1 2 ψΨU = S 1 2 div ψ∇U -1 4 ψ Ψ 2 2 + 1 2 div ψΨ U = S 1 2 div ψ∇U -ρU (14)
Since in view of (3) we have P ∇S • ν = 0 on ∂Ω, it follows that P ∇u • ν = S 1 2 ψ∇U • ν on ∂Ω. Hence, the problem (1)-( 2) is equivalent to the following system:

         ∂ t U -div ψ∇U + ρU = S -1 2 f in Ω × (0, T ) U (•, 0) = 0 in Ω ψ∇U • ν = 0 on ∂Ω × (0, T ) (15)
Furthermore, as the set of normalized eigenfunctions {ξ n } solutions of the eigenvalue problem (10) forms a complete orthonormal family of L 2 (Ω), the solution U of the system (15) using the pumping source f given in (4), can be written under the form

U (x 1 , x 2 , t) = n≥0 e n (t)ξ n (x 1 , x 2 )
, where

   e n (t) + µ n e n (t) = ξn(a) √ S(a) (t), ∀t ∈ (0, T ) e n (0) = 0 (16)
Therefore, from [START_REF] Yeh | Stochastic Modelling of Groundwater Flow and Solute Transport in Aquifers[END_REF] it follows that: For all (x 1 , x 2 , t) ∈ Ω × (0, T ),

U (x 1 , x 2 , t) = 1 S(a) n≥0 ξ n (a)ξ n (x 1 , x 2 ) t 0 (η)e -µn(t-η) dη (17)
Then, we establish the following technical result that leads to express the value of the unknown storativity function S at the measuring wells b i=1,...,I ∈ Ω \ {a} in terms of its value at the pumping well a ∈ Ω:

Lemma 2.4 Let a ∈ Ω, ∈ L 2 (0, T ) and f (x 1 , x 2 , t) = (t)δ a (x 1 , x 2 ) be the pumping source employed in the system [START_REF] Garca | A null controllability data assimilation methodology applied to a large scale ocean circulation model[END_REF]. For all b i ∈ Ω \ {a}, it holds

S(b i ) = 1 S(a) B i
where:

B i = n≥0 ξ n (a)ξ n (b i ) µ n T 0 (t) 1 -e -µn(T -t) dt T 0 u(b i , t)dt (18)
Moreover, if the forcing function leads the solution U of the system (15) to satisfy U (•, T ) = 0 a.e. in Ω then, the form in [START_REF] Hoeksema | An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling[END_REF] of the coefficients B i is reduced to

B i = T 0 (t)dt T 0 u(b i , t)dt G b i (a) (19)
where G b i is the impulse response solution of the problem [START_REF] Cooper | A generalized graphical method for evaluating formation constants and summarizing well-field history[END_REF].

Proof. From multiplying the first equation of the system (11) by the solution U of the problem [START_REF] Garca | A null controllability data assimilation methodology applied to a large scale ocean circulation model[END_REF], where f is the pumping source defined in (4), and integrating by parts over Ω using Green's formula, we get: For all t ∈ (0, T ),

U (b i , t) = -div(ψ∇U ) + ρU, G b i = G b i (a) S(a) (t) - d dt U, G b i L 2 (Ω) (20)
where , is the product in the distribution sense. Since

U (b i , t) = S(b i )u(b i , t) for all t ∈ (0, T ) and U (•, 0) = 0 in Ω, it follows from integrating (20) over (0, T ) that S(b i ) = G b i (a) S(a) T 0 (t)dt -U (•, T ), G b i L 2 (Ω) T 0 u(b i , t)dt (21)
If the time-dependent forcing function ∈ L 2 (0, T ) employed in [START_REF] Garca | A null controllability data assimilation methodology applied to a large scale ocean circulation model[END_REF] yields U (•, T ) = 0 a.e. in Ω then, from [START_REF] Meier | A numerical study on the relationship between transmissivity and specific capacity in heterogeneous aquifers[END_REF] we find the form of the coefficients B i announced in [START_REF] Lin | A uniqueness theorem for parabolic equations[END_REF]. Otherwise, from using [START_REF] Hamdi | Detection and identification of multiple unknown time-dependent point sources occurring in 1D evolutionj transport equations[END_REF] to compute the final state U (•, T ) and by replacing the impulse response G b i by its value given in [START_REF] Jai | Capteurs et Actionneurs dans l'Analyse des systèmes distribués[END_REF], we obtain

U (•, T ), G b i L 2 (Ω) = 1 S(a) n≥0 ξ n (a) T 0 (t)e -µn(T -t) dt ξ n , k≥0 ξ k (b i ) µ k ξ k L 2 (Ω) = 1 S(a) n≥0 ξ n (a)ξ n (b i ) µ n T 0 (t)e -µn(T -t) dt (22)
The second equality in [START_REF] Nelson | In-place measurement of permeability in heterogeneous media: Experimental and computational considerations[END_REF] holds since {ξ n } forms an orthonormal family of L 2 (Ω). Afterwards, replacing in ( 21) the term U (•, T ), G b i L 2 (Ω) by its value obtained in [START_REF] Nelson | In-place measurement of permeability in heterogeneous media: Experimental and computational considerations[END_REF] and G b i (a) using ( 13) leads to the result announced in [START_REF] Hoeksema | An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling[END_REF].

Remark 2.5 Since u(x 1 , x 2 , t) = U (x 1 , x 2 , t)/ S(x 1 , x 2 ) for all (x 1 , x 2 , t) ∈ Ω × (0, T ),
it follows from ( 17) and according to Lemma 2.4 that: For all b i ∈ Ω \ {a},

u(b i , t) = 1 B i n≥0 ξ n (a)ξ n (b i ) t 0 (η)e -µn(t-η) dη, ∀t ∈ (0, T ) (23)
Therefore, given time records d i (t) in (0, T ) of the state u solution of the system (1)-( 4) at a measuring well b i ∈ Ω \ {a} and since (µ n ) is an increasing sequence that tends to

+∞, for minimising u(b i , t) -d i (t) 2
L 2 (0,T ) with respect to ψ and ρ, one could truncate the series in ( 23) using a sufficiently large number N of first terms and solve:

min ψ,ρ R N i (ψ, ρ) := 1 2 1 B i N n=0 ξ n (a)ξ n (b i ) t 0 (η)e -µn(t-η) dη -d i (t) 2 L 2 (0,T ) (24)
The solutions ψ and ρ obtained from minimising the sum over all measuring wells b i=1,...,I of the residuals R N i in [START_REF] Richter | An inverse problem for the steady state diffusion equation[END_REF] lead to determine the coefficients B i in ( 18)-( 19) that yield, according to Lemma 2.4, S(a)S(b i ) = B i 2 for i = 1, . . . , I. However, as far as the form of B i used in [START_REF] Richter | An inverse problem for the steady state diffusion equation[END_REF] is concerned, the two main advantages of using ( 19) consist of: 1. Avoiding the approximation of U (•, T ) and G b i (a) done by truncating the series defining 19) can be computed from solving numerically only one time the system [START_REF] Cooper | A generalized graphical method for evaluating formation constants and summarizing well-field history[END_REF] with a instead of b i to determine G a (x 1 , x 2 ), for all (x 1 , x 2 ) ∈ Ω.

B i in (18) 2. Since from (12) it holds G b i (a) = G a (b i ) for all b i ∈ Ω \ {a}, the coefficients G b i (a) in (
Nevertheless, using [START_REF] Lin | A uniqueness theorem for parabolic equations[END_REF] to compute the coefficients B i requires the solution U of the system (15) to fulfill U (•, T ) = 0 a.e. in Ω. To this end, let ϕ m (t) = 2/T sin(mπt/T ) for all t ∈ (0, T ) and m ∈ IN * . It is well known that the set {ϕ m } forms a complete orthonormal family of L 2 (0, T ). Then, we establish the following technical result:

Proposition 2.6 Let a ∈ Ω, M ∈ IN * and (t) = M m=1 m ϕ m (t), ∀t ∈ (0, T ).
From using in the system (15) the pumping source f (x 1 , x 2 , t) = (t)δ a (x 1 , x 2 ), it follows that its final state U N (•, T ) given by truncating the series in [START_REF] Hamdi | Detection and identification of multiple unknown time-dependent point sources occurring in 1D evolutionj transport equations[END_REF] to the order N ∈ IN * is subject to:

U N (•, T ) L 2 (Ω) = 1 S(a) AX 2 (25)
where X = 1 , . . . , M ) ∈ IR M and A is the (N + 1) × M matrix defined by

A nm = ξ n (a) √ 2T mπ e -µnT -(-1) m 1 + T µ n mπ 2 , for n = 0, . . . , N ; m = 1, . . . , M (26)
In ( 26), ξ n and µ n are the normalized eigenfunctions and eigenvalues solutions of [START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF].

Proof. For all n ∈ {0, . . . , N } and m ∈ {1, . . . , M }, let

I nm = T 0 ϕ m (t)e -µn(T -t) dt.
Then, using twice an integration by parts, we get

I nm = √ 2T mπ e -µnT -(-1) m 1 + T µ n mπ 2 (27)
Besides, by employing in the system (15) the pumping source

f (x 1 , x 2 , t) = (t)δ a (x 1 , x 2 )
where (t) = M m=1 m ϕ m (t), ∀t ∈ (0, T ), it follows from truncating the series in [START_REF] Hamdi | Detection and identification of multiple unknown time-dependent point sources occurring in 1D evolutionj transport equations[END_REF] to the order N ∈ IN * that the final state U N (•, T ) of ( 15) satisfies: where A nm = ξ n (a)I nm are the entries of the (N + 1) × M matrix A involved in [START_REF] Rogelio | The differential system method for the identification of transmissivity and storativity[END_REF]. Moreover, from replacing I nm by its value obtained in [START_REF] Schwartz | Théorie des distributions[END_REF], we find the entries A nm announced in [START_REF] Shapoori | Estimating aquifer properties using groundwater hydrograph modelling[END_REF]. Notice that the third equality in ( 28) is obtained since the set {ξ n } forms an orthonormal family of L 2 (Ω).

U N (•, T ) 2 L 2 (Ω) = Ω N n=0 ξ n (x 1 , x 2 ) ξ n (a) S (a)
Therefore, according to Proposition 2.6, in order to drive the solution U of the system (15) towards U (•, T ) = 0 a.e. in Ω, we employ the pumping source f (x 1 , x 2 , t) = (t)δ a (x 1 , x 2 ) defined by (t) = M m=1 m ϕ m (t), ∀t ∈ (0, T ), where the coefficients m are the components of the vector X = 1 , . . . , M solution of the following minimisation problem:

min X∈IR M 1 2 AX 2 2
subject to:

M m=1 m = 1 (29)
In [START_REF] Sterrett | Groundwater and wells, 3rd Eddition[END_REF], A = A nm is the (N + 1) × M matrix defined in [START_REF] Shapoori | Estimating aquifer properties using groundwater hydrograph modelling[END_REF].

Identifiability

We study the identifiability of the two unknown auxiliary variables ψ, ρ in (8) and of the unknown storativity function S using the observation operator M [S, P ] introduced in (6). We start by proving that M [S, P ] yields uniqueness of ρ when this latest is a real number, and of the unknown values of the storativity function at the measuring wells times its value at the pumping well. Then, under two additional assumptions, we establish a second result that yields identifiability for a wider class of unknown functions S and P .

Theorem 3.1 Let ∈ L 2 (0, T ) be such that (t) = 0 a.e. in (0, T ), a ∈ Ω be a pumping well and b i=1,...,I ∈ Ω \ {a} be I ∈ IN * measuring wells. For all unknown functions S and P occurring in the problem ( 1)-( 4) that are elements of the admissible set (3) and for which ρ in ( 8) is a real number, the observation operator M [S, P ] introduced in ( 6) determines uniquely the unknown values of ρ and of S(a)S(b i ), for i = 1, . . . , I.

Proof. Let S (k=1,2) , P (k=1,2) be elements of the admissible set (3) and u (k) be the solution of the system (1)-( 4) with S (k) and P (k) instead of S and P . Besides, let ψ (k) be the function and ρ (k) be the real number defined from S (k) and P (k) as in [START_REF] Castagna | A Bayesianapproach for inversion of hydraulic tomographic data[END_REF]. We denote by (µ

(k)
n ) the eigenvalues and by {ξ

(k)
n } the normalized eigenfunctions of the problem (10) with ψ (k) and ρ (k) instead of ψ and ρ. Since (µ (k) n ) is an increasing sequence of real numbers that tends to +∞, it follows that the series in [START_REF] Hamdi | Detection and identification of multiple unknown time-dependent point sources occurring in 1D evolutionj transport equations[END_REF] defining

U (k) = √ S (k) u (k) converges uniformly in ]τ, +∞[,
for all τ > 0. Therefore, u (k) can be written under the form

u (k) (x 1 , x 2 , t) = t 0 (η)Φ (k) (x 1 , x 2 , t -η)dη, ∀(x 1 , x 2 , t) ∈ Ω × (0, T ) (30)
where the kernel Φ (k) is defined in Ω × (0, T ) by

Φ (k) (x 1 , x 2 , t) = n≥0 ξ(k) n (a) ξ(k) n (x 1 , x 2 )e -µ (k) n t with: ξ(k) n (x 1 , x 2 ) = ξ (k) n (x 1 , x 2 ) S (k) (x 1 , x 2 ) ( 31
)
Let M [S (k) , P (k) ] be the observation operator defined as in [START_REF] Berestycki | Elliptic eigenvalue problems with large drift and applications to nonlinear propagation phenomena[END_REF] from recording in (0, T ) the state u (k) at the measuring wells b i=1,...,I . Thus, we have M [S (2) , P (2) ] = M [S (1) , P (1) ] =⇒ u (2)

(b i , t) = u (1) (b i , t), ∀t ∈ (0, T ), for i = 1, . . . , I (32)
Afterwards, according to ( 30)-( 31), the assertion (32) yields:

For i = 1, . . . , I t 0 (η) Φ (2) (b i , t -η) -Φ (1) (b i , t -η) dη = 0, ∀t ∈ (0, T ) (33)
Since the forcing function is such that (t) = 0 a.e. in (0, T ) and using Titchmarsh's Theorem on convolution of L 1 functions [START_REF] Titchmarsh | Introducton to the theory of Fourier Integrals[END_REF], it follows from (33) that Φ (2) (b i , t)-Φ (1) (b i , t) = 0 a.e. in (0, T ). Hence, in view of [START_REF] Titchmarsh | Introducton to the theory of Fourier Integrals[END_REF], that leads to:

For i = 1, . . . , I n≥0 ξ(2) n (a) ξ(2) n (b i )e -µ (2) n t -ξ(1) n (a) ξ(1) n (b i )e -µ (1) n t = 0, a.e. in (0, T ) (34)
Moreover, since µ (k=1,2) n are both increasing sequences of real numbers that tend to infinity then, the series in [START_REF] Yeh | Time to change the way we collect and analyze data for aquifer characterization[END_REF] converges uniformly in ]τ, +∞[, for all τ > 0. Thus, this series defines a real-valued analytic function of t ∈]0, +∞[. Therefore, in view of [START_REF] Yeh | Time to change the way we collect and analyze data for aquifer characterization[END_REF] and by analytic extension, we get: For i = 1, . . . , I e -µ (2) 0 t ξ( 2)

0 (a) ξ(2) 0 (b i ) - ξ(1) 0 (a) ξ(1) 0 (b i )e -(µ (1) 0 -µ (2) 0 )t + e -µ (2) 0 t n≥1 ξ(2) n (a) ξ(2) n (b i )e -(µ (2) n -µ (2) 0 )t -ξ(1) n (a) ξ(1) n (b i )e -(µ (1)
n -µ

(2)

0 )t = 0, ∀t > 0 (35)
Furthermore, in view of Remark 2.2, it follows that the principal eigenvalue of the problem ( 10) is unique i.e., µ

0 < µ (1)
n and µ

(2)

0 < µ (2)
n , for all n ∈ IN * . Suppose that µ

(1) 0 = µ (2)
0 . Say, for example, it holds µ

(1) 0 > µ (2)
0 which implies that we have also µ

(1) n > µ (2)
0 , ∀n ∈ IN * . Otherwise, we put in (35) rather e -µ (1) 0 t in factor. Afterwards, in (35) from cancelling out e -µ (2) 0 t and setting the limit when t tends to +∞, we obtain ξ( 2)

0 (a) ξ(2) 0 (b i ) = 0. That is absurd since in Ω, S (2)
> 0 and, from Remark 2.2, we have ξ

0 = µ (2)
0 , cancelling out e -µ (2) 0 t and reevaluating the limit when t tends to +∞, we find

   µ (2) 0 = µ (1) 0 ⇔ ρ (2) = ρ (1)
ξ( 2)

0 (a) ξ(2) 0 (b i ) = ξ(1) 0 (a) ξ(1) 0 (b i ) ⇔ S (2) (a)S (2) (b i ) = S (1) (a)S (1) (b i ), for i = 1, . . . , I (36)
The two equivalence results in (36) are obtained from Remark 2.2.

Remark 3.2 Provided the unknown functions S and P occurring in the problem ( 1)-( 4) are elements of the admissible set (3) and such that their ρ in ( 8) is a real number, it follows from Theorem 3.1 that the observation operator M [S, P ] determines uniquely ρ and the coefficients B i in (18) that yield S(a)S(b i ) = B i 2 , for i = 1, . . . , I. Besides, if ψ in ( 8) is also a real number then, provided µ 1 is of multiplicity 1 and ξ 1 (a)ξ 1 (b i 0 ) = 0 where i 0 ∈ {1, . . . , I}, it follows by employing similar techniques as in (35)-(36) that µ

(2)

1 = µ (1)
1 . Since ρ (2) = ρ (1) , that implies ψ (2) = ψ (1) and thus, M [S, P ] yields also uniqueness of ψ.

As far as the identifiability of the unknown auxiliary variables ψ and ρ in (8) for a wider class of unknown functions S and P is concerned, we establish the following result: Theorem 3.3 In the problem (1)-( 4), provided ∈ L 2 (0, T ) satisfying (t) = 0 a.e. in (0, T ) and the unknown functions S and P are elements of (3) generating auxiliary variables ψ and ρ in [START_REF] Castagna | A Bayesianapproach for inversion of hydraulic tomographic data[END_REF] such that the eigenpairs solutions of the system (10) fulfill:

1. The eigenvalues µ n are distinct, for all n ∈ IN . That implies (µ n ) is a strictly increasing sequence of real numbers that tends to +∞. Proof. We use the same notations employed in the proof of Theorem 3.1 and assume that the two sequences (µ

(k=1,2) n
) are both strictly increasing whereas the pumping position a and the measuring position b i 0 are both strategic with respect to {ξ (k=1,2) n }.

• For i = i 0 in (35): Suppose that µ ) are strictly increasing implies that µ

(1) n > µ (2) 0 , ∀n ≥ 1. If µ (1) 0 < µ (2)
0 then, we put in (35) rather e -µ (1) 0 t in factor. From cancelling out the term e -µ (2) 0 t then, setting the limit when t tends to +∞ in (35), we get ξ(2) 0 (a) ξ( 2) 0 (b i 0 ) = 0. That is absurd since a and b i 0 are both strategic with respect to {ξ 0 and cancelling out the term e -µ (2) 0 t then, reevaluating the limit when t tends to +∞, we find ξ( 2)

0 (a) ξ(2) 0 (b i 0 ) = ξ(1) 0 (a) ξ(1)
0 (b i 0 ). Therefore, for i = i 0 the term associated with n = 0 in the series (34) vanishes. By iterating the same process for all n ≥ 1, we obtain

∀n ∈ IN, µ (1)
n = µ (2) n and ξ( 1) n (a) ξ( 1)

n (b i 0 ) = ξ(2) n (a) ξ(2) n (b i 0 ) (37)
• Setting µ

(1) [START_REF] Yeh | Time to change the way we collect and analyze data for aquifer characterization[END_REF]: It follows by analytic extension and putting in factor e -µ 0 t that: For i = 1, . . . , I and all t > 0,

n = µ (2) n = µ n , ∀n ∈ IN in
e -µ 0 t ξ(2) 0 (a) ξ(2) 0 (b i ) - ξ(1) 0 (a) ξ(1) 0 (b i ) + n≥1 ξ(2) n (a) ξ(2) n (b i ) -ξ(1) n (a) ξ(1) n (b i ) e -(µn-µ 0 )t = 0(38)
Afterwards, in (38), from cancelling out the term e -µ 0 t then, setting the limit when t tends to +∞, we find ξ( 2) 1) 0 (b i ). Moreover, iterating the same process for all n ≥ 1 leads to

0 (a) ξ(2) 0 (b i ) = ξ(1) 0 (a) ξ(
∀n ∈ IN, ξ (1)
n (a) ξ( 1)

n (b i ) = ξ(2) n (a) ξ(2) n (b i ), for i = 1, . . . , I (39)
Therefore, from (37) and (39), it follows that

∀n ∈ IN, µ (1

Remark 3.4

The analysis of the results in Theorem 3.3 leads to point out the two following remarks that yield uniqueness of S(a)S(b i ), for i = 1, . . . , I as well as of ψ and ρ for a wide class of unknown functions S and P : k) in Ω and the second equality in (40) holds for all n ∈ IN , we believe that this equality would imply that S (2) (a)S (2) (b i ) = S (1) (a)S (1) (b i ), for i = 1, . . . , I in a large class of auxiliary variables ψ and ρ. Notice that according to Theorem 3.1, the implication holds true when ρ is a real number, for all ψ.

1. Since ξ(k) n = ξ (k) n / √ S (
2. Provided S (2) (a)S (2) (1) (b i ), for i = 1, . . . , I and using a sufficiently large number I of measuring wells, the second equality in (40) implies that there exists a subset N ⊂ IN such that: ξ

(b i ) = S (1) (a)S
n = ξ (1)
n a.e. in Ω, for all n ∈ N . Afterwards, by setting ξ

n = µ (2)
n , it follows from [START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF] that

Ω ψ (2) -ψ (1) ∇ξ n 2 2 + Ω ρ (2) -ρ (1) ξ 2 n = 0, ∀n ∈ N (41)
For example, when ρ (k=1,2) are two real numbers, it follows from Theorem 3.1 that ρ (2) = ρ (1) . In this case, provided N is a non-empty set and ∇ξ n

Identification method

In this section, using the observation operator M [S, P ] introduced in (6), we develop an identification method that leads to determine the two unknown functions S and P occurring in the problem (1)-( 4). Under some conditions on the number and the locations of the employed measuring wells, the developed method starts by determining the two unknown intermediate variables ψ and ρ in (8) from minimising the sum over all measuring wells of the residuals R N i defined by [START_REF] Richter | An inverse problem for the steady state diffusion equation[END_REF]. Then, the determined ψ and ρ lead to compute the coefficients B i in ( 18)-( 19) that, according to Lemma 2.4, yield S(a)S(b i ) = B i 2 for i = 1, . . . , I. Afterwards, we reconstruct the unknown function S and thus, in view of (8), deduce P = Sψ. To this end, we propose the two following ways for the reconstruction of the unknown storativity function S in Ω:

Local determination of S

This first way of determining the unknown function S is based on assuming that it holds div(ψΨ) = 0 in Ω. In fact, from dividing the equations defining the problem ( 1)-( 4) by S, it follows that the state u solves also the system:

           ∂ t u -div ψ∇u -ψΨ∇u = (t) S(a) δ a (x 1 , x 2 ) in Ω × (0, T ) u(•, 0) = 0 in Ω ψ∇u • ν = 0 on ∂Ω × (0, T ) (42)
where the vector field ψΨ stands for the advection term. Therefore, the sense of the assumption div(ψΨ) = 0 in Ω follows from the fact that water is an incompressible fluid. Moreover, in view of [START_REF] Castagna | A Bayesianapproach for inversion of hydraulic tomographic data[END_REF], this assumption reduces the equation satisfied by the unknown function S into the following first order nonlinear partial differential equation:

∂ x 1 S S 2 + ∂ x 2 S S 2 = 4 ρ ψ in Ω (43)
Since the two variables ψ and ρ have been already identified, it follows that the equation (43) could inform us about the local distribution in Ω of the unknown function S. Thus, once the coefficients B i in ( 18)-( 19) are computed from the identified ψ and ρ, we establish a local determination procedure of the function S that combines the knowledge on its distribution obtained from (43) with the knowledge in (18) of S(a)S(b i ) = (B i ) 2 , for i = 1, . . . , I. This local determination proceeds as follows:

Algorithm: Local determination 1.
Find an open subset ω 0 ⊆ Ω that contains a measuring well b i 0 , the pumping well a and within which ρ = 0. It follows from (43) that ∇S = 0 in ω 0 . Since S(a)S(b i 0 ) = (B i 0 ) 2 then, S(x 1 , x 2 ) = B i 0 in ω 0 and thus, S(a) = B i 0 . Furthermore, from dividing the known products S(a)S(b i ) by S(a), we deduce the value of S(b i ), for i = 1, . . . , I.

2.

For all open subset ω ⊂ Ω containing a measuring well b i = (b i 1 , b i 2 ) and within which the identified function ρ/ψ admits symmetric variations: There exists two real numbers α, β and a non-negative derivable real-valued function h :

IR -→ IR + such that ρ x 1 , x 2 ψ x 1 , x 2 = h α(x 1 + x 2 ) + β in ω (44)
which implies that it holds

∂ x 1 ρ(x 1 , x 2 ) ψ(x 1 , x 2 ) = ∂ x 2 ρ(x 1 , x 2 ) ψ(x 1 , x 2 ) in ω (45)
the solution in ω of the equation ( 43) is given by

S(x 1 , x 2 ) = S(b i ) exp x 1 b i 1 2 ρ(η, x 2 ) ψ(η, x 2 ) dη + x 2 b i 2 2 ρ(b i 1 , ζ) ψ(b i 1 , ζ) dζ in ω ( 46
)
where exp stands for the exponential function. According to (44)-( 45), it follows that the solution S obtained in (46) applies, in particular, for all region ω ⊂ Ω where it holds either ρ = 0 or ρ/ψ is equal to a constant. Moreover, in these two particular cases, from (46) it comes that: For all ω ⊂ Ω containing a measuring well b i = (b i 1 , b i 2 ) within which

• ρ = 0 in ω =⇒ S(x 1 , x 2 ) = S(b i ) in ω • ∇ ρ ψ = 0 in ω =⇒ S(x 1 , x 2 ) = S(b i ) exp 2 ρ ψ x 1 -b i 1 + x 2 -b i 2 in ω (47)
However, if (45) doesn't apply then, in view of ( 18) and (43), we search for S from solving:

     min S>0 1 2 ∂ x 1 S S 2 + ∂ x 2 S S 2 -4 ρ ψ 2 L 2 (ω)
Subject to: S(a)S(b i ) = B i 2 , for all b i ∈ ω (48)

Remark 4.1 The existence of the subset ω 0 in 1. could be ensured by setting a measuring well b i 0 as close as possible to the pumping well a in order to get these two wells lying in a small region of Ω where the storativity S remains constant.

Global determination of S

By searching for a solution to the equation ( 9) under the form S(x 1 , x 2 ) = e G(x 1 ,x 2 ) in Ω, it follows that the unknown function G should solve the second order elliptic nonlinear partial differential equation with gradient terms defined by:

∆G + 1 2 ∇G 2 2 + 1 ψ ∇ψ • ∇G = 2 ρ ψ in Ω (49)
In the litterature, solving (49) appears to be a challenging task since the existence and the behaviour near the boundary ∂Ω of its solutions rely on the growth of (1/ψ)∇ψ in Ω, the regularity of ∂Ω and of 2ρ/ψ, see [START_REF] Alarcon | Keller-Osserman type conditions for some elliptic problems with gradient terms[END_REF][START_REF] Radulescu | Singular phenomena in nonlinear elliptic problems: from boundary blow-up solutions to equations with singular nonlinearities, Handbook of Differential Equations: Stationary Partial Differential Equations[END_REF]. Therefore, in view of ( 18)-( 19), we determine rather an approximation g of the unknown function G based on the knowledge of e G(a)+G(b i ) = (B i ) 2 , for i = 1, . . . , I. Thus, we consider that S(x 1 , x 2 ) ≈ e g(x 1 ,x 2 ) in Ω, where the unknown function g is a real polynomial subject to:

g(a) + g(b i ) = 2 ln(B i ), for i = 1, . . . , I (50)
Hence, (50) yields a system of I linear equations on the N g ∈ IN * unknown coefficients defining the sought polynomial g. Provided I ≥ N g and the measuring positions b i=1,...,I are set in Ω \ {a} such that N g equations of (50) are linearly independent, the unknown polynomial g is uniquely determined from (50). For example, in the case when I ≥ N g = 3 i.e., g is a real polynomial of degree 1: g(x 1 , x 2 ) = g 1 x 1 + g 2 x 2 + g 0 , it follows from (50) that the unknown coefficients g 0 , g 1 and g 2 defining g are subject to:

     a 1 + b 1 1 a 2 + b 1 2 2 a 1 + b 2 1 a 2 + b 2 2 2 a 1 + b 3 1 a 2 + b 3 2 2           g 1 g 2 g 0      = 2      ln(B 1 ) ln(B 2 ) ln(B 3 )      (51)
The determinant of the 3 × 3 matrix in (51) is:

b 3 2 b 2 1 -b 1 1 + b 3 1 b 1 2 -b 2 2 + b 1 1 b 2 2 -b 1 2 b 2 1 .
Selecting the positions of the measuring wells b 1 , b 2 and b 3 in a way that affects this determinant to be non-null leads to uniquely determine the unknown coefficients g 0 , g 1 and g 2 from solving the linear system (51). Thus, the global determination proceeds as follows:

Algorithm: Global determination Begin 1. Let g be a desired polynomial defined by N g ∈ IN * unknown coefficients.

Set the measuring wells in the domain Ω \ {a} such that:

i) Their number I ≥ N g .

ii) Their positions yield: N g equations of (50) are linearly independent.

2. Use the identified ψ and ρ to compute B i in ( 18)-( 19), for i = 1, . . . , I.

3. Determine g from solving N g linearly independent equations of (50).

4. Set S(x 1 , x 2 ) ≈ e g(x 1 ,x 2 ) in Ω. End Therefore, provided the number I and the positions of the measuring wells b i=1,...,I fulfill i) and ii), the global determination gives an approximation of the unknown storativity S i.e., S(x 1 , x 2 ) ≈ e g(x 1 ,x 2 ) , where g is a real polynomial of degree up to the user.

Procedure for the identification of S and P

For the clarty of our presentation, in this subsection we summarize the main steps defining the developed identification method. Let {ζ j } be a complete orthonormal family of L 2 (Ω) and (J, M, N ) ∈ (IN * ) 3 be sufficiently large numbers of first terms. The identification of the two unknown functions S and P proceeds in the following four steps:

• Step 1. Select the pumping source and the measuring wells by fulfilling:

1. Hydraulic forcing function ∈ L 2 (0, T ) such that Theorem 3.1 applies.

2. Pumping position a ∈ Ω and a measuring position b i 0 ∈ Ω \ {a} such that Remark 4.1 applies. This is required for using Local Determination of S.

3. The number I of measuring wells b i and their positions in Ω\{a} such that Algorithm Global Determination applies.

Employ the pumping source f (x 1 , x 2 , t) = (t)δ a (x 1 , x 2 ) in Ω × (0, T ). Record in (0, T ) the resulting drawdown d i (t) = u(b i , t) at the measuring well b i , for i = 1, . . . , I.

• Step 2. In the eigenvalue problem [START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF], set ψ(

x 1 , x 2 ) = J j=0 ψ j ζ j (x 1 , x 2 )
. Determine the coefficients ψ 0 , . . . , ψ J and ρ that solve the following minimisation problem:

min ψ 0 ,...,ψ J ,ρ I i=1 R N i (ψ, ρ) subject to: J j=0 ψ j ζ j (x 1 , x 2 ) > 0 in Ω (52)
where R N i is the residual, associated to the measuring well b i , introduced in (24). • Step 3. Identification of the two unknown functions S and P : Computation of the coefficients B i that yield S(a)S(b i ) = B i 2 , for i = 1, . . . , I: Option 1. Use ψ and ρ determined in step 2 to compute the coefficients B i for i = 1, . . . , I from (18).

Option 2. Using ψ and ρ identified in step 2 do:

1. Apply Proposition 2.6 to select ∈ L 2 (0, T ) such that (t) = 0 a.e. in (0,T) and yields U (•, T ) = 0 a.e. in Ω. 2. Reforce the system with and record d i (t) = u(b i , t), ∀t ∈ (0, T ) for i = 1, . . . , I. 3. Compute the coefficients B i=1,...,I from [START_REF] Lin | A uniqueness theorem for parabolic equations[END_REF].

Determine the unknown function S in Ω using the local or the global determination.

• Step 4. Deduce the unknown function P (x 1 , x 2 ) = S(x 1 , x 2 ) J j=0 ψ j ζ j (x 1 , x 2 ) in Ω.
Remark 4.2 Regarding the developed identification method, we point out the following:

1. In step 3., computing the coefficients B i=1,...,I from (19) enables to avoid the approximation done by truncating the series in [START_REF] Hoeksema | An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling[END_REF]. Moreover, according to (12), the coefficients G b i (a) in ( 19) can be computed for i = 1, . . . , I from solving only once the problem (11) with a instead of b i to compute G a (x 1 , x 2 ), ∀(x 1 , x 2 ) ∈ Ω. Then, the symmetric property [START_REF] Badia | On an inverse source problem for the heat equation: Application to a pollution detection problem[END_REF] implies that G b i (a) = G a (b i ), for i = 1, . . . , I.

2.

To apply the developed identification method for determining ρ as an unknown function, one could consider in the eigenvalue problem [START_REF] Cleveland | Type-Curve Matching Using a Computer Spreadsheet[END_REF]

that ρ(x 1 , x 2 ) = J j=0 ρ j ζ j (x 1 , x 2 )
and solves (52) with respect to the unknown coefficients ψ 0 , . . . , ψ J and ρ 0 , . . . , ρ J .

Numerical experiments

We apply the developed identification method to the case of a rectangular aquifer represented by the domain Ω := (0, L 1 ) × (0, L 2 ), where 0 < L 2 ≤ L 1 . This aquifer is characterized by unknown hydraulic storativity S and transmissivity P functions whose the intermediate variables ψ and ρ defined in [START_REF] Castagna | A Bayesianapproach for inversion of hydraulic tomographic data[END_REF] are two real unknown numbers.

For numerical purposes, due to the different ranges of the two unknown functions S and P which would lead to a significant range difference between the two optimisation variables ψ and ρ in (52), we derive the non-dimensional version of the results established in this paper. To this end, for all (x 1 , x 2 , t) ∈ Ω × (0, T ), let:

x = x 1 L 1 , y = x 2 L 2 and s = t T (53)
That reduces the domain of study from Ω × (0, T ) into the unit cube 0, 1 3 . Moreover,

let ū(x, y, s) = u(xL 1 , yL 2 , sT ) = u(x 1 , x 2 , t), S(x, y) = S(xL 1 , yL 2 ) = S(x 1 ,
x 2 ) and P (x, y) = P (xL 1 , yL 2 ) = P (x 1 , x 2 ), for all (x, y, s) ∈ (0, 1) 3 . Then, u solves the problem (1)-( 4) is equivalent to ū satisfies the following system:

         S∂ s ū -div D P ∇ū = T ¯ (s)δ ā(x, y) in 0, 1 3 ū(x, y, 0) = 0 in 0, 1 2 P ∇ū • ν = 0 on ∂ 0, 1 2 × 0, 1 (54)
where ¯ (s) = (sT ) = (t), δ ā is the dirac mass at ā = ( a 1 L 1 , a 2 L 2 ) ∈ 0, 1 2 and D is the diagonal 2 × 2 matrix defined by

D =    T L 2 1 0 0 T L 2 2    (55)
Afterwards, it follows that Ū (x, y, s) = S(x, y)ū(x, y, s) solves: Then, we introduce the following associated eigenvalue problem:

           ∂ s Ū -div D ψ∇ Ū + ρ Ū = T ¯ (s) S(ā) δ ā(x, y) in 0, 1 3 Ū (x, y, 0) = 0 in 0, 1 2 ψ∇ Ū • ν = 0 on ∂ 0, 1 2 × 0, 1 (56)
   -div D ψ∇ξ n + ρξ n = µ n ξ n in 0, 1 2 ψ∇ξ n • ν = 0 on ∂ 0, 1 2 (58)
Since in the case of our numerical experiments ψ and ρ are two unknown real numbers, it follows that the eigenfunctions {ξ nm } and eigenvalues (µ nm ) solutions of (58) are ξ nm (x, y) = c nm cos nπx cos mπy and

µ nm = T ψ nπ L 1 2 + mπ L 2 2 + ρ (59)
for all (x, y) ∈ 0, 1 2 , where (n, m) ∈ IN 2 and c nm are normalizing coefficients:

c nm =          1 if n = m = 0 √ 2 if nm = 0 and n + m > 0 2 if nm = 0 (60) Remark 5.2 From (59), µ n 1 m 1 = µ n 2 m 2 implies that L 2 2 /L 2 1 = (m 2 2 -m 2 1 )/(n 2 1 -n 2
2 ) ∈ Q. Hence, by contraposition it follows that selecting L 1 and L 2 such that L 2 2 /L 2 1 ∈ IR \ Q implies that all eigenvalues µ nm in (59) are of multiplicity equal to 1.

Besides, to generate synthetic measurements, we use: For all (x, y) ∈ 0, 1 2 , S(x, y) = e -γ 1 L 1 x-γ 2 L 2 y-γ 0 and P (x, y) = ψ S(x, y)

In (61), the coefficients γ 0 , γ 1 , γ 2 and 0 < ψ are four real unknown numbers. Therefore, Ψ = -γ 1 L 1 , γ 2 L 2 which, according to (57), implies that

ρ = T ψ 4 γ 2 1 + γ 2 2 =⇒ ρ = T ρ (62)
The implication in (62) is obtained since ψ = ψ and Ψ = -γ 1 , γ 2 . Given the pumping position ā = (ā 1 , ā2 ) ∈ 0, 1 2 , we use the following definition of the dirac mass:

δ ā(x, y) = L 1 L 2 π lim η-→0 + 1 η e -L 2 1 (x-ā 1 ) 2 +L 2 2 (y-ā 2 ) 2 η (63) Afterwards, we determine ŪF (x, y, s) = L 1 L 2 4π ψT s H(s)e - L 2 1 (x-ā 1 ) 2 +L 2 2 (y-ā 2 ) 2 4 ψT s -ρs that solves ∂ s ŪF -ψdiv D∇ ŪF + ρ ŪF = δ 0 (s)δ ā(x, y) in IR 2 × IR (64)
where δ 0 (s) is the dirac mass at s = 0 and H is the Heaviside function [START_REF] Schwartz | Théorie des distributions[END_REF]. Hence, the solution Ū of the system (56), where ψ > 0 and ρ are two real numbers, is given by Ū (x, y, s) = H(s)T ¯ (s) √

S(ā)

s ŪF (x, y, s) + Ū0 (x, y, s), ∀(x, y, s) ∈ 0, 1

3 = L 1 L 2 4π ψ S(ā) s 0 ¯ (s -η) 1 η e - L 2 1 (x-ā 1 ) 2 +L 2 2 (y-ā 2 ) 2 4 ψT η -ρη dη + Ū0 (x, y, s) (65)
where s represents the convolution product with respect to the variable s and Ū0 solves

         ∂ s Ū0 -ψdiv D∇ Ū0 + ρ Ū0 = 0 in 0, 1 3 Ū0 (x, y, 0) = 0 in 0, 1 2 ∇ Ū0 • ν = -H(s)T ¯ (s) √ S(ā) s ∇ ŪF • ν on ∂ 0, 1 2 (66)
We employ the source forcing function ¯ (s) = 0 sin kπs , ∀s ∈ 0, 1 , where the coefficients 0 ∈ IR * and k ∈ IN * . Furthermore, to compute the non-dimensional version of the residuals R N i in [START_REF] Richter | An inverse problem for the steady state diffusion equation[END_REF] and their partial derivatives with respect to the optimisation variables ψ and ρ, we verify that: For all s ∈ 0, 1 , s 0 (η)e -µnm(s-η) dη = 0 µ 2 nm + kπ 2 µ nm sin kπs + kπ e -µnms -cos kπs (67)

Then, for all X ∈ { ψ, ρ}, it follows from (67) that

∂ X s 0 (η)e -µnm(s-η) dη = 0 ∂ X µ nm µ 2 nm + kπ 2 2 Q k n (s), ∀s ∈ 0, 1 (68)
where

Q k n (s) = kπ 2 -µ 2 nm sin kπs + kπ 2µ nm cos kπs -2µ nm + µ 2 nm + kπ 2 s e -µnms (69)
Moreover, we have

1 0 (η) 1 -e -µnm(1-η) dη = 0 1 kπ 1 -(-1) k - kπ µ 2 nm + kπ 2 e -µnm -(-1) k (70)
Hence, in view of [START_REF] Radulescu | Singular phenomena in nonlinear elliptic problems: from boundary blow-up solutions to equations with singular nonlinearities, Handbook of Differential Equations: Stationary Partial Differential Equations[END_REF] and provided U (x, T ) doesn't vanish a.e. in Ω, we get

∂ X 1 µ nm 1 0 (η) 1 -e -µnm(1-η) dη = - ∂ X µ nm µ 2 nm    1 0 (η) 1 -e -µnm(1-η) dη + 0 µ nm µ 2 nm + kπ 2 2 Q k n (1)    (71)
We solved numerically the problem (66) using the five-point finite difference Crank-Nicolson scheme and generated state time records di (s), for all s ∈ 0, 1 and i = 1, . . . , I from (65) such that di (s) = Ū ( bi , s)/ S( bi ), where bi = (

b i 1 L 1 , b i 2 L 2 )
. Then, we solved the non-dimensional version of the minimisation problem (52) using the BFGS quasi-Newton method combined with Wolfe line search. In the sequel, we present the numerical results obtained from solving the non-dimensional version of (52).

• Part1: Identification of the auxiliary variables ψ and ρ

We carried out numerical experiments using a pumping well located at ā = (0.5, 0.5) in the non-dimensional domain 0, 1 2 and forcing the domain with ¯ (s) = 0 sin(kπs) for all s ∈ 0, 1 , where 0 = 1 and k = 4. We used the first N = M = 5 eigenpairs of (58)-( 60) and a total number of N t = 60 measurements taken regularly i.e., with the uniform time step ∆t = T /N t during the monitoring time T , at each of the measuring wells:

Measuring wells b1 b2 b3 b4 b5 b6

Position in 0, 1 2 (0.4, 0.6) (0.6, 0.4) (0.5, 0.2) (0.6, 0.6) (0.4, 0.4) (0.5, 0.8)

Table 1: Positions of the measuring wells in the non-dimensional domain 0, 1 2 .

We initialized the two optimisation variables to ψj = 1 and ρj = 10 -6 . Then, we solved the non-dimensional version of the minimisation problem (52) using measurements, taken at the I first measuring wells of Table 1, generated from the storativity function S introduced in (61) and different values of ψ. Given S, ψ and T , the variable ρ is determined from (62). The obtained numerical results are presented in the following table:

L 2 = 100m, L 1 = √ πL 2 , γ 1 = γ 2 = 10 -2 , γ 0 = 5, T = 1800s.
The analysis of the numerical results in Table 2 shows that the developed identification method leads to identify the two auxiliary variables ψ and ρ with accuracy. This latest seems to be improved by adding more measuring wells in the case where the sought ψ is far away from the initial iterate ψj = 1. Besides, we carried out numerical experiments by considering a constant storativity function i.e., S(x, y) = e -γ 0 , for all (x, y) ∈ 0, 1 2 which implies that ρ = 0. The obtained results are regrouped in the following table:

: L 2 = 100m, L 1 = √ πL 2 , γ 1 = γ 2 = 0, γ 0 = 5 and T = 2400s.
In the case of a constant storativity function, we were able to identify ψ and ρ = 0 only by increasing the final monitoring time T . Indeed, starting from T = 2400s the developed identification method determines with accuracy the variable ψ whereas ρ = 0 is obtained with an opposite sign and relatively small values. In addition, we carried out numerical experiments in the case of a wider domain Ω. Thus, we considered the domain Ω = 0, L 1 × 0, L 2 , where L 2 = 100m and L 1 = π 2 L 2 . We generated measurements using ψ = 137 and the storativity function S in (61) with γ 1 = γ 2 = 10 -2 and γ 0 = 5.

For these experiments, we employed two different final monitoring times and studied the behaviour of the relative errors on the identified ψj i.e., | ψψj |/ ψ and on ρj i.e., |ρ-ρj |/ρ with respect to the used total number N t of measurements taken during the considered monitoring time at each of the measuring wells in Table 1. The results obtained using I = 6 are given by: The behaviour of the relative errors presented in Figure 1 shows that the identified results are improved by increasing the total number N t of measurements taken in each used measuring well during the monitoring time T . The difference between each two curves in Figure 1 observed for relatively small number N t of measurements is due to the numerical method used to approximate the integrals with respect to the time defining the cost function. In our experiments, we used the trapezoidal rule whose the approximation error depends on the time step size ∆t = T /N t .

• Part2: Identification of the storativity function S

As far as the identification of the storativity function is concerned, we apply the algorithm of global determination developed in section 4 to identify the function of (61) that generated the used measurements i.e., S(x, y) = e -γ 1 L 1 x-γ 2 L 2 y-γ 0 in 0, 1 2 . We carried out numerical experiments on identifying S for different values of the coefficients γ 0 , γ 1 and γ 2 . The obtained results are presented in the following two tables: We give for each experiment E the values of γ 0 , γ 1 , γ 2 used to generate the measurements, the associated coefficients B i=1,...,I computed from [START_REF] Hoeksema | An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling[END_REF] and the identified storativity Sident (x, y) = e g(x,y) .

Measurements with

B 1 B 2 B 3 B 4 B 5 B 6
Sident (x, y) = e g 0 E 1 :

γ 1 = 0 γ 2 = 0 γ 0 = 5
10 -3 × 6.67 6.68 6.67 6.68 6.68 6.69 g 0 = -5.01 For each of the two experiments E 1,2 presented in Table 4, the computed coefficients B i=1,...,I have about the same value. Thus, using the first step of the global determination algorithm, let N g = 1 i.e., g(x, y) = g 0 in Ω. It follows that setting g 0 = ln(B 1 ) satisfies with respect to a certain tolerance all equations in (50). Therefore, we set the identified storativity function to Sident (x, y) = e g 0 = B 1 . Moreover, using the results given in Table 4, we determine g 0 for the experiment E 1 from g 0 = ln 6.67 × 10 -3 = -5.01 and for E 2 from g 0 = ln 4.41 × 10 -5 = -10.03.

The numerical results obtained for the identification of non-constant storativity are:

Measurements with B 1 B 2 B 3 B 4 B 5 B 6
Sident (x, y) = e g(x,y) Since the coefficients B i=1,...,I associated to each of the two experiments E 3,4 presented in Table 5 don't have the same value, it follows that a polynomial g with N g = 1 cann't satisfy all the equations in (50). However, as the positions of the three measuring wells b 1 , b 2 and b 3 are such that the 3 × 3 matrix of the linear system (51) is invertible, let N g = 3 i.e., g(x, y) = g 1 L 1 x + g 2 L 2 y + g 0 . The coefficients g 0 , g 1 and g 2 presented in Table 5 have been determined from solving for each experiment the linear system in (51).

Conclusion

We developed an identification approach that leads to reconstruct unknown storativity and transmissivity functions occurring in 2D groundwater equations. Using an appropri-ate change of variables, we transformed the groundwater equation into a diffusion-reaction one, where the diffusion term is the fraction transmissivity/storativity whereas the reaction coefficient yields the right hand side of a second order nonlinear elliptic partial differential equation satisfied by the unknown storativity function. Under some conditions on the used pumping source as well as on the number and the locations of the employed measuring wells, the developed approach starts by identifying the introduced diffusion and reaction variables. Then, it proposes local and global determination procedures for reconstructing the unknown storativity function. The unknown transmissivity is then deduced from the product of the already determined storativity and fraction transmissivity/storativity functions. The numerical results carried out in this paper show that the developed approach determines accurately the used storativity and fraction functions.

1 /

 1 S(Ω). Hence, it follows that µ

2 .

 2 There exists a pumping position a ∈ Ω and a measuring position b i 0 ∈ Ω \ {a} that are both strategic with respect to {ξ n }.The observation operator M [S, P ] introduced in (6) yields uniqueness for all n ∈ IN of µ n and of ξn (a) ξn (b i ) for i = 1, . . . , I, where ξn = ξ n / √ S in Ω.

  ) n = µ(2) n and ξ[START_REF] Abdel-Aziz | Parameter estimation of pumping test data using genetic algorithm[END_REF] n (a) ξ(1)n (b i ) = ξ(2) n(a) ξ(2) n (b i ), for i = 1, . . . , I (40) The results in (40) are those announced in Theorem 3.3.

Remark 5 . 1

 51 Since ψ = ψ, it follows from (57) that the non-dimensional version ρ of the auxiliary variable ρ in (8) could indicate how do the final monitoring time T should be selected to keep the two optimisation variables ψ and ρ having about the same range. That leads to enhance the minimisation of the non-dimensional version of R N i in (52).

Figure 1 :

 1 Figure 1: (a) Relative errors on ψ (b) Relative errors on ρ

Table 4 :

 4 Constant storativity: L 2 = 100m, L 1 = √ πL 2 , T = 2400s, I = 6 and ψ = 71.

5 10 - 4 × 2 g 2 = 2 g 0 = - 4 5 10 - 4 × 4 . 3 g 2 =Table 5 :

 5422204544325 10.17 8.32 10.72 7.53 11.24 7.94 g 1 = -3.02 × 10 --1.01 × 10 -51 7.44 12.66 4.51 7.44 2.82 g 1 = -2.10 × 10 --5.21 × 10 -2 g 0 = -4.74 Varying storativity: L 2 = 100m, L 1 = √ πL 2 , T = 1800s, I = 6 and ψ = 71.

Table 2 :

 2 Measurements:

Table 3 :

 3 Measurements

= 0 a.e. in Ω, the result (41) yields uniqueness of ψ in the class of functions whose the difference is a function that keeps a constant sign a.e. in Ω.


Ba Chien
email: tbchien@dut.udn.vn

Anissa Mokraoui
email: anissa.mokraoui@univ-paris13.fr

Tone mapped HDR images contrast enhancement using piecewise linear perceptual transformation

Keywords: High Dynamic Range image, Low Dynamic Range image, Tone mapping, Histogram, Contrast enhancement

This paper addresses the conversion problem of High Dynamic Range (HDR) images into Low Dynamic Range (LDR) images. The aim of this conversion is to ensure a good visual rendering of the displayed Tone Mapped (TM) HDR images in accordance with the observers' assessment. To do so, the proposed algorithm adjusts a piecewise linear function to the logarithm luminance distribution of the HDR image to adapt this transformation to the perceptual quantizer in accordance with the Human Visual System (HVS). The computation of the slope value of the piecewise linear transformation is simplified. Indeed, it is given by the ratio between the probability and the logarithm luminance distance in the sampling bin. This leads to histograms which tend to be flat, thus inducing a contrast enhancement of the tone mapped HDR images in both under-exposed and overexposed areas. Simulation results provide good results, both in terms of visual quality and TMQI metric, compared to existing competitive TM approaches.

I. INTRODUCTION

Acquiring visual information (e.g. images, videos) of real scenes representing the full range of luminance values of real world (i.e. including both under-and over-exposed areas) is possible through High Dynamic Range (HDR) capture devices. However, the cost of these devices remain too expensive. Therefore to continue viewing HDR content on conventional Low Dynamic Range (LDR) display devices, requires the reduction of the Dynamic Range (DR) of the HDR content using specific tone mapping algorithms. Indeed the main purpose of HDR image Tone Mapping (TM) approach is to find a trade-off between the relevant information (e.g. details, contrast, brightness...) to be preserved or discarded in the image so as to ensure a good visual quality of the displayed image on LDR devices that would be appreciated by observers.

A state of the art on HDR image TM approaches is fairly complete in [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF], [START_REF] Dufaux | High Dynamic Range Video 1st Edition : From Acquisition, to Display and Applications[END_REF] and [START_REF] Reinhard | High Dynamic Range Imaging 2nd Edition: Acquisition, Display, and Image-Based Lighting[END_REF]. Among the developed TM strategies, this paper quickly reviews those that caught our attention because of their performance. Some approaches aim only at preserving detail or enhancing contrast, while others prefer to improve both. In what follows, the described approaches aimed to highlight mainly the details. A second generation of wavelets based on the edge content of the image avoiding having pixels from both sides of an edge has been proposed in In [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF]. In [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF], a separable non-linear multiresolution approach based on essentially non-oscillatory interpolation strategy has been developed. These approaches take into account the singularities (e.g. edge points) in their mathematical models preserving then the structural information of the HDR images. A non-separable version has been investigated and proposed in [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF]. The results provided in [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF], [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF] and [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF] show that the decomposition of the HDR image on different resolution levels would seem to be a good strategy. However, the choice of the decomposition filters is extremely important and decisive in the extractive power of the approximation and detail information.

However other TM strategies aimed both to preserve detail and enhance contrast. Indeed in [START_REF] Durand | Fast Bilateral Filtering for The Display of High-dynamic-range Images[END_REF], an edge-preserving bilateral filter is introduced to decompose the HDR image into two layers: a base layer encoding large-scale variations and a detail one. Contrast is then reduced only in the first layer while the details are kept unchanged. In [START_REF] Drago | Adaptive Logarithmic Mapping for Displaying High Contrast Scenes[END_REF], an adaptive logarithmic mapping method of luminance values is presented. It concerns the adjustment of the logarithmic basis depending on the radiance of the pixels. In [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF], a subband architecture related on an oversampled Haar pyramid representation is proposed. Subband coefficients are re-scaled according to a gain control function reducing the high frequency magnitudes and boosting low ones.

Below are described methods that enhance contrast using in most cases histograms. In [START_REF] Duan | Tone-mapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF], a TM optimization approach using a histogram adjustment between linear mapping and the equalized histogram mapping is developed. A modification of this approach is made in [START_REF] Husseis | Revisited Histogram Equalization as HDR Images Tone Mapping Operators[END_REF]. The latter considers both histogram equalization and human sensitivity to the light function. A piecewise linear perceptual quantizer modeling the "s-shaped" curve of the Human Visual Adaptation with getting more optional contrast by changing the norm space is proposed in [START_REF] Thai | Piecewise linear perceptual quantizer as a non-uniform histogram equalization adjustment for contrast enhancement of tone mapped HDR images[END_REF]. Later, in [START_REF] Thai | Contrast Enhancement and Details Preservation of Tone Mapped High Dynamic Range Images[END_REF], authors proposed a combination of a near optimal lifting scheme and perceptual quantizer to enhance the contrast and preserve more details in the tone mapped HDR images.

This paper is organized as follows. Section II presents the the proposed TM HDR image approach. The later is based on the adjustment of the contrast according to a perceptual piecewise linear quantizer. Section III discusses the simulation results. Section IV concludes the paper.

II. PROPOSED TONE MAPPING HDR IMAGE APPROACH

The purpose of this section is to find the TM HDR values. For this, the proposed approach adjusts locally the distribution of the HDR image logarithm luminance values according to the HVS to enhance the contrast using a simplified piecewise linear function compared to [START_REF] Thai | Piecewise linear perceptual quantizer as a non-uniform histogram equalization adjustment for contrast enhancement of tone mapped HDR images[END_REF].

The original HDR image is assumed to be of size N × M . Denote l HDR the HDR image luminance. It is considered in the logarithm domain since the latter is well adapted to the HVS. It is denoted I HDR and defined as follows:

I HDR := {I HDR (x n , y m ) = log 10 ( l HDR (x n , y m )) (1)
for

1 ≤ n ≤ N and 1 ≤ m ≤ M }, (2)
where I HDR (x n , y m ) is the HDR logarithm luminance value of the pixel located at position (x n , y m ) on the image.

A. Non-uniform distribution of the HDR logarithm luminance

The I HDR values are first sorted and classified into equal B bins defined by cutting points denoted c i uHDR (with 1 ≤ i ≤ B). The index i refers to the bin number. A nonuniform histogram equalization is also performed. c i nuHDR (with 1 ≤ i ≤ B) cutting points, defining the bounds of the non-uniform consecutive B bins, are deduced. From uniform and non-uniform cutting points, new cutting points locally adapted to the logarithm luminance values are deduced as described below.

The new lower bound (i.e. new cutting point) of each bin, denoted l i HDR [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF], is then adjusted by the introduction of the parameter β i as follows:

l i HDR (1) = c i uHDR + β i (c i nuHDR -c i uHDR ), (3)
where β i is a positive parameter smaller than 1 depending on the sub-interval [c i uHDR , c i nuHDR ]. Since the lower cutting point l i HDR (1) is deduced as the average of the logarithm luminance values on each sub-interval [c i uHDR , c i nuHDR ] as in [START_REF] Husseis | Revisited Histogram Equalization as HDR Images Tone Mapping Operators[END_REF], the parameter β i is then given by:

β i = mean(I HDR [c i uHDR , c i nuHDR ]) -c i uHDR c i nuHDR -c i uHDR . (4)
Note that this strategy avoids to get empty bins. The I HDR values are then classified into non-uniform B bins as follows:

I HDR = {l HDR (k) for k = 1, ..., N × M } = {[l 1 HDR (1), ..., l 1 HDR (K 1 )], ..., [l i HDR (1), ..., l i HDR (K i )], ..., [l B HDR (1), ..., l B HDR (K B )]}, (5)
depending on the quantization level set, where K i is the number of values in the i-th bin (i.e. 1 ≤ i ≤ B; K i > 0) and satisfying the following relation

B i=1 K i = N × M. l 1 (1) =l LDRmin l B (K B ) =l LDRmax l 1 (1) =l HDRmin l B (K B ) =l HDRmax l i (1) l i (k) l i (1)
l i+1 (1)

l i (k)
HDR logarithm luminance values LDR luminance values

d 1 d B l i+1 (1)

Normalized histogram

Fig. 1: TM approach using piecewise linear curve.

B. Piecewise linear perceptual quantizer

The "s-shaped" TM perceptual curve, as discussed in [START_REF] Dowling | The Retina: An Approachable Part of the Brain[END_REF] and [START_REF] Geisler | Effects of Bleaching and Backgrounds on the Flash Response of the Cone System[END_REF], is modelled by a piecewise linear curve on each bin (see Fig. 1). Consider the i-th bin, defined by [l i HDR (1), ..., l i HDR (K i )], the TM values are then modeled as follows:

l i LDR (k) = a i l i HDR (k) + b i with k ∈ [1, K i ], (6)
where a i (with a i = 0) and b i are two unknown parameters depending on the i-th bin. Introduce δ i as the difference between HDR luminance in two consecutive bins (see Fig. 1):

δ i = l i+1 HDR (1) -l i HDR (1). (7)
Denote l LDRmax (respectively l LDRmin ) the maximum (respectively minimum) LDR luminance value. The LDR dynamic range (i.e. l LDRmax -l LDRmin ) must satisfy the following relation:

l LDRmax -l LDRmin = B i=1 a i • δ i . (8)
The unknown slope a i has been computed in a sophisticated way in [START_REF] Thai | Piecewise linear perceptual quantizer as a non-uniform histogram equalization adjustment for contrast enhancement of tone mapped HDR images[END_REF]. In this paper, a simplified solution is proposed. Remember that the slope a i is proportional to (l i LDR (k)l i LDR (1)) which is also unknown. The latter can be linked to the probability p i of the l i HDR (k) value in the i-th bin and is given by p i = Ki B i=1 Ki . The probability associated to the LDR dynamic range (i.e. l LDRmax -l LDRmin ) is assumed to be equal to one. Therefore one can deduce that l i LDR (k)l i LDR (1) = p i × (l LDRmax -l LDRmin ). The slope a i is then computed as follows:

a i = l LDRmax -l LDRmin δ i • p i . ( 9
)
The more probability p i , the more value a i , and vice versa.

The set of slopes satisfy the constraint condition given by ( 8). The unknown parameter b i is computed by the following formula: b i = l i LDR (1) -a i × l i HDR [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF]. The HDR mapped values are deduced according to equation [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF]. The global piecewise linear curve is continuous and strictly monotonic increasing according to the positive slopes (i.e a i > 0, or angles 0 • < atan(a i ) < 90 • ).

III. SIMULATION RESULTS

This section discusses the performance of the proposed TM HDR image approach. The TM image objective quality is measured with the "Tone-Mapped image Quality Index" (TMQI) metric [START_REF] Yeganeh | Objective Quality Assessment of Tonemapped Images[END_REF]. Simulations have been conducted under Matlab environnement using the HDR Toolbox ( [START_REF] Banterle | Advanced High Dynamic Range Imaging: Theory and Practice[END_REF]) with 274 test HDR images. For lack of space, we only present the results obtained with 8 HDR images ("Bottle Small", "Office", "Oxford Church", "Memorial", "Light", "ClarkCenter", "Doll-Doll", "ClockBuilding", "AtriumNight" and "PeaceRocks") having different dynamic range (or contrast ratio) from 8 fstops to 29 f-stops.

The results are compared to: (i) CEDP [START_REF] Thai | Contrast Enhancement and Details Preservation of Tone Mapped High Dynamic Range Images[END_REF] with norm 1 and level 1; (ii) NUHA [START_REF] Thai | Piecewise linear perceptual quantizer as a non-uniform histogram equalization adjustment for contrast enhancement of tone mapped HDR images[END_REF] with norm 1; (iii) NONSEP ENO-CA [START_REF] Thai | Image Tone Mapping Approach Using Essentially Non-Oscillatory Bi-quadratic Interpolations Combined with a Weighting Coefficients Strategy[END_REF]; (iv) SEP ENO-CA [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF] with parameters α 1 = 0.3, α 2 = 0.7; (v) Li TMO [START_REF] Li | Compressing and Companding High Dynamic Range Images with Subband Architectures[END_REF] with Haar multiscale; (vi) Fattal [START_REF] Fattal | Edge-Avoiding Wavelets and their Applications[END_REF] using RBW method with parameters α = 0.8, β = 0.3, γ = 0.8; (vii) Duan [START_REF] Duan | Tone-mapping High Dynamic Range Images by Novel Histogram Adjustment[END_REF] using β = 0.5; (viii) TMOs in HDR Toolbox: Drago [START_REF] Drago | Adaptive Logarithmic Mapping for Displaying High Contrast Scenes[END_REF], Reinhard [START_REF] Reinhard | Dynamic Range Reduction Inspired by Photoreceptor Physiology[END_REF], Ward [START_REF] Ward | A Visibility Matching Tone Reproduction Operator for High Dynamic Range Scenes[END_REF], Durand [START_REF] Durand | Fast Bilateral Filtering for The Display of High-dynamic-range Images[END_REF], Schlick [START_REF] Schlick | Quantization Techniques for Visualization of High Dynamic Range Pictures[END_REF] using default parameters as given in the HDR Toolbox. The parameters are chosen so as to give the best results in terms of TMQI metric in all methods.

Table I provides the TMQI metrics. The proposed TM approach namely "Proposed" is deployed with B = 256, l LDRmax = 255 and l LDRmin = 0. The results show that our TM approach is competitive to those developed in the literature. Fig. 2 provides :(i) the TM "ClarkCenter" image of good visual quality in adequacy with the TMQI metric; (ii) normalized histogram with flatness; and (iii) tone mapping curve closest to the "s-shaped" curve. Fig. 3 compares the visual quality of the TM "AtriumNight" images using "Li" and our approach. Although "Li" approach added some post-processing to improve the visual rendering such as cutting off the brightest and darkest parts as well as adding 15% of a histogram equalized layer to the result, our approach is of better visual quality and higher TMQI. Fig. 4 concerns the TM "DollDoll" images using "Fattal" and our approach. These images confirm the performance of our proposal.

The visual quality of the TM "ClockBuilding" images, given by Fig. 5, using "SEPENO" and our approach show that some details are lost on "SEPENO" TM mapped image compared to our approach. Moreover, our TM image is of better contrast. A similar result is provided by Fig. 6 where the "PeaceRock" HDR image has been mapped using "NUHA" and our method. One can observe that the brightness of our TM image is better.

For TM "church" images in Fig. 7, although TMQI metric is slightly different for "Duan" method (0.986) and our approach (0.987), the visual rendering of our TM image is better.

The performance of our approach is confirmed on more than 274 test HDR images where the details and contrast are better represented than other competitive methods for a wide range of f-stops.

IV. CONCLUSION

This paper proposed a simplified piecewise linear perceptual quantizer for tone mapped HDR images contrast enhancement. The HDR logarithm luminance values distribution are adjusted to fit the perceptual "s-shape" curve. Simulation results confirm the relevance of the proposed approach both in terms of the TMQI objective metric and the visual quality of the displayed image.

Fig. 2 :

 2 Fig. 2: TM "ClarkCenter" image (19.42 f-stops).

Fig. 3 :

 3 Fig. 3: TM "AtriumNight" image (28.68 f-stops): (Left) Proposed (TMQI=0.958), (Right) Li (TMQI=0.888).

Fig. 4 :

 4 Fig. 4: TM "DollDoll" image (13.89 f-stops).

Fig. 5 :

 5 Fig. 5: TM "ClockBuilding" image (14.19 f-stops): (Left) Proposed (TMQI=0.962), (Right) SEPENO (TMQI=0.753).

Fig. 6 :

 6 Fig. 6: TM "PeaceRocks" image (24.14 f-stops).

TABLE I :

 I Tone Mapped Image Quality Index (TMQI)

	TMOs	Anturium Bottle Office Church Memorial Light
	DR f-stops	8.73	16.03	16.29	15.46	18.38	17.46
	Drago [10]	0.874	0.801	0.800	0.814	0.800	0.800
	Reinhard [17]	0.778	0.807	0.826	0.789	0.791	0.794
	Ward [18]	0.806	0.783	0.775	0.817	0.795	0.789
	Durand [9]	0.811	0.892	0.825	0.929	0.814	0.800
	Tumblin [19]	0.715	0.713	0.735	0.675	0.759	0.750
	Schlick [20]	0.770	0.835	0.926	0.970	0.787	0.780
	Duan [12]	0.964	0.916	0.955	0.986	0.935	0.916
	Fattal [4]	0.889	0.928 0.943	0.889	0.927	0.971
	Li [11]	0.964	0.954 0.854	0.877	0.834	0.888
	Anas [13]	0.900	0.837	0.896	0.850	0.823	0.832
	SEPENO [5]	0.896	0.934	0.943	0.895	0.932	0.970
	NONSEP [6]	0.938	0.873	0.935	0.820	0.832	0.932
	NUHA [7]	0.929	0.850	0.927	0.897	0.922	0.916
	CEDP [8]	0.957	0.878	0.943	0.980	0.960	0.971
	Proposed	0.975	0.959 0.960	0.987	0.961	0.985