Adel Hamdi

Antoine Tonnoir

Dispersion-Current adjoint functions for monitoring accidental sources in 3D transport equations

Keywords: Nonlinear inverse source problems, Dispersion-Current adjoint functions, Data assimilation, 3D Advection-Dispersion-Reaction equations, Surface water pollution

The paper deals with the identification of multiple unknown time-dependent point sources occurring in 3D dispersion-advection-reaction equations. Based on developed appropriate adjoint functions, we establish a constructive identifiability result depending on the flow nature that yields guidelines leading to a quasi-direct Detection-Identification method. In practice, assuming to be available within a monitored domain some interfaces subdividing it into suspected sections, the developed method goes throughout the entire domain to detect the presence of all active sources. If an activity is detected within a suspected section, the method identifies the total amount discharged in this section and determines whether it is done by a single or multiple unknown occurring sources. Moreover, it localizes the sought position of a detected source as the unique root of a Dispersion-Current vector function defined from the developed adjoint functions. Application to different types of flow and some numerical experiments on surface water pollution are presented.

Introduction

Inverse source problems are usually tasked with identifying unknown/hidden sources that drived the solution of an associated mathematical model to the measured response. Since the achievement of such task leads to illuminate the unknown causes of certain observed effects, we have seen over the last few decades inverse source problems covering a wide range of applications in science and engineering: In medicine, for example, inverse source problems are used to detect and localize a tumour as well as to identify some biological parameters, like the chemotaxis, that lead to simulate numerically the growth of a tumour in order to predict the adapted treatement to heal the patient [START_REF] Hawkins-Daarud | Numerical simulation of a thermodynamically consistent four-species tumor growth model[END_REF][START_REF] Quiroga | Adjoint method for a tumor invasion Pde-constrained optimization problem in 2D using adaptive finite element method[END_REF][START_REF] Gholami | An inverse problem formulation for parameter estimation of a reactiondiffusion model of low grade gliomas[END_REF][START_REF] Fokas | The unique determination of neuronal currents in the brain via magnetoencephalogrphy[END_REF]. For environmental monitoring, inverse source problems are employed to identify pollution sources in surface water [START_REF] Andrle | Identification of Moving Pointwise Sources in an Advection-Dispersion-Reaction Equation[END_REF][START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF][START_REF] Belgacem | Uniqueness for an ill-posed reaction-dispersion model. Application to organic pollution in stream-waters (One-dimensional model)[END_REF][START_REF] Hamdi | Detection-identification of multiple unknown time-dependent point sources in a 2D transport equation: Application to accidental pollution[END_REF], in groundwater [START_REF] Huang | An inverse problem in estimating the strength of contaminant source for groundwater systems[END_REF][START_REF] Gurarslan | Solving inverse problems of groundwaterpollution-source identification using a differential evolution algorithm[END_REF] and in atmosphere [START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF]. Some important applications include earthquake source localization [START_REF] Larmat | Time reversal location of glacial earthquakes[END_REF][START_REF] Kawakatsu | Time reversal seismic imaging and moment tensor[END_REF][START_REF] Koketsu | Inverse Problems. Inverse Problems in Seismology[END_REF], source identification in electromagnetics [START_REF] Benoit | Source identification in time domain electromagnetics[END_REF] and acoustics [START_REF] Idemen | On an inverse source problem connected with photoacoustic and thermo-acoustic tomographies[END_REF], among many others.

As far as the identifiability of unknown time-dependent point sources is concerned, in the case of one dimensional transport equations, the authors in [START_REF] Badia | Identification of a point source in a linear advection-dispersion-reaction equation: application to a pollution source problem[END_REF] proved that for a single unknown occurring source, its position and its time-dependent intensity are determined in a unique way from time records of the generated state and of its flux taken at two positions framing the source region. This identifiability result has been improved in [START_REF] Hamdi | The recovery of a time-dependent point source in a linear transport equation: application to surface water pollution[END_REF] to hold using only state time records then, extended firstly in [START_REF] Hamdi | Inverse source problem in a one-dimensional evolution linear transport equation with spatially varying coefficients: application to surface water pollution[END_REF] to treat the case of transport equations with spatially varying coefficients and secondly in [START_REF] Hamdi | Identification of a time-varying point source in a system of two coupled linear diffusion-advection-reaction equations: application to surface water pollution[END_REF] to address the case of two coupled transport equations. In [START_REF] Belgacem | Uniqueness for an ill-posed reaction-dispersion model. Application to organic pollution in stream-waters (One-dimensional model)[END_REF], the author proved the identifiability of a moving time-dependent point source. Besides, regarding two/three dimensional transport equations, authors in [START_REF] Andrle | Identification of Moving Pointwise Sources in an Advection-Dispersion-Reaction Equation[END_REF][START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF][START_REF] Isakov | Inverse Source Problems[END_REF] established identifiability results for multiple unknown time-dependent point sources. These results based on the unique continuation theorem say that the unknown elements defining all finite number of occurring time-dependent point sources are uniquely determined from state/flux time records taken on whatever nonempty part of the boundary. In our view, these results yield an ideal theoretical framework which in practice does neither take into account the flow nature i.e., transport/dispersion dominance in selecting the appropriate measuring part of the boundary nor gives visibility on how to proceed for determining the unknown elements defining the occurring sources.

This study is motivated by the following two points: 1. We present in section 8 of this paper a Baby Example showing that for the one dimensional transport equation, it is not possible to uniquely identify multiple unknown point sources using measurements taken only upstream and/or downstream all occurring point sources. To ensure identifiability, all two distinct point sources should be separated by either two state measuring points or by one state and flux measuring point 2. The identification method proposed in [START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF] that consists of minimizing two objective functions (Least squares and Kohn-Vogelius) defined using measurements taken on the whole boundary, encountered serious difficulties starting from the case of two occurring sources: The least squares approach doesn't enable to identify more than one active source whereas the Kohn-Vogelius approach identifies two sources only if they are well separated and the diffusion coefficient is very small in order to avoid rapid intermixing. Otherwise, the authors reported in [START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF] that sources effects become quickly indistinguishable. The originality of the results established in the present paper lies in the developement of direct mathematical techniques to determine the unknown elements defining the occurring sources. Indeed, we develop appropriate adjoint functions that lead to establish a constructive identifiability result depending on the flow nature and transforming the localization of an unknown source position into determining the unique root of a Dispersion-Current vector function defined from the developed adjoint functions. Comparing to the minimization of objective functions, the identification approach established in this paper is expected to improve the accuracy of the identified results as well as to reduce the total identification cost since we do not have to solve two PDEs (State and adjoint state) for each optimization iteration and also do neither face the problem of local minimizers nor of adding regularization terms.

The paper is organized as follows: In section 2, we present the problem statement and introduce some assumptions defining the framework of this study. Section 3 is reserved to develop appropriate adjoint functions and to prove some of their properties for later use. In section 4, we establish a constructive identifibility result based on the developed adjoint functions. Section 5 is devoted to follow the guidelines given by the established identifiability result to develop a Detection-Identification method. Applications, Baby Example on the non-identifiability in the one dimensional case and numerical experiments on the surface water BOD pollution model are presented in the remaining sections.

Mathematical modelling and problem statement

Let T > 0 be a final monitoring time and Ω be a bounded and connected open subset of IR 3 with a sufficiently smooth boundary ∂Ω := Γ in ∪ Γ F S ∪ Γ F A ∪ Γ out , where Γ in is the inflow boundary, Γ out is the outflow boundary, Γ F S regroups the Fluid-Solid boundaries and Γ F A is the Fluid-Air boundary. For example, in surface water pollution, the BOD concentration u within a portion Ω of a river is governed by [START_REF] Bear | Introduction to Modeling of Transport Phenomena in Porous Media[END_REF][START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Cox | A review of currently available in-stream water-quality models and their applicability for simulating dissolved oxygen in lowland rivers[END_REF][START_REF] Linfield | The enhanced stream water quality models QUAL2E and QUAL2E-UNCAS: Documentation and user manual[END_REF]:

L[u](x, t) = F (x, t) in Ω × (0, T ), (1)
where F represents the set of all occurring pollution sources and L is the second-order linear partial differential operator defined as follows:

L[u](x, t) := ∂ t u(x, t) -div D(x)∇u(x, t) + V (x) • ∇u(x, t) + Ru(x, t).

(

) 2
The dot product designates the inner product in IR 3 , D is the hydrodynamic dispersion tensor, V = V 1 , V 2 , V 3 is the flow velocity field and R is a real number that represents the reaction coefficient. The tensor D is a 3 × 3 real matrix symmetric, uniformly elliptic and uniformly bounded in Ω. Moreover, D and V satisfy

div(V ) = 0 in Ω, V = 0 on Γ F S and V • ν = 0 on Γ F A . (3)
ν is the unit outward vector normal to ∂Ω. In (3), the first condition stands for the incompressibility of water whereas the second one is the so-called no-slip condition [START_REF] Lauga | Microfluidics: The No-Slip Boundary Condition[END_REF].

Besides, for appending boundary and initial conditions to (1)-( 2), we use [START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF][START_REF] Hamdi | Detection-identification of multiple unknown time-dependent point sources in a 2D transport equation: Application to accidental pollution[END_REF]:

         u(•, 0) = 0 in Ω, u = 0 on in = Γ in × (0, T ), D∇u • ν = 0 on F S ∪ F A ∪ out , (4)
where F S = Γ F S × (0, T ), F A = Γ F A × (0, T ) and out = Γ out × (0, T ). Notice that due to the linearity of the operator L introduced in (2) and according to the superposition principle, the use of a non-zero initial condition and/or inhomogeneous boundary conditions do not affect the results established in this paper.

In the present study, we are interested in the case of multiple time-dependent point sources occurring in the problem (1)-( 4) i.e., F is defined as follows:

F (x, t) = N n=1 λ n (t)δ(x -S n ) in Ω × (0, T ), (5)
where N ∈ IN * , δ denotes the Dirac mass, S n=1,...,N are N distinct interior locations in Ω that represent the positions of the occurring sources and λ n=1,...,N ∈ L 2 (0, T ) designate their associated time-dependent source intensity functions that satisfy:

Then, using the transposition method introduced by Lions [START_REF] Lions | Pointwise control for distributed systems in control and estimation in distributed parameters systems[END_REF][START_REF] Lions | Some methods in the mathematical analysis of systems and their control[END_REF], it follows that the forward problem (1)-( 5) admits a unique solution u that belongs to:

L 2 0, T ; L 2 (Ω) ∩ C 0 0, T ; H 1 (Ω) . (7)
We assume to be available I ∈ IN * sufficiently smooth interfaces Γ i=1,...,I subdividing the domain Ω into I + 1 subdomains Ω i=1,...,I+1 defined by

∂Ω i = Γ i-1 ∪ Γ i F S ∪ Γ i F A ∪ Γ i , with Γ 0 = Γ in and Γ I+1 = Γ out whereas Γ i
F S and Γ i F A are the parts of the boundaries Γ F S and Γ F A situated between Γ i-1 and Γ i . Moreover, the interface Γ I is set in a way such that the last subdomain Ω I+1 remains a "No source region of Ω". Therefore, we have

Ω = ∪ I+1 i=1 Ω i , where Ω I+1 ⊂ Ω \ S 1 , . . . , S N . (8)
In our study, Ω i=1,...,I represent the suspected sections within the monitored domain Ω wherein sources could occur and the interfaces Γ i = ∂Ω i ∩∂Ω i+1 are available for measuring the state u and its flux D∇u • ν crossing each intersection. Since the source positions S n=1,...,N are interior locations in Ω i=1,...,I ⊂ Ω, it follows that the state u solution of the problem (1)-( 5) is smooth enough on ∂Ω as well as on the interfaces Γ i=1,...,I . That allows us to define the following observation operator:

M [F ] := D∇u • ν on in ; u, D∇u • ν on i=1,...,I ; uDν on F A , (9)
where i = Γ i × (0, T ). The inverse source problem in which we are interested here is: Given the records dd in of D∇u•ν on in , the records d i , dd i of u, D∇u•ν on i=1,...,I and the records d F A of u on F A , detect whether there is or not active sources occurring in each suspected subdomain Ω i ⊂ Ω, for i = 1, . . . , I. If the presence of active sources is detected within Ω i , determine the total amount discharged in Ω i and whether it is done by a single or multiple occurring unknown sources. In the case of a single active source occurring within Ω i , localize its sought position i.e., find S n ∈ Ω i that yields

d i , dd i = u, D∇u • ν on i-1 ∪ i and d F A = u on i F A . ( 10
)
Remark 2.1 The last term in [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF] defining the observation operator M [F ] means that if the dispersion tensor D along the boundary Γ F A is such that Dν either null or of norm small enough that can be neglected then, the underlined inverse source problem can be solved without any need of state records taken on the Fluid-Air part Γ F A of the boundary.

Appropriate adjoint functions

Let i ∈ {1, . . . , I}. From multiplying the equations ( 1)-( 2) by a sufficiently regular test function v and integrating by parts over Ω i × (0, T ) using Green's formula, we obtain

T 0 F, v D(Ω i ),D (Ω i ) = L a [v], u L 2 (Ω i ×(0,T )) + A v i , (11)
where , D(Ω i ),D (Ω i ) designates the product in the distribution sense and L a is the adjoint operator associated to the operator L in (2) i.e.:

• L a [v] := -∂ t v -div D∇v -V • ∇v + Rv, • A v i = Ω i u(•, T )v(•, T ) + ∂Ω i ×(0,T ) u D∇v + vV -vD∇u • ν. ( 12
)
Since the state u is subject to only knowledge of M [F ] then, to keep in the right hand side of [START_REF] Huang | An inverse problem in estimating the strength of contaminant source for groundwater systems[END_REF] only known terms, we aim to develop adjoint functions v i that solve the system

         L a [v i ](x, t) = 0 in Ω i × (0, T ), D∇v i + v i V • ν = 0 on Γ i F S × (0, T ), u D∇v i + v i V • ν known from M [F ] on Γ i F A × (0, T ). (13)
From ( 3), the boundary conditions in [START_REF] Badia | Identification of a point source in a linear advection-dispersion-reaction equation: application to a pollution source problem[END_REF] are equivalent to D∇v i • ν = 0 on Γ i F S × (0, T ) and uD∇v i • ν known on Γ i F A × (0, T ). Hence, by searching for v i under the form

v i (x, t) = e Rt v 0i (x) in Ω i × (0, T ) =⇒ -∂ t v i + Rv i = 0 in Ω i × (0, T ), (14)
it follows according to (13) that the spatial part v 0i in ( 14) should solve

         div(D∇v 0i ) + V • ∇v 0i = 0 in Ω i , D∇v 0i • ν = 0 on Γ i F S , uD∇v 0i • ν known from M [F ] on Γ i F A . (15)
Therefore, the first kind of adjoint functions is what we refer to as free space adjoint function defined from v 0i = 1 by v i (x, t) = e Rt in Ω i × (0, T ). In the following two subsections, we develop two other kinds of spatial adjoint functions v 0i solving (15):

Full adjoint function

Since from (3) we have div

(V ) = 0 in Ω i , it follows that div(v 0i V ) = V • ∇v 0i in Ω i .
Afterwards, the first equation in ( 15) is equivalent to: div(D∇v 0i + v 0i V ) = 0 in Ω i . By looking for v 0i under the form v 0i (x) = e ψ i (x) for all x ∈ Ω i , we get

D∇v 0i + v 0i V = 0 in Ω i ⇔ D∇ψ i + V = 0 in Ω i . (16)
Thus, the scalar potential ψ i is well defined from the second equation in [START_REF] Hamdi | Identification of a time-varying point source in a system of two coupled linear diffusion-advection-reaction equations: application to surface water pollution[END_REF] provided the following compatibility condition on D and V holds true:

-→ rot D -1 V = 0 in Ω i . (17)
For the simplicity of our notations, in the remainder of this paper we denote for k = 1, 2, 3 by (D -1 ) k the k th row of the matrix D -1 . Then, we have

-→ rot D -1 V =      ∂ x 2 (D -1 ) 3 • V -∂ x 3 (D -1 ) 2 • V ∂ x 3 (D -1 ) 1 • V -∂ x 1 (D -1 ) 3 • V ∂ x 1 (D -1 ) 2 • V -∂ x 2 (D -1 ) 1 • V      , (18)
where (D -1 ) k • V designates the inner product of the k th row of the matrix D -1 and the velocity field V . Hence, provided [START_REF] Hamdi | Inverse source problem in a one-dimensional evolution linear transport equation with spatially varying coefficients: application to surface water pollution[END_REF] holds true, the scalar potential ψ i fulfilling the second equation in ( 16) is defined for all x = (x 1 , x 2 , x 3 ) ∈ Ω i by

ψ i (x) = - x 1 a i (D -1 ) 1 • V (η, x 2 , x 3 )dη - x 2 b i (D -1 ) 2 • V (a i , ξ, x 3 )dξ - x 3 c i (D -1 ) 3 • V (a i , b i , ζ)dζ, (19)
where (a i , b i , c i ) ∈ Ω i . In addition, for later use we establish the following property of ψ i : Lemma 3.1 Provided ( 17) holds, the scalar potential ψ i defined in [START_REF] Hasanov | Fourier collocation algorithm for identifying the spacewise-dependent source in the advection-diffusion equation from boundary data measurements[END_REF] satisfies: For all x = (x 1 , x 2 , x 3 ) and y = (y 1 , y 2 , y 3 ) in Ω i , ψ i (x) = ψ i (y) implies that

x 1 y 1 (D -1 ) 1 •V (η, x 2 , x 3 )dη+ x 2 y 2 (D -1 ) 2 •V (y 1 , ξ, x 3 )dξ+ x 3 y 3 (D -1 ) 3 •V (y 1 , y 2 , ζ)dζ = 0.
Proof. See the Appendix.

Therefore, in view of (3) and from [START_REF] Hamdi | Identification of a time-varying point source in a system of two coupled linear diffusion-advection-reaction equations: application to surface water pollution[END_REF], it follows that v 0i (x) = e ψ i (x) solves the two first equations of the system (15) and reduces to null its last boundary condition. Thus, we obtain the so-called full adjoint function defined by

v i (x, t) = e Rt+ψ i (x) in Ω i × (0, T ). ( 20
)

Separated adjoint functions

The third way of developing v 0i that solves the system (15) consists of reducing to zero separately each of the two terms defining the first equation in [START_REF] Hamdi | The recovery of a time-dependent point source in a linear transport equation: application to surface water pollution[END_REF]. To this end, let

V = V 1 , V 2 , V 3 be the velocity field in Ω and V ⊥ = V ⊥ 1 , V ⊥ 2 , V ⊥ 3
be a vector field defined from the same components of V and perpendicular to it i.e.,

V ⊥ k=1,2,3 ∈ 0, ±V 1 , ±V 2 , ±V 3 and V • V ⊥ = 0 in Ω. (21)
Afterwards, we define v 0i such that its gradient is collinear to V ⊥ i.e.,

∇v 0i = w i V ⊥ in Ω i =⇒ V • ∇v 0i = 0 in Ω i , (22)
where w i is a weighting function. Besides, in order to satisfy also div D∇v 0i = 0 in Ω i as well as to ensure the compatibility condition -→ rot ∇v 0i = 0 in Ω i , we require the weighting function w i involved in [START_REF] Gurarslan | Solving inverse problems of groundwaterpollution-source identification using a differential evolution algorithm[END_REF] to solve the following system:

   div w i DV ⊥ = 0 in Ω i , -→ rot w i V ⊥ = 0 in Ω i , ⇔    w i div DV ⊥ + ∇w i • DV ⊥ = 0 in Ω i , w i -→ rot V ⊥ + ∇w i ∧ V ⊥ = 0 in Ω i , (23)
where ∧ designates the vector product.

Remark 3.2 In view of [START_REF] Larmat | Time reversal location of glacial earthquakes[END_REF], if the dispersion tensor D and the velocity field V satisfy div DV ⊥ = 0 and -→ rot V ⊥ = 0 in Ω i which is the case, for example, when D and V are defined in Ω i by mean values, then the weighting function can be taken

w i = 1 in Ω i .
Moreover, from searching for a weighting function w i that solves [START_REF] Larmat | Time reversal location of glacial earthquakes[END_REF] under the form: w i = e g i , it follows that the unknown function g i should fulfill in Ω i :

   ∇g i • DV ⊥ = -div DV ⊥ , ∇g i ∧ V ⊥ = - -→ rot V ⊥ , ⇔                ∇g i • DV ⊥ = -div DV ⊥ "scalar equation",      0 V ⊥ 3 -V ⊥ 2 -V ⊥ 3 0 V ⊥ 1 V ⊥ 2 -V ⊥ 1 0      ∇g i = -      ∂ x 2 V ⊥ 3 -∂ x 3 V ⊥ 2 ∂ x 3 V ⊥ 1 -∂ x 1 V ⊥ 3 ∂ x 1 V ⊥ 2 -∂ x 2 V ⊥ 1      . ( 24
)
Notice that the 3 × 3 matrix involved in ( 24) is antisymmetric and thus, its determinant is equal to zero. Then, we consider the following three particular choices of V ⊥ :

• Let V 1 ⊥ = 0, -V 3 , V 2 .
From replacing in the system (24) V ⊥ by V 1 ⊥ and g i by g 1 i , it follows under the condition:

V 3 ∂ x 1 V 2 = V 2 ∂ x 1 V 3 in Ω i , (25)
that the scalar equation can be taken into account in the 3 × 3 linear system of (24) as follows: If V 2 = 0 or V 3 = 0 in Ω i then, we substitute the null equation in the 3 × 3 linear system of [START_REF] Kawakatsu | Time reversal seismic imaging and moment tensor[END_REF] with the scalar equation. If V 2 V 3 = 0 in Ω i then, since [START_REF] Koketsu | Inverse Problems. Inverse Problems in Seismology[END_REF] implies that the two last equations in the 3 × 3 linear system of ( 24) are equivalent, we replace one of those two equations by the scalar equation. For example, if V 2 = 0 in Ω i , replacing the last equation of the 3 × 3 linear system in [START_REF] Kawakatsu | Time reversal seismic imaging and moment tensor[END_REF] by the scalar equation gives

     0 V 2 V 3 -V 2 0 0 (DV 1 ⊥ ) 1 (DV 1 ⊥ ) 2 (DV 1 ⊥ ) 3      ∇g 1 i =      ∂ x 1 V 1 ∂ x 1 V 2 -div DV 1 ⊥      , (26)
where (DV 1 ⊥ ) k is the k th component of the vector DV 1 ⊥ and the first component in the right hand side vector is obtained from applying div(V ) = 0 in Ω i .

• Let V 2 ⊥ = -V 3 , 0, V 1 . We set in (24) V ⊥ = V 2 ⊥ and g i = g 2 i .
Afterwards, provided the following condition holds:

V 3 ∂ x 2 V 1 = V 1 ∂ x 2 V 3 in Ω i , (27)
the scalar equation can be included in the 3 × 3 linear system of (24) as follows: If V 1 = 0 or V 3 = 0 in Ω i , we substitute the null equation in the 3 × 3 linear system of (24) with the scalar equation. [START_REF] Lions | Pointwise control for distributed systems in control and estimation in distributed parameters systems[END_REF] implies that the first and the last equations in the 3 × 3 linear system of ( 24) are equivalent, we replace one of those two equations by the scalar equation. For example, if V 1 = 0 in Ω i , replacing the last equation of the 3 × 3 linear system in [START_REF] Kawakatsu | Time reversal seismic imaging and moment tensor[END_REF] by the scalar equation leads to

If V 1 V 3 = 0 in Ω i , then since
     0 V 1 0 -V 1 0 -V 3 (DV 2 ⊥ ) 1 (DV 2 ⊥ ) 2 (DV 2 ⊥ ) 3      ∇g 2 i = -      ∂ x 2 V 1 ∂ x 2 V 2 div DV 2 ⊥      . ( 28
)
The second component of the right hand side in ( 28) is obtained from div(V ) = 0 in Ω i .

• Let V 3 ⊥ = -V 2 , V 1 , 0 .
From substituting in (24) V ⊥ with V 3 ⊥ and g i with g 3 i , it follows under the condition:

V 2 ∂ x 3 V 1 = V 1 ∂ x 3 V 2 in Ω i , (29)
that the scalar equation can be taken into account within the 3 × 3 linear system of ( 24) as follows: If V 1 = 0 or V 2 = 0 in Ω i , we substitute the null equation in the 3 × 3 linear system of ( 24) with the scalar equation. [START_REF] Lin | A uniqueness theorem for parabolic equations[END_REF] implies that the two first equations in the 3 × 3 linear system of ( 24) are equivalent, we replace one of those two equations by the scalar equation. For example, if V 1 = 0 in Ω i , replacing the second equation of the 3 × 3 linear system by the scalar equation gives

If V 1 V 2 = 0 in Ω i , then since
     0 0 -V 1 (DV 3 ⊥ ) 1 (DV 3 ⊥ ) 2 (DV 3 ⊥ ) 3 V 1 V 2 0      ∇g 3 i =      ∂ x 3 V 1 -div DV 3 ⊥ ∂ x 3 V 3      . ( 30
)
The third component of the right hand side in ( 30) is obtained from div(V ) = 0 in Ω i .

Remark 3.3 For example, if the velocity field V (x 1 , x 2 , x 3 ) = V 1 (x 1 , x 2 ), V 2 (x 1 , x 2 ), 0
in Ω i , then all of the three conditions ( 25), ( 27) and ( 29) are well fulfilled.

Thus, according to [START_REF] Gurarslan | Solving inverse problems of groundwaterpollution-source identification using a differential evolution algorithm[END_REF], we determine the function v 0i associated to g i from solving

∇v 0i = e g i V ⊥ in Ω i , for = 1, 2, 3. (31)
Moreover, (31) leads to: For all

x = (x 1 , x 2 , x 3 ) ∈ Ω i , • v 1 0i (x) = - x 2 b i e g 1 i V 3 (a i , η, x 3 )dη + x 3 c i e g 1 i V 2 (a i , b i , ξ)dξ, • v 2 0i (x) = - x 1 a i e g 2 i V 3 (η, b i , x 3 )dη + x 3 c i e g 2 i V 1 (a i , b i , ξ)dξ, • v 3 0i (x) = - x 1 a i e g 3 i V 2 (η, x 2 , c i )dη + x 2 b i e g 3 i V 1 (a i , ξ, c i )dξ, (32)
where (a i , b i , c i ) ∈ Ω i . In view of ( 22)-( 23), since V ⊥ are defined from the components of V and as according to (3), V = 0 on Γ F S , it follows that the separated adjoint functions v =1,2,3 0i in (32) solve the two first equations of (15) and, according to [START_REF] Machado | A new One-Shot Pointwise Source Reconstruction Method[END_REF], satisfy

D∇v 0i • ν = e g i V ⊥ • Dν on Γ i F A . (33)
(33) holds since D is a symmetric matrix. Hence, all the developed adjoint functions i.e.,

v i (x, t) = e Rt v 0i (x) in Ω i × (0, T ), ∀v 0i (x) ∈ 1, e ψ i (x) , v 1 0i (x), v 2 0i (x), v 3 0i (x) , (34)
solve the two first equations of the system [START_REF] Badia | Identification of a point source in a linear advection-dispersion-reaction equation: application to a pollution source problem[END_REF]. In addition, the free space adjoint function v i (x, t) = e Rt and the full adjoint function v i (x, t) = e Rt+ψ i (x) reduce to null the last boundary condition of (13) whereas the three separated adjoint functions

v i (x, t) = e Rt v 0i (x) yield u D∇v i + V v i • ν = e g i +Rt V ⊥ • uDν on Γ i F A × (0, T ).

Dispersion-Current vector function

Using the developed adjoint functions in (34) and given two distinct integers m and n from {1, 2, 3}, we introduce what we will refer to as Dispersion-Current vector function:

Ψ m,n i : Ω i -→ IR 3 , x → Ψ m,n i (x) =      e ψ i (x) v m 0i (x) v n 0i (x)      , (35)
where ψ i is the scalar potential obtained in [START_REF] Hasanov | Fourier collocation algorithm for identifying the spacewise-dependent source in the advection-diffusion equation from boundary data measurements[END_REF] and v m 0i , v n 0i are two functions from [START_REF] Nara | A projective method for an inverse source problem of the Poisson equation[END_REF]. Moreover, to prove the injectivity of Ψ m,n i , we establish the following technical result:

Lemma 3.4 Let V = V 1 , V 2 , V 3 be a vector field such that V k=1,2,3 ∈ W 1,∞ (Ω)
. For all two elements x = (x 1 , x 2 , x 3 ) and y = (y 1 , y 2 , y 3 ) of Ω i , we have

• v 1 0i (x) = v 1 0i (y) =⇒ x 2 y 2 e g 1 i V 3 (a i , η, x 3 )dη = x 3 y 3 e g 1 i V 2 (a i , y 2 , ξ)dξ, • v 2 0i (x) = v 2 0i (y) =⇒ x 1 y 1 e g 2 i V 3 (η, b i , x 3 )dη = x 3 y 3 e g 2 i V 1 (y 1 , b i , ξ)dξ, • v 3 0i (x) = v 3 0i (y) =⇒ x 1 y 1 e g 3 i V 2 (η, x 2 , c i )dη = x 2 y 2 e g 3 i V 1 (y 1 , ξ, c i )dξ, (36)
where v =1,2,3 0i are the three adjoint functions obtained in [START_REF] Nara | A projective method for an inverse source problem of the Poisson equation[END_REF].

Proof. See the Appendix.

That leads to prove the following theorem on the injectivity in a suspected subdomain Ω i ⊂ Ω of the Dispersion-Current function Ψ m,n i introduced in (35):

Theorem 3.5 Let V = (V 1 , V 2 , V 3 )
and D be such that Lemmas 3.1 and 3.4 apply. If there exists k ∈ {1, 2, 3} for which D and V yield a.e. in Ω i one of the two assertions:

1. V k > 0 and V m,n ≥ 0 with D -1 k • V > 0 and D -1 m,n • V ≥ 0, 2. V k < 0 and V m,n ≥ 0 with D -1 k • V < 0 and D -1 m,n • V ≥ 0, (37)
where m, n are the two distinct elements of {1, 2, 3} \ {k}, then Ψ m,n i is injective in Ω i .

Proof. Let k ∈ {1, 2, 3} and m, n be the two distinct elements of {1, 2, 3} \ {k}. In view of (35), it follows that for all x = (x 1 , x 2 , x 3 ) and y = (y 1 , y 2 , y 3 ) elements of Ω i , we have

Ψ m,n i (x) = Ψ m,n i (y) ⇔          ψ i (x) = ψ i (y), v m 0i (x) = v m 0i (y), v n 0i (x) = v n 0i (y).
(38)

• Assertion 1. If V k > 0 and V m,n ≥ 0 a.e.
in Ω i , then since Lemma 3.4 applies, the two last equations in (38) imply, in view of (36), that each element of {x 1 -y 1 , x 2 -y 2 , x 3 -y 3 } is either null or admits the same sign as the other non-null element(s) of this set. Therefore, because according to Lemma 3.1,

ψ i (x) = ψ i (y) is equivalent to x 1 y 1 (D -1 ) 1 • V (η, x 2 , x 3 )dη + x 2 y 2 (D -1 ) 2 • V (y 1 , ξ, x 3 )dξ + x 3 y 3 (D -1 ) 3 • V (y 1 , y 2 , ζ)dζ = 0, (39)
and it holds D -1 j • V ≥ 0 a.e. in Ω i , for j = 1, 2, 3, then each of the three integrals in (39) is either null or admits the same sign as the other non-null integral(s). Hence, those three integrals are all null. Since D -1 k • V > 0, it follows that x k = y k . Afterwards, using x k = y k in the two last equations of (38) implies in view of (36) and as V k > 0 that x m = y m and x n = y n . That means x = y.

• Assertion 2. If V k < 0 and V m,n ≥ 0 a.e. in Ω i then, the two last equations in (38) imply in view of (36) that: 1. If x k -y k = 0, then each of the two terms x m -y m and x n -y n is either null or admits the opposite sign of

x k -y k . 2. If x k = y k , then x m = y m and x n = y n . Since D -1 k • V < 0 whereas D -1 m,n • V ≥ 0 a.e.
in Ω i , it follows that each of the three integrals in (39) is either null or admits the same sign as the other non-null integral(s). Thus, each of those integrals is null which implies that x k = y k . Because V k < 0, using x k = y k in the two last equations of (38) gives x m = y m and x n = y n .

Identifiability

We start this section by proving that within all suspected subdomain Ω i ⊂ Ω, there is no a single point source that can generate on the outflow boundary Γ i state and flux time measurements same as those generated by multiple distinct point sources. Then, we establish that the Dispersion-Current vector function Ψ m,n i introduced in (35) leads to localize in a unique manner the position of an unknown point source occurring in Ω i .

Theorem 4.1 Let i ∈ {1, . . . , I} and Ω i ⊂ Ω, where ∂Ω i = Γ i-1 ∪ Γ i F S ∪ Γ i F A ∪ Γ i .
We denote by u and u (P ) the solutions of the problem (1)-( 5) with F (x, t) = λ(t)δ(x -S) and F (x, t) = P p=1 λ p (t)δ(x -S p ) in Ω i × (0, T ), where P ∈ IN * , S p=1,...,P , S are interior points in Ω i and λ p=1,...,P , λ are functions of L 2 (0, T ) fulfilling [START_REF] Quiroga | Adjoint method for a tumor invasion Pde-constrained optimization problem in 2D using adaptive finite element method[END_REF].

   u = u (P ) on Γ i × (0, T ), D∇u • ν = D∇u (P ) • ν on Γ i × (0, T ). =⇒ P = 1 and S 1 = S. (40)
Proof. Let w = u -u (P ) in Ω i × (0, T ) and Ω 0 i = Ω i \ {S, S p=1,...,P }. In view of the boundary conditions in (40) and the two main equations ( 1)-( 2), it follows that w solves

       L[w](x, t) = λ(t)δ(x -S) - P p=1 λ p (t)δ(x -S p ) in Ω i × (0, T ), w = D∇w • ν = 0 on Γ i × (0, T ). (41) Let ω i be a connected open subset of Ω such that: ∅ = ω i ∩ Ω i ⊂ Ω 0 i , ω i ∩ Γ i = ∅ and ∅ = ω i ∩ Ω i+1 ⊂ Ω 0 i+1
, where Ω 0 i+1 is an open subset of Ω i+1 that doesn't contain any source position. Since w fulfills (41), then its extension by zero in ω i ∩ Ω i+1 solves in ω i × (0, T ) the homogeneous parabolic equation obtained from replacing by zero the right hand side of the first equation in (41). Thus, by applying twice the unique continuation Theorem [START_REF] Lin | A uniqueness theorem for parabolic equations[END_REF], we get w = 0 in ω i × (0, T ) and then, w = 0 in Ω 0 i × (0, T ). That implies w is either null a.e. in Ω i × (0, T ) or it's a linear combination of the dirac mass and some of its derivatives at the points S, S p=1,...,P . This second option it's not possible due to the regularity of w. Hence, by setting w = 0 a.e. in Ω i × (0, T ), it follows from the first equation of (41) that P = 1 and S 1 = S.

Besides, we establish the following constructive identifiability result. This result is said to be constructive since it sets guidelines that lead later on to localize in a unique manner the sought position of an unknown point source occurring in Ω i ⊂ Ω. Theorem 4.2 Let i ∈ {1, . . . , I}, Ω i ⊂ Ω and u be the solution of the problem (1)-( 5) with F (x, t) = λ(t)δ(x -S) in Ω i × (0, T ), where S is an interior point in Ω i . Provided 1. D and V are such that Theorem 3.5 applies in Ω i , 2. λ ∈ L 2 (0, T ) and fulfills [START_REF] Quiroga | Adjoint method for a tumor invasion Pde-constrained optimization problem in 2D using adaptive finite element method[END_REF], boundary records of u, D∇u • ν on Γ i-1 ∪ Γ i × 0, T and of uDν on Γ i F A × (0, T ) determine uniquely the unknown source elements S and λ = T 0 λ(t)e Rt dt. Proof. Let u (k=1,2) be the solution of the problem (1)- [START_REF] Hawkins-Daarud | Numerical simulation of a thermodynamically consistent four-species tumor growth model[END_REF] with

F in Ω i × (0, T ) equal to F (k) (x, t) = λ (k) (t)δ(x -S (k) )
, where λ (k) ∈ L 2 (0, T ) fulfilling ( 6) and S (k) is an interior point in Ω i . We denote by w = u (2) -u (1) . Then, from assuming

w = D∇w • ν = 0 on Γ i-1 ∪ Γ i × 0, T and wDν = 0 on Γ i F A × (0, T ), ( 42
)
it follows that the variable w solves the system:

               L[w](x, t) = λ (2) (t)δ(x -S (2) ) -λ (1) (t)δ(x -S (1) ) in Ω i × (0, T ), w(•, 0) = 0 in Ω i , w = 0 on Γ i-1 × (0, T ), D∇w • ν = 0 on Γ i F S ∪ Γ i F A ∪ Γ i × (0, T ). ( 43
)
From the last equation in (42) and according to [START_REF] Oelkers | Physical and chemical properties of rocks and fluids for chemical mass transport calculations[END_REF], it follows that all of the adjoint functions v i obtained in (34) solve the desired system [START_REF] Badia | Identification of a point source in a linear advection-dispersion-reaction equation: application to a pollution source problem[END_REF], where its last boundary condition is replaced by w D∇v i + v i V ) • ν = 0 on Γ i F A × (0, T ). Hence, as in ( 11)-( 12), from multiplying the first equation in (43) by an adjoint function v i (x, t) = e Rt v 0i (x) of (34) and integrating by parts over Ω i × (0, T ) using Green's formula, we get λ(2) v 0i (S (2) ) -λ(1) v 0i (S (1)

) = A v 0i i , ∀v 0i (x) ∈ 1, e ψ i (x) , v 1 0i (x), v 2 0i (x), v 3 0i (x) , (44)
where λ(k) =

T 0 e Rt λ (k) (t)dt, for k = 1, 2 and

A v 0i i = e RT Ω i w(•, T )v 0i + (Γ i-1 ∪Γ i )×(0,T ) e Rt w D∇v 0i + v 0i V -v 0i D∇w • ν. (45)
Moreover, using (42) in ( 44)-(45), gives: For all

v 0i (x) ∈ 1, e ψ i (x) , v 1 0i (x), v 2 0i (x), v 3 0i (x) , λ(2) v 0i (S (2) ) -λ(1) v 0i (S (1) ) = e RT Ω i w(•, T )v 0i . (46)
Besides, according to (8), w solves in Ω I+1 × (0, T ) a system similar to (43), where the first equation becomes homogeneous. That implies w = 0 in Ω I+1 × (0, T ) which, in view of (6) and using the problem satisfied by w in Ω × (T 0 , T ), gives from applying the unique continuation Theorem [START_REF] Lin | A uniqueness theorem for parabolic equations[END_REF] that w(•, T ) = 0 in Ω. Hence, (46

) leads to λ(2) v 0i (S (2) ) = λ(1) v 0i (S (1) ), ∀v 0i (x) ∈ {1, e ψ i (x) , v 1 0i (x), v 2 0i (x), v 3 0i (x)}. ( 47
)
Using v 0i = 1 in (47) gives λ(2) = λ(1) . Afterwards, from (47) and (35), we get

Ψ m,n i (S (2) ) = Ψ m,n i (S (1) ) =⇒ S (2) = S (1) . ( 48
)
The implication in (48) is obtained from applying Theorem 3.5.

Detection-Identification of sources

Let i ∈ {1, . . . , I}, Ω i be a suspected section of Ω and N i ∈ {0, . . . , N } be the number of point sources occurring in Ω i i.e., S n ∈ Ω i , for n = n i , . . . , n i + N i -1, where n i ∈ {1, . . . , N }. Thus, in the problem (1)-( 5), F (x, t) 11)-( 12), from multiplying the equations ( 1)-( 2) by an adjoint function v i (x, t) = e Rt v 0i (x) from (34) and integrating by parts in Ω i ×(0, T ) using Green's formula, we obtain: For all v 0i (x) ∈ {1,

= n i +N i -1 n=n i λ n (t)δ(x-S n ) in Ω i ×(0, T ). Therefore, likewise (
e ψ i (x) , v 1 0i (x), v 2 0i (x), v 3 0i (x)}, A v 0i i =        n i +N i -1 n=n i λn v 0i (S n ) if N i ≥ 1, 0 if N i = 0, (49)
where λn =

T 0 e Rt λ n (t)dt and the coefficient

A v 0i i = e RT Ω i u(•, T )v 0i + ∂Ω i ×(0,T ) e Rt u D∇v 0i + v 0i V -v 0i D∇u • ν. (50)
Afterwards, using in (49)-(50) the free space adjoint function i.e., v 0i = 1 leads to

A 1 i = e RT Ω i u(•, T ) + (Γ i-1 ∪Γ i )×(0,T ) e Rt uV -D∇u • ν =        n i +N i -1 n=n i λn if N i ≥ 1, 0 if N i = 0. (51)
Hence, employing the records of u, D∇u•ν taken on the intersections (Γ i-1 ∪Γ i )×(0, T ) to compute the coefficient A 1 i in (51) enables to detect the presence of all active sources within Ω i . Indeed, if A 1 i = 0 then, its value corresponds to the total suspended amount discharged by all sources occurring in Ω i . Moreover, for example, provided the reaction coefficient R ≥ 0 and n i +N i -1 n=n i λ n (t) ≥ 0 a.e. in (0, T ), it follows from (51) that

e -RT A 1 i ≤ n i +N i -1 n=n i T 0 λ n (t)dt ≤ A 1 i . (52)
Remark 5.1 Since in practice R is usually small, (52) could provide an approximation/estimation of the total amount discharged by all sources occurring in Ω i without having to identify a priori the historic in (0, T ) of each time-dependent intensity λ n .

Localization of a detected source

Assume the detection coefficient A 1 i = 0 and there exists one source occurring in Ω i i.e., N i = 1 and thus, F (x, t) = λ n i (t)δ(x -S n i ) in Ω i × (0, T ). From (49)-(50), it follows that the two unknown elements S n i and λ n i defining the occurring source are subject to:

               λn i = A 1 i , e ψ i (Sn i ) = A e ψ i i A 1 i , v 0i (S n i ) = A v 0i i A 1 i , for = 1, 2, 3. =⇒ Ψ m,n i (S n i ) = 1 A 1 i      A e ψ i i A v m 0i i A v n 0i i      , (53)
where Ψ m,n i is the Dispersion-Current function in (35) defined from the two distinct integers m and n selected in {1, 2, 3} such that Theorem 3.5 yields injectivity of Ψ m,n i in Ω i . Therefore, the unknown source position S n i is the unique point of Ω i that solves the last equation in (53). Besides, according to (50) and in view of (34), we get

• A 1 i = e RT Ω i u(•, T ) + (Γ i-1 ∪Γ i )×(0,T )
e Rt uV -D∇u • ν,

• A e ψ i i = e RT Ω i u(•, T )e ψ i - (Γ i-1 ∪Γ i )×(0,T ) e Rt+ψ i D∇u • ν, • A v =1,2,3 0i i = e RT Ω i u(•, T )v 0i + (Γ i-1 ∪Γ i )×(0,T ) e Rt u e g i DV ⊥ + v 0i V -v 0i D∇u • ν + Γ i F A ×(0,T ) e Rt+g i V ⊥ • uDν. (54)
From computing the coefficients in (54), we localize the unknown source position S n i as the unique point within the suspected section Ω i that solves the last equation in (53).

Determination of the final state u(•, T )

For the computation of the coefficients in (54), we propose the following two different ways to determine the unknown data u(•, T ):

• First way. Least squares: Since from (8) the subdomain Ω I+1 is a "No source region", we use the state records taken on Γ I × (0, T ) to solve the forward problem satisfied by the state u in Ω I+1 × (0, T ). Then, because all sources become inactive in (T 0 , T ), we identify u(•, T ) as the final condition that leads the solution of the Cauchy problem satisfied by u in Ω × (T 0 , T ) to fit the forward problem solution in Ω I+1 × (T * , T ), where T * ∈ (T 0 , T ).

• Second way. Data assimilation: For i = 1, . . . , I, use the state records d i-1 on Γ i-1 × (0, T ) and d i on Γ i × (0, T ) to determine an approximation û(•, T ) in Ω i of u(•, T ) from solving the following system:

                     L[û i ](x, t) = 0 in Ω i × (0, T ), ûi (•, 0) = 0 in Ω i , ûi = d i-1 on Γ i-1 × (0, T ), ûi = d i on Γ i × (0, T ), D∇û i • ν = 0 on Γ i F S ∪ Γ i F A × (0, T ). ( 55
)
Then, we establish the following Lemma on how does ûi (•, T ) approximate u(•, T ) in Ω i :

Lemma 5.2 Provided (6) holds, there exists 0 < α ∈ IR such that the state ûi solution of the system (55) fulfills: For all t ∈ (T 0 , T ),

ûi (•, t) -u(•, t) 2 L 2 (Ω i ) ≤ ûi (•, T 0 ) -u(•, T 0 ) 2 L 2 (Ω i ) e -α(t-T 0 ) , ( 56
)
where T 0 is the final active time introduced in (6) and u is the solution of ( 1)- [START_REF] Hawkins-Daarud | Numerical simulation of a thermodynamically consistent four-species tumor growth model[END_REF].

Proof. See the Appendix.

For the clarity of our presentation, we summarize in the following algorithm the main steps defining the Detection-Identification method developed in the present paper:

Algorithm. Detection-Identification method

For i = 1 to I do • Compute the detection coefficient A 1 i from (54). • If | A 1 i | ≤ ε 0 then,

Application to three types of flow

We apply the developed Detection-Identification method to the following three different types of flow crossing a suspected subdomain Ω i ⊂ Ω: 1. Flow defined by a mean velocity vector 2. Flow defined by a unidirectional velocity field 3. Flow defined by a bidirectional velocity field. We employ Bear's hydrodynamic dispersion tensor, see [START_REF] Bear | Introduction to Modeling of Transport Phenomena in Porous Media[END_REF][START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Oelkers | Physical and chemical properties of rocks and fluids for chemical mass transport calculations[END_REF]:

D = D M + α T V 2 I + α L -α T V 2 V • V , (57)
where D M > 0 is a real number that represents the molecular diffusion, I is the 3 × 3 identity matrix and 0 < α T < α L are two real numbers that represent the longitudinal and transverse dispersivities of the porous medium. From (57), it follows that

• DV = D M + α L V 2 V ⇔ D -1 V = 1 D M + α L V 2 V, • DV ⊥ = D M + α T V 2 V ⊥ , for = 1, 2, 3. (58)
Therefore, the scalar potential ψ i in ( 19) is given by: For all x = (x 1 , x 2 , x 3 ) ∈ Ω i ,

ψ i (x) = - x 1 a i V 1 D M + α L V 2 (η, x 2 , x 3 )dη - x 2 b i V 2 D M + α L V 2 (a i , ξ, x 3 )dξ - x 3 c i V 3 D M + α L V 2 (a i , b i , ζ)dζ.
(59)

Flow of mean velocity vector

In

Ω i ⊂ Ω, assume V = V 1 , V 2 , V 3
to be a mean velocity vector and D to be the associated dispersion tensor in (57). Then, ( 17) is satisfied and from (59), we get

ψ i (x) = - 1 D M + α L V 2 V • x -A i , ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i , (60)
where

A i = (a i , b i , c i ) ∈ Ω i .
Besides, according to Remark 3.2, we define v =1,2,3 0i from solving ∇v 0i = V ⊥ in Ω i . Thus, setting g i = 0 in ( 31)- [START_REF] Nara | A projective method for an inverse source problem of the Poisson equation[END_REF] gives: For = 1, 2, 3

v 0i (x) = V ⊥ • x -A i , ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (61)
Provided the mean velocity vector V fulfills in Ω i , for example, V 1 > 0, V 2 ≥ 0 and V 3 ≥ 0, it follows from (58) that (D -1 ) 1 V > 0 and (D -1 ) 2,3 V ≥ 0 and thus, the first condition of Theorem 3.5 is fulfilled for k = 1, m = 2 and n = 3. Therefore, from (35), the function

Ψ 2,3 i : Ω i -→ IR 3 , x → Ψ 2,3 i (x) =      e ψ i (x) v 2 0i (x) v 3 0i (x)      , (62)
is injective in Ω i . Afterwards, assuming the detection coefficient A 1 i = 0 and there exists one unknown source λ n i (t)δ(x -S n i ) occurring in Ω i , it comes from (53)-(54) that the sought source position S n i is the unique point within Ω i that satisfies

Ψ 2,3 i (S n i ) = 1 A 1 i      A e ψ i i A v 2 0i i A v 3 0i i      , (63)
which, in view of (60)-(61), can be written under the following matrix form:

     V 1 V 2 V 3 -V 3 0 V 1 -V 2 V 1 0      S n i -A i =          -D M + α L V 2 ln A e ψ i i A 1 i A v 2 0i i A 1 i A v 3 0i i A 1 i          . ( 64
)
The determinant of the 3 × 3 matrix involved in the linear system (64) is -V 1 V 2 2 < 0. That confirms the result announced by Theorem 3.5 on the injectivity of the Dispersion-Current function Ψ 2,3 i in (62). Therefore, the unknown position S n i defining the detected source is the unique point within Ω i that solves the linear system in (64).

Flow of unidirectional velocity field

Assume within Ω i ⊂ Ω, there exists two functions ϕ and ζ such that the velocity field

V (x) = V 1 (x), 0, 0 ⊥ , where V 1 (x) = ϕ(x 2 )ζ(x 3 ) > 0, for all x = (x 1 , x 2 , x 3 ) ∈ Ω i .
Since the molecular diffusion D M is usually small enough [START_REF] Bear | Introduction to Modeling of Transport Phenomena in Porous Media[END_REF][START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF], for the simplicity of our presentation we set D M = 0 in (58). That leads to

D -1 V = 1 α L 1, 0, 0 in Ω i =⇒ -→ rot D -1 V = 0 in Ω i . (65)
Thus, the condition (17) on the existence of the scalar potential ψ i defined from ( 16) is well fulfilled. Moreover, from using D M = 0 in (59), we obtain

ψ i (x) = - 1 α L (x 1 -a i ), ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (66)
Since V 1 > 0 and V 2 = V 3 = 0 in Ω i , it follows that the first condition of Theorem 3.5 is satisfied for k = 1, m = 2 and n = 3. Therefore, the function Ψ 2,3 i defined from (35) is injective in Ω i . To define Ψ 2,3 i , we determine the two adjoint functions v 2 0i and v 3 0i .

According to (31), v 2 0i is defined by ∇v 2 0i = e g 2 i V 2 ⊥ in Ω i , where V 2 ⊥ = 0, 0, V 1 and g 2 i solves [START_REF] Lions | Some methods in the mathematical analysis of systems and their control[END_REF]. Moreover, using (58) in the linear system [START_REF] Lions | Some methods in the mathematical analysis of systems and their control[END_REF] implies that g 2 i is subject to:

∇g 2 i = - 1 V 1 G 2 , where G 2 =      0 ∂ x 2 V 1 2∂ x 3 V 1      . ( 67
) Because V 1 (x) = ϕ(x 2 )ζ(x 3 ) in Ω i , it follows that -→ rot 1 V 1 G 2 = 0
in Ω i and thus, we get

g 2 i (x) = -ln |ϕ(x 2 )| -2ln |ζ(x 3 )| = ln 1 |ϕ(x 2 )|ζ 2 (x 3 ) , ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (68)
That leads, in view of (32), to determine v 2 0i as follows:

v 2 0i (x) = x 3 c i e g 2 i V 1 (b i , ξ)dξ = sign ϕ(b i ) x 3 c i 1 ζ(ξ) dξ, ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (69)
Besides, v 3 0i is defined from (31) by ∇v 3 0i = e g 3 i V 3 ⊥ in Ω i , where V 3 ⊥ = 0, V 1 , 0 and g 3 i solves [START_REF] Linfield | The enhanced stream water quality models QUAL2E and QUAL2E-UNCAS: Documentation and user manual[END_REF]. From employing (58) in the linear system [START_REF] Linfield | The enhanced stream water quality models QUAL2E and QUAL2E-UNCAS: Documentation and user manual[END_REF], it comes that g 3 i solves in Ω i

∇g 3 i = - 1 V 1 G 3 , where G 3 =      0 2∂ x 2 V 1 ∂ x 3 V 1      . ( 70
) Since V 1 (x) = ϕ(x 2 )ζ(x 3 ) in Ω i , it follows that -→ rot 1 V 1 G 3 = 0 in Ω i .
From (70), we get

g 3 i (x) = -2ln |ϕ(x 2 )| -ln |ζ(x 3 )| = ln 1 ϕ 2 (x 2 )|ζ(x 3 )| , ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (71)
Afterwards, from (32), we find

v 3 0i (x) = x 2 b i e g 3 i V 1 (η, c i )dη = sign ζ(c i ) x 2 b i 1 ϕ(η) dη, ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (72)
Hence, given (a i , b i , c i ) ∈ Ω i and in view (66)-(72), the function Ψ 2,3 i is defined by

Ψ 2,3 i : Ω i -→ IR 3 , x → Ψ 2,3 i (x) =           exp -1 α L (x 1 -a i ) sign ϕ(b i ) x 3 c i 1 ζ(ξ) dξ sign ζ(c i ) x 2 b i 1 ϕ(η) dη           . (73)
From ( 73) and since the two functions ϕ and ζ are such that ϕ(x 2 )ζ(x 3 ) > 0 in Ω i , we verify readily the injectivity of Ψ 2,3 i in Ω i announced by Theorem 3.5. Thus, assuming the detection coefficient A 1 i = 0 and one unknown source λ n i (t)δ(x -S n i ) is occurring in Ω i , it follows that the sought source position S n i is the unique point within Ω i that solves the last equation in (53), where Ψ m,n i = Ψ 2,3 i obtained in (73).

Flow of bidirectional velocity field

The third type of flow crossing a suspected section Ω i ⊂ Ω that we consider for the application of the developed Detection-Identification method is a flow defined by a velocity field: V (x) = V 0 (x 3 ) 1, β, 0 , where the function V 0 (x 3 ) > 0 for all x = (x 1 , x 2 , x 3 ) ∈ Ω i and β ≥ 0 is a real number. From setting D M = 0 in (58), it follows that

• D -1 V = 1 α L 1 + β 2 1, β, 0 in Ω i =⇒ -→ rot D -1 V = 0 in Ω i , • DV ⊥ = α T 1 + β 2 V 0 V ⊥ in Ω i , for = 1, 2, 3. (74)
The second equation in (74) implies that the condition [START_REF] Hamdi | Inverse source problem in a one-dimensional evolution linear transport equation with spatially varying coefficients: application to surface water pollution[END_REF] on the existence of the scalar potential ψ i defined from ( 16) is well fulfilled. Moreover, using D M = 0 in (59) gives

ψ i (x) = - 1 α L 1 + β 2 x 1 -a i + β(x 2 -b i ) , ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (75)
Besides, in view of (74), since the function V 0 > 0 and the real number β ≥ 0, the first condition of Theorem 3.5 is satisfied for k = 1, m = 2 and n = 3. Thus, the function Ψ 2,3 i in (35) is injective in Ω i . To define Ψ 2,3 i , we determine the two functions v 2 0i and v 3 0i . From [START_REF] Machado | A new One-Shot Pointwise Source Reconstruction Method[END_REF], the adjoint function v 2 0i is defined by ∇v 2 0i = e g 2 i V 2 ⊥ in Ω i , where V 2 ⊥ = 0, 0, V 0 and g 2 i satisfies [START_REF] Lions | Some methods in the mathematical analysis of systems and their control[END_REF]. Moreover, using (74) in the linear system (28) leads to

∇g 2 i =       0 0 -2 V 0 V 0       =⇒ g 2 i (x) = ln 1 V 2 0 (x 3 ) , ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . ( 76
)
Furthermore, from replacing in (32) g 2 i by its value obtained in (76), we find

v 2 0i (x) = x 3 c i 1 V 0 (ξ) dξ, ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . ( 77
)
The adjoint function v 3 0i is defined from (31) by ∇v 3 0i = e g 3 i V 3 ⊥ in Ω i , where V 3 ⊥ = -βV 0 , V 0 , 0 and g 3 i solves [START_REF] Linfield | The enhanced stream water quality models QUAL2E and QUAL2E-UNCAS: Documentation and user manual[END_REF]. Then, using (74) in the linear system (30) leads to

∇g 3 i =       0 0 - V 0 V 0       =⇒ g 3 i (x) = ln 1 V 0 (x 3 ) , ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . ( 78
)
Afterwards, from subtituting in [START_REF] Nara | A projective method for an inverse source problem of the Poisson equation[END_REF] the function g 3 i by its value in (78), we get 75), ( 77) and ( 79) as follows:

v 3 0i (x) = -β(x 1 -a i ) + (x 2 -b i ), ∀x = (x 1 , x 2 , x 3 ) ∈ Ω i . (79) Therefore, given (a i , b i , c i ) ∈ Ω i , Ψ 2,3 i is defined in Ω i from (
Ψ 2,3 i : Ω i -→ IR 3 , x → Ψ 2,3 i (x) =         exp - 1 α L 1 + β 2 x 1 -a i + β(x 2 -b i ) x 3 c i 1 V 0 (ξ) dξ -β(x 1 -a i ) + (x 2 -b i )         . ( 80
)
Since the function V 0 > 0 in Ω i , we check readily from (80) the injectivity of Ψ 2,3 i in Ω i announced by Theorem 3.5. Hence, if the detection coefficient A 1 i = 0 due to the presence in Ω i of an unknown source λ n i (t)δ(x -S n i ), the sought source position S n i is the unique point within Ω i that solves the last equation in (53), where Ψ m,n i = Ψ 2,3 i given in (80).

Numerical experiments

We carry out numerical experiments in the case of a domain Ω defined by

Ω = x = (x 1 , x 2 , x 3 ) ∈ 0, 10 × 0, 1 × 0, 1 , (81)
wherein, we assume to be available the following I = 4 interfaces:

Γ i = x = (2i, x 2 , x 3 ) ∈ 0, 1 2 , for i = 1, . . . , I . (82)
We generate synthetic measurements on the interfaces Γ i=1,...,I × (0, T ) from solving the problem (1)-( 5), where D is Bear's dispersion tensor introduced in (57) for two different types of flow: Flow of mean velocity vector V = 0.5, 0, 0 and flow of unidirectional velocity field V (x) = x 2 (1 -x 2 )x 3 , 0, 0 . We used the coefficients R = 0, α T = 0.02, α L = 0.3 and D M = 0. Regarding the source F in (5), we considered N = 2 point sources located at S 1 = 1.2, 0.4, 0.3 ∈ Ω 1 , S 2 = 3.0, 0.6, 0.2 ∈ Ω 2 and loading in 0, T :

λ 1 (t) = sin 2π T t χ 0, T 2 (t) and λ 2 (t) = 5 2 sin 2π T (t - T 4 ) χ T 4 , 3T 4 (t). ( 83
)
Since along Γ F A = x = (x 1 , x 2 , 1) ∈ 0, 10 × 0, 1 we have ν = 0, 0, 1 and thus, V • ν = 0, it follows from (57) that Dν = α T V 2 ν on Γ F A . Therefore, in the carried out numerical experiments, we neglected the term Dν which means we do not use state measurements on Γ F A × 0, T in the computation of the coefficients v 0i from (54).

We computed from (54) the detection coefficient A 1 i associated to all suspected sections Ω i=1,...,I . If A 1 i = 0, we localize the unknown position defining the detected source occurring in Ω i from: 1. Solving the linear system in (64), for the flow of mean velocity vector. 2. Solving the last equation in ( 53), where Ψ m,n i = Ψ 2,3 i given in (73), for the flow defined by the unidirectional velocity field. Besides, since we use R = 0, it follows from (53) that if there exists a source

λ n i (t)δ(x -S n i ) occurring in Ω i , then A 1 i = T 0 λ n i (t)dt.
However, if there is no source occurring in Ω i , then A 1 i = 0. We set T = 20 and using (83), we calculate the exact values of A 1 i in all suspected sections Ω i=1,...,I :

A 1 1 = T π ≈ 6.37, A 1 2 = 5T 2π ≈ 15, 92 and A 1 3 = A 1 4 = 0. ( 84
)
As far as the unknown data u(•, T ) involved in ( 54) is concerned, we used data assimilation to determine an approximation ûi (•, T ) of u(•, T ) in Ω i . To this end, we solved the forward problem (55) whose, according to Lemma 5.2, the solution ûi (•, t) converges to u(•, T ) in Ω i when t tends to T . We start by presenting the behaviour of the detection coefficient 54), where we replaced u(•, T ) by ûi (•, t) for t ∈ (0, T ). The results in Figure 1 show that for the flow of mean velocity vector, the detection coefficient A 1 i calculated from (54) converges quickly to its exact values given in (84) for all suspected sections Ω i=1,...,I . For the flow of velocity field, the calculated A 1 i=1,2 in the two sections Ω i=1,2 , wherein there is an occurring source, converge to their exact values when t is sufficiently close to T = 20. This observation could be explained, in view of Lemma 5.2, by the convergence rate of ûi (•, t) to the unknown final state u(•, T ).

A 1 i computed from (
For the two considered types of flow, we computed the coefficients in (54) and localized the positions of the detected unknown sources. In the sequel, we present the obtained results, where the dashed lines give the exact coordinates of each sought source position: The analysis of the numerical results presented in Figures 23shows that for t sufficiently close to T , the used ûi (•, t) yields a good approximation in Ω i of u(•, T ) which leads the developed method to localize accurately the unknown positions of all detected sources.

Baby Example

Here, we illustrate in the one dimensional case of the underlined inverse source problem the non-identifiability of multiple unknown point sources from measurements taken only upstream and/or downstream all occurring point sources.

Let N ∈ IN * , 0 < ∈ IR and u be the solution of the following problem:

-Du (x) + V u (x) + Ru(x) = N n=1 λ n δ(x -S n ) in 0, and u(0) = u ( ) = 0. ( 85
)
Provided DV = 0 and V 2 + 4RD > 0, the characteristic equation -Dr 2 + V r + R = 0 admits the roots:

r 1 = V /(2D) 1+ 1 + 4RD/V 2 and r 2 = V /(2D) 1-1 + 4RD/V 2 . Let ψ(x) = -(e r 1
x -e r 2 x )/ D(r 1 -r 2 ) , ∀x ∈ (0, ). The solution u of (85) is given by

u(x) = N n=1 λ n ψ(x -S n )H(x -S n ) - ψ ( -S n ) ψ ( ) ψ(x) , ∀x ∈ 0, , (86)
where H is the Heaviside function. We introduce the following observation operator:

M : IR * × (0, ) N -→ IR 2N , F = λ 1 , S 1 , . . . , λ N , S N → M [F ] = u(p 1 ), u(P 1 ), . . . , u(p N ), u(P N ) , (87)
where p i=1,...,N and P i=1,...,N are 2N distinct state measuring points within (0, ).

Theorem 8.1 Let D, V and R be real numbers such that DV = 0 and

V 2 + 4RD > 0. The 2N × 2N Jacobian matrix J M of the observation operator M [F ] in (85)-(87) fulfills 1. For N = 1: H(p 1 -S 1 ) = H(P 1 -S 1 ) ⇐⇒ det(J M ) = 0.
2. For N = 2: Let 0 < p 1 < P 1 < p 2 < P 2 < and, for example, S 1 < S 2 . We have

(i) Cardinality p i=1,2 , P i=1,2 ∩ S 1 , S 2 ≤ 1 =⇒ det(J M ) = 0, (ii) 0 < p 1 < S 1 < P 1 < p 2 < S 2 < P 2 < =⇒ det(J M ) = 0. (88)
Proof. See the Appendix.

Theorem 8.1 indicates that for N = 1 point source λδ(x -S) occurring in (0, ), two state measuring points p 1 , P 1 framing the source region i.e., 0 < p 1 < S < P 1 < yield local injectivity of the observation operator M [F ] in (85)-(87) and thus, local identifiability of the two source elements λ, S. However, for N = 2 point sources occurring in (0, ), the local identifiability of the four source elements holds only if in addition to the two state measuring points framing the sources region, two other state measuring points are set between the two source positions. Moreover, using similar computations, we prove that the two additional state measuring points that should separate the two occurring source positions can be replaced by one state and flux measuring point whereas the two measuring points p 1 and P 1 framing the sources region can be placed at the two boundary points i.e., p 1 = 0 and P 1 = if at these points both the state and its flux are known.

Conclusion, discussion and comparaison

In this paper, we developed appropriate adjoint functions that led to establish a constructive Detection-Identification method for solving the nonlinear inverse source problem of identifying multiple unknown time-dependent point sources occurring in 3D transport equations. In the literature, many authors have addressed similar inverse source problems in different PDEs, for instance [START_REF] Andrle | Identification of Moving Pointwise Sources in an Advection-Dispersion-Reaction Equation[END_REF][START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF][START_REF] Belgacem | Uniqueness for an ill-posed reaction-dispersion model. Application to organic pollution in stream-waters (One-dimensional model)[END_REF][START_REF] Hamdi | Detection-identification of multiple unknown time-dependent point sources in a 2D transport equation: Application to accidental pollution[END_REF][START_REF] Hasanov | Fourier collocation algorithm for identifying the spacewise-dependent source in the advection-diffusion equation from boundary data measurements[END_REF][START_REF] Isakov | Inverse Source Problems[END_REF][START_REF] Machado | A new One-Shot Pointwise Source Reconstruction Method[END_REF]. In these works, the developed identification approaches are mainly either iterative based on the minimisation of cost functions such as least squares and Kohn-Vogelius or quasi-direct such as the algebraic method especially used for elliptic equations [START_REF] Abdelaziz | Direct algorithms for solving some inverse source problems in 2D elliptic equations[END_REF][START_REF] Nara | A projective method for an inverse source problem of the Poisson equation[END_REF]. Besides, the underlined inverse source problem becomes more challenging in n = 2, 3 dimensions where the involved PDE admits an advection term. Indeed, this asymmetric term creates hydrodynamic dispersion as a consequence of molecular diffusion and mechanical dispersion caused by non-uniform velocities [START_REF] Bear | Introduction to Modeling of Transport Phenomena in Porous Media[END_REF][START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF]. Consequently, in such PDEs the effectiveness of a developed identification approach relies on how much it takes into account the impact of these physical phenomena on the used measurements: For example, in the case of an advection dominant flow, using measurements taken only upstream all occurring unknown sources do not yield identifiability as well as in the case of high dispersion coefficients, the signals emmitted by different sources might get rapidly intermixed and thus, without using measurements taken between the occurring sources these latest become indistinguishable. These assertions are illustrated in the one dimensional case by the Baby Example of section 8. Moreover, the authors in [START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF] used measurements taken on the whole boundary to address the identification of multiple time-dependent point sources in 2D advection-diffusion equation from minimising two different cost functions i.e., least squares and Kohn-Vogelius. They reported that the least squares approach doesn't enable to identify more than one active source whereas the Kohn-Vogelius approach identifies two sources only if they are well separated and the diffusion coefficient is very small in order to avoid rapid intermixing.

Therefore, comparing with identification approaches from the literature, it follows that:

1. The constructive Detection-Identification method developed in the present paper takes into account the velocity field effects on the flow nature, is quasi-direct and thus easy to implement and does not have to deal always with the classic issues of iterative approaches i.e., questions related to the choices of an initial iterate and a regularisation term as well as to ensure convergence.

2. The herein developed method is based on a direct detection procedure that determines whether within the monitored domain a single or rather multiple unknown point sources occur and provides a framing/approximation of the total amount discharged by all occurring sources. These results lead to i) Reducing the total identification cost: In the case of no significant unknown sources i.e., sources discharging a total amount smaller than a certain tolerance or the case where the immediate aim is rather to know how important is the total discharged amount in order to take prompt appropriate actions, we can choose to not go ahead with the identification of the remaining unknown elements defining the occurring sources. ii) Avoiding misinterpretations: In the case of multiple unknown point sources, the procedure detects the presence of more than one source and suggests to add measurements taken between the occurring sources in order to ensure their identifiability. However, as reported in [START_REF] Andrle | Identification of multiple moving pollution sources in surface waters or atmospheric media with boundary observations[END_REF], in this case usually the iterative approaches give identification results that could be inaccurate and thus, leading to misinterpretations.

3. For the identification approach developed in this paper, the construction of the appropriate adjoint functions defining the dispersion-current function Ψ m,n i relies on some conditions that should be fulfilled by the velocity field V and the dispersion tensor D. In practice, failing those conditions, one could employ rather approximations of V and D that ensure the applicability of the developed approach.

Appendix.

A. We start by establishing the proof of the results announced in Lemmas 3.1, 3.4 and 5.2.

I. Proof of Lemma 3.1: Let x = (x 1 , x 2 , x 3 ) and y = (y 1 , y 2 , y 3 ) be two elements of Ω i . In view of [START_REF] Hasanov | Fourier collocation algorithm for identifying the spacewise-dependent source in the advection-diffusion equation from boundary data measurements[END_REF], from writing ψ i (x) = ψ i (y) under the form ψ i (y) -ψ i (x) = 0 and treating in this latest equation separately each two integrals following the same axis, we obtain:

• The two integrals with respect to the variable η in the equation ψ i (y) -ψ i (x) = 0:

x 1 a i (D -1 ) 1 • V (η, x 2 , x 3 )dη - y 1 a i (D -1 ) 1 • V (η, y 2 , y 3 )dη = x 1 y 1 (D -1 ) 1 • V (η, x 2 , x 3 )dη + y 1 a i (D -1 ) 1 • V (η, x 2 , x 3 )dη - y 1 a i (D -1 ) 1 • V (η, y 2 , y 3 )dη. ( 89
)
The first integral in the right hand side of the equation ( 89) is obtained from using Chasles's relation by introducing y 1 on the integral over (a i , x 1 ) in the left hand side of (89). Moreover, by adding and substracting the following integral:

y 1 a i (D -1 ) 1 • V (η, y 2 , x 3 )dη, (90)
to the right hand side of the equation in (89), we get

x 1 a i (D -1 ) 1 • V (η, x 2 , x 3 )dη - y 1 a i (D -1 ) 1 • V (η, y 2 , y 3 )dη = x 1 y 1 (D -1 ) 1 • V (η, x 2 , x 3 )dη + y 1 a i (D -1 ) 1 • V (η, ξ, x 3 ) ξ=x 2 ξ=y 2 dη + y 1 a i (D -1 ) 1 • V (η, y 2 , ζ) ζ=x 3 ζ=y 3 dη. (91)
Thus, using Fubini's Theorem in the two last integrals of the right hand side in (91) gives

x 1 a i (D -1 ) 1 • V (η, x 2 , x 3 )dη - y 1 a i (D -1 ) 1 • V (η, y 2 , y 3 )dη = x 1 y 1 (D -1 ) 1 • V (η, x 2 , x 3 )dη + x 2 y 2 y 1 a i ∂ ξ (D -1 ) 1 • V (η, ξ, x 3 )dηdξ + x 3 y 3 y 1 a i ∂ ζ (D -1 ) 1 • V (η, y 2 , ζ)dηdζ. (92)
In addition, since

-→ rot D -1 V = 0 in Ω i , it follows from (18) that * ∂ ξ (D -1 ) 1 • V (η, ξ, x 3 ) = ∂ η (D -1 ) 2 • V (η, ξ, x 3 ), * ∂ ζ (D -1 ) 1 • V (η, y 2 , ζ) = ∂ η (D -1 ) 3 • V (η, y 2 , ζ). (93)
Hence, from replacing the partial derivatives in (92) by their values in (93), we obtain

x 1 a i (D -1 ) 1 • V (η, x 2 , x 3 )dη - y 1 a i (D -1 ) 1 • V (η, y 2 , y 3 )dη = x 1 y 1 (D -1 ) 1 • V (η, x 2 , x 3 )dη + x 2 y 2 (D -1 ) 2 • V (η, ξ, x 3 ) η=y 1 η=a i dξ + x 3 y 3 (D -1 ) 3 • V (η, y 2 , ζ) η=y 1 η=a i dζ. (94)
• The two integrals with respect to the variable ξ in the equation ψ i (y) -ψ i (x) = 0:

x 2 b i (D -1 ) 2 • V (a i , ξ, x 3 )dξ - y 2 b i (D -1 ) 2 • V (a i , ξ, y 3 )dξ = x 2 y 2 (D -1 ) 2 • V (a i , ξ, x 3 )dξ + y 2 b i (D -1 ) 2 • V (a i , ξ, ζ) ζ=x 3 ζ=y 3 dξ, = x 2 y 2 (D -1 ) 2 • V (a i , ξ, x 3 )dξ + x 3 y 3 y 2 b i ∂ ζ (D -1 ) 2 • V (a i , ξ, ζ)dξdζ. (95)
The right hand side of the first equality in (95) is obtained from its left hand side term by introducing, according to Chasles's relation, y 2 on the integral over (b i , x 2 ). Moreover, since -→ rot D -1 V = 0 in Ω i and in view of [START_REF] Hamdi | Detection-identification of multiple unknown time-dependent point sources in a 2D transport equation: Application to accidental pollution[END_REF], it follows that

∂ ζ (D -1 ) 2 • V (a i , ξ, ζ) = ∂ ξ (D -1 ) 3 • V (a i , ξ, ζ). (96)
Thus, using (96) to replace the partial derivative in the last integral of (95) leads to

x 2 b i (D -1 ) 2 • V (a i , ξ, x 3 )dξ - y 2 b i (D -1 ) 2 • V (a i , ξ, y 3 )dξ = x 2 y 2 (D -1 ) 2 • V (a i , ξ, x 3 )dξ + x 3 y 3 (D -1 ) 3 • V (a i , ξ, ζ) ξ=y 2 ξ=b i dζ. (97)
• The two integrals with respect to the variable ζ in the equation ψ i (y) -ψ i (x) = 0:

x 3 c i (D -1 ) 3 • V (a i , b i , ζ)dζ - y 3 c i (D -1 ) 3 • V (a i , b i , ζ)dζ = x 3 y 3 (D -1 ) 3 • V (a i , b i , ζ)dζ.( 98
)
Therefore, since ψ i (y) -ψ i (x) is equal to the sum of the left hand side terms in the three equations ( 94), ( 97) and (98) then, from setting the sum of the right hand side terms of those equations to null, we find the result announced in Lemma 3.1.

II. Proof of Lemma 3.4: Let x = (x 1 , x 2 , x 3 ) and y = (y 1 , y 2 , y 3 ) be two elements of Ω i . In view of [START_REF] Nara | A projective method for an inverse source problem of the Poisson equation[END_REF], it follows that:

• The first equation v 1 0i (x) = v 1 0i (y) is equivalent to x 3 y 3 e g 1 i V 2 (a i , b i , ξ)dξ = x 2 b i e g 1 i V 3 (a i , η, x 3 )dη - y 2 b i e g 1 i V 3 (a i , η, y 3 )dη, = x 2 y 2 e g 1 i V 3 (a i , η, x 3 )dη + y 2 b i e g 1 i V 3 (a i , η, ξ) ξ=x 3 ξ=y 3 dη, = x 2 y 2 e g 1 i V 3 (a i , η, x 3 )dη + x 3 y 3 y 2 b i ∂ ξ e g 1 i V 3 (a i , η, ξ) dηdξ, = x 2 y 2
e g 1 i V 3 (a i , η, x 3 )dη -

x 3 y 3 y 2 b i ∂ η e g 1 i V 2 (a i , η, ξ) dη dξ. (99)
The second equality in (99) is obtained using Chasles's relation by adding the point y 2 on the interval of integration (b i , x 2 ). The third equality is found using Fubini's theorem by integrating first with respect to the variable η then, with respect to ξ. The Fourth equality is obtained from the compatibility condition -→ rot e

g 1 i V 1 ⊥ = 0 in Ω i which implies that ∂ x 3 e g 1 i V 3 = -∂ x 2 e g 1 i V 2 in Ω i .
Therefore, by computing the integral over dη in the right hand side of the last equation of (99) and then, cancelling out the two same terms occurring in both sides of this equation, we obtain

v 1 0i (x) = v 1 0i (y) ⇔ x 2 y 2 e g 1 i V 3 (a i , η, x 3 )dη = x 3 y 3 e g 1 i V 2 (a i , y 2 , ξ)dξ. (100)
• The second equation v 2 0i (x) = v 2 0i (y) is equivalent to

x 3 y 3 e g 2 i V 1 (a i , b i , ξ)dξ = x 1 a i e g 2 i V 3 (η, b i , x 3 )dη - y 1 a i e g 2 i V 3 (η, b i , y 3 )dη, = x 1 y 1 e g 2 i V 3 (η, b i , x 3 )dη + y 1 a i e g 2 i V 3 (η, b i , ξ) ξ=x 3 ξ=y 3 dη, = x 1 y 1 e g 2 i V 3 (η, b i , x 3 )dη + x 3 y 3 y 1 a i ∂ ξ e g 2 i V 3 (η, b i , ξ) dηdξ, = x 1 y 1 e g 2 i V 3 (η, b i , x 3 )dη - x 3 y 3 y 1 a i ∂ η e g 2 i V 1 (η, b i , ξ) dη dξ. (101)
The second equality in (101) is obtained using Chasles's relation by adding the point y 1 on the interval of integration (a i , x 1 ). The third equality is found using Fubini's theorem by integrating first with respect to the variable η then, with respect to ξ. The Fourth equality is obtained from the compatibility condition -→ rot e

g 2 i V 2 ⊥ = 0 in Ω i which implies that ∂ x 3 e g 2 i V 3 = -∂ x 1 e g 2 i V 1 in Ω i .
Hence, by computing the integral over dη in the right hand side of the last equation of (101) and then, cancelling out the two same terms occurring in both sides of this equation, we get

v 2 0i (x) = v 2 0i (y) ⇔ x 1 y 1 e g 2 i V 3 (η, b i , x 3 )dη = x 3 y 3 e g 2 i V 1 (y 1 , b i , ξ)dξ. (102)
• The third equation v

3 0i (x) = v 3 0i (y) is equivalent to x 2 y 2 e g 3 i V 1 (a i , ξ, c i )dξ = x 1 a i e g 3 i V 2 (η, x 2 , c i )dη - y 1 a i e g 3 i V 2 (η, y 2 , c i )dη, = x 1 y 1 e g 3 i V 2 (η, x 2 , c i )dη + y 1 a i e g 3 i V 2 (η, ξ, c i ) ξ=x 2 ξ=y 2 dη, = x 1 y 1 e g 3 i V 2 (η, x 2 , c i )dη + x 2 y 2 y 1 a i ∂ ξ e g 3 i V 2 (η, ξ, c i ) dηdξ, = x 1 y 1 e g 3 i V 2 (η, x 2 , c i )dη - x 2 y 2 y 1 a i ∂ η e g 3 i V 1 (η, ξ, c i ) dη dξ. (103)
The second equality in (103) is obtained using Chasles's relation by adding the point y 1 on the interval of integration (a i , x 1 ). The third equality is found using Fubini's theorem by integrating first with respect to the variable η then, with respect to ξ. The Fourth equality is obtained from the compatibility condition -→ rot e

g 3 i V 3 ⊥ = 0 in Ω i which implies that ∂ x 2 e g 3 i V 2 = -∂ x 1 e g 3 i V 1 in Ω i .
Thus, by computing the integral over dη in the right hand side of the last equation of (103) and then, cancelling out the two same terms occurring in both sides of this equation, we find

v 3 0i (x) = v 3 0i (y) ⇔ x 1 y 1 e g 3 i V 2 (η, x 2 , c i )dη = x 2 y 2 e g 3 i V 1 (y 1 , ξ, c i )dξ. (104)
Hence, (100), ( 102) and (104) lead to the result announced in (36).

III. Proof of Lemma 5.2: Let ũi = ûi -u in Ω i × (0, T ). In view of (6), ũi solves

               L[ũ i ] = 0 in Ω i × (T 0 , T ), ũi (•, T 0 ) = ûi (•, T 0 ) -u(•, T 0 ) in Ω i , ũi = 0 on Γ i-1 ∪ Γ i × (T 0 , T ), D∇ũ i • ν = 0 on ∂Ω i \ (Γ i-1 ∪ Γ i ) × (T 0 , T ). ( 105
)
Multiplying the first equation of (105) by ũi and integrating by parts in Ω i using Green's formula and (3) gives: For all t ∈ (T 0 , T ),

1 2 d dt ũi (•, t) 2 L 2 (Ω i ) + R ũi (•, t) 2 L 2 (Ω i ) + Ω i D∇ũ i (•, t) -ũi (•, t)V • ∇ũ i (•, t) = 0. (106)
Since, from (3), div(V ) = 0 in Ω, it follows that div(V ũi ) = V • ∇ũ i in Ω. Afterwards, by applying an integration by parts and using Green's formula, we get

       Ω i ũi (•, t)V • ∇ũ i (•, t) = Ω i div V ũi (•, t) ũi (•, t) = - Ω i ũi (•, t)V • ∇ũ i (•, t) + ∂Ω i ũ2 i (•, t)V • ν, ∂Ω i ũ2 i (•, t)V • ν = Γ i-1 ∪Γ i ũ2 i (•, t)V • ν + ∂Ω i \(Γ i-1 ∪Γ i ) ũ2 i (•, t)V • ν = 0, (107)
=⇒

Ω i ũi (•, t)V • ∇ũ i (•, t) = 0, ∀ t ∈ (T 0 , T ). ( 108
)
Because the matrix D is uniformly elliptic in Ω and using Poincaré inequality, we obtain

     ∃γ > 0, γ ∇ũ i (•, t) 2 L 2 (Ω i ) ≤ Ω i D∇ũ i (•, t) • ∇ũ i (•, t), ∃C p > 0, ũi (•, t) L 2 (Ω i ) ≤ C p ∇ũ i (•, t) L 2 (Ω i ) , (109)
=⇒ γ C 2 p ũi (•, t) 2 L 2 (Ω i ) ≤ Ω i D∇ũ i (•, t) • ∇ũ i (•, t), ∀ t ∈ (T 0 , T ). (110)
Therefore, in view of ( 108)-( 110), it follows from (106) that

d dt ũi (•, t) 2 L 2 (Ω i ) ≤ -2 R + γ C 2 p ũi (•, t) 2 L 2 (Ω i ) , ∀t ∈ (T 0 , T ). (111)
Afterwards, by applying Gronwall's Lemma on (111), we get ũi (•, t)

2 L 2 (Ω i ) ≤ ũi (•, T 0 ) 2 L 2 (Ω i ) exp -2 R + γ C 2 p (t -T 0 ) , ∀t ∈ (T 0 , T ). (112)
That is the result announced in (56

) for α = 2 R + γ C 2 p .

B. Baby Example:

To establish the proof of Theorem 8.1, we prove the following two technical lemmas:

Lemma 9.1 Let ψ be the function involved in (86) and (n, m) ∈ {1, . . . , N } 2 . The following three functions defined in (0, ) 2 from ψ by

• ϕ n,m λS (x, y) = ψ(y) ψ ( ) ψ(x -S n )ψ ( -S m ) -ψ (x -S m )ψ ( -S n ) -ψ(x -S n )ψ (y -S m ), • ϕ n,m SS (x, y) = ψ(y) ψ ( ) ψ (x -S m )ψ ( -S n ) -ψ (x -S n )ψ ( -S m ) + ψ (x -S n )ψ (y -S m ), • ϕ n,m λλ (x, y) = ψ(y) ψ ( ) ψ(x -S m )ψ ( -S n ) -ψ(x -S n )ψ ( -S m ) + ψ(x -S n )ψ(y -S m ), are symmetric: ϕ n,m λS (x, y) = ϕ n,m λS (y, x), ∀(x, y) ∈ (0, ) 2 . The same holds for ϕ n,m SS , ϕ n,m λλ . Proof. Let ζ n,m i,j (x, y) = D 2 (r 1 -r 2 ) 2 ψ (i) (x -S n )ψ (j) (y -S m )
, where (n, m) ∈ {1, . . . , N } 2 , (i, j) ∈ IN 2 and ψ (k=i,j) is the k th derivative of the function ψ in (86). That leads to ζ n,m i,j (x, y) = r i 1 e r 1 (x-Sn) -r i 2 e r 2 (x-Sn) r j 1 e r 1 (y-Sm) -r j 2 e r 2 (y-Sm) .

Afterwards, from setting y = in (113), it follows that: For all i ∈ IN and j ∈ IN * ,

ζ n,m i,j (x, ) -ζ m,n j-1,i+1 (x, ) = β n,m i,j r 2 e r 1 x+r 2 -r 1 e r 2 x+r 1 , (114)
where β n,m i,j = r j-1 1 r i 2 e -r 1 Sm-r 2 Sn -r i 1 r j-1 2 e -r 1 Sn-r 2 Sm . Then, using (114), we get ψ(y) ψ ( ) ζ n,m i,j (x, ) -ζ m,n j-1,i+1 (x, ) = βn,m i,j

r 2 e r 1 (x+y)+r 2 + r 1 e r 2 (x+y)+r 1 -βn,m i,j

r 1 e r 1 y+r 2 x+r 1 + r 2 e r 1 x+r 2 y+r 2 , (

with βn,m i,j = β n,m i,j / r 1 e r 1 -r 2 e r 2 . Furthermore, from (113), we verify that ζ n,m i,j-1 (x, y) = r i 1 e r 1 (x-Sn) -r i 2 e r 2 (x-Sn) r j-1 1 e r 1 (y-Sm) -r j-1 2 e r 2 (y-Sm) . = r i+j-1 1 e r 1 (x+y-Sn-Sm) + r i+j-1 2 e r 2 (x+y-Sn-Sm)

-r i 1 r j-1 2 e r 1 x+r 2 y-r 1 Sn-r 2 Sm -r j-1 1 r i 2 e r 1 y+r 2 x-r 1 Sm-r 2 Sn . + r i+1 1 r j-1 2 e r 1 -r 1 Sn-r 2 Sm -r j-1 1 r i+1 2 e r 2 -r 1 Sm-r 2 Sn r 1 e r 1 -r 2 e r 2 e r 1 x+r 2 y + e r 1 y+r 2 x .

(117)

Besides, using ζ n,m i,j , the functions introduced in Lemma 9.1 can be rewritten as follows: * ϕ n,m λS (x, y) =

Hence, the symmetry of the functions ϕ n,m λS , ϕ n,m SS and ϕ n,m λλ in (118) is an immediate consequence of the symmetry with respect to x and y of the right hand side in (117). Lemma 9.2 Let (n, m) ∈ {1, . . . , N } 2 and (p i , P i ) ∈ (0, ) 2 . If p i and P i lie in a same side with respect to {S n , S m } i.e., H(p i -S n ) = H(p i -S m ) = H(P i -S n ) = H(P i -S m ), then we have

• ∂ λp u(p i )∂ Sq u(P i ) -∂ Sq u(p i )∂ λp u(P i ) = 0, ∀(p, q) ∈ {n, m} 2 , • ∂ Sn u(p i )∂ Sm u(P i ) -∂ Sm u(p i )∂ Sn u(P i ) = 0, • ∂ λn u(p i )∂ λm u(P i ) -∂ λm u(p i )∂ λn u(P i ) = 0, ( 119
)
where H is the Heaviside function and u is the solution of (85) given in (86).

Proof. Let (p i , P i ) ∈ (0, ) 2 be such that H(p i -S n ) = H(p i -S m ) and H(P i -S n ) = H(P i -S m ). Then, from (86), it follows that: For all (p, q) ∈ {n, m} 2 , * 1 λ q ∂ λp u(p i )∂ Sq u(P i ) -∂ Sq u(p i )∂ λp u(P i ) = ϕ p,q λS (p i , P i ) + ψ(p i -S p )ψ (P i -S q ) 1 -H(P i -S p ) H(p i -S p ) -ϕ p,q λS (P i , p i ) + ψ(P i -S p )ψ (p i -S q ) 1 -H(p i -S p ) H(P i -S p ),

where ϕ n,m λS , ϕ n,m SS and ϕ n,m λλ are the three symmetric functions introduced in Lemma 9.1. Therefore, if p i and P i are both upstream S n and S m i.e., H(p i -S n,m ) = H(P i -S n,m ) = 0 then, all terms in (120)-(122) vanish. Moreover, if p i and P i are rather both downstream S n and S m i.e., H(p i -S n,m ) = H(P i -S n,m ) = 1 then, the terms in (120)-(122) vanish too due the symmetry of the three functions ϕ n,m λS , ϕ n,m SS and ϕ n,m λλ .

Proof of Theorem 8.1:

1. For N = 1, the determinant of the 2 × 2 Jacobian matrix J M is det J M = ∂ λ 1 u(p 1 )∂ S 1 u(P 1 ) -∂ S 1 u(p 1 )∂ λ 1 u(P 1 ). ( 123)

From (123) and by applying (120) for i = p = q = 1, it follows due to the symmetry of the function ϕ 1,1 λS introduced in Lemma 9.1 that if H(p 1 -S 1 ) = H(P 1 -S 1 ) = 0 or H(p 1 -S 1 ) = H(P 1 -S 1 ) = 1 then, det(J M ) = 0. However, if H(p 1 -S 1 ) = H(P 1 -S 1 ) i.e., for example 0 < p 1 < S 1 < P 1 < then, from applying (120) with i = p = q = 1 and since r 1 > 0 whereas r 2 < 0, we obtain det J M = -λ 1 D 2 (r 1 -r 2 ) e r 1 p 1 -e r 2 p 1 r 1 e r 1 -r 2 e r 2 r 2 e r 1 P 1 +r 2 -r 1 e r 2 P 1 +r 1 e -S 1 (r 1 +r 2 ) = 0. (124) 2. For N = 2: From computing the determinant of the 4 × 4 Jacobian matrix J M by developing firstly with respect to its last row then, with respect to its third row, we find

det J M = ∂ λ 2 u
• (i) If Cardinality p i=1,2 , P i=1,2 ∩ S 1 , S 2 ≤ 1 and since 0 < p 1 < P 1 < p 2 < P 2 < then, there exists i ∈ {1, 2} such that p i and P i lie both in a same side with respect to {S 1 , S 2 }. Hence, by applying Lemma 9.2 for n = 1 and m = 2, it follows from (119) that all terms in (125) vanish which leads to det(J M ) = 0.

• (ii) From setting 0 < p 1 < S 1 < P 1 < p 2 < S 2 < P 2 < and using (86), it follows that the last term in (125) vanishes. Besides, by computing the other five terms, we obtain det(J M ) = λ 1 λ 2 ψ(p 1 ) ψ ( ) ψ (P 1 -S 1 )ψ(p 2 -S 1 ) -ψ(P 1 -S 1 )ψ (p 2 -S 1 )

ψ ( -S 2 )ψ (P 2 -S 2 ) -ψ ( -S 2 )ψ(P 2 -S 2 ) .

= λ 1 λ 2 ψ(p 1 ) D 2 ψ ( ) e -(r 1 +r 2 )(S 1 +S 2 ) r 1 e r 1 +r 2 P 2 -r 2 e r 2 +r 1 P 2 e r 2 P 1 +r 1 p 2 -e r 1 P 1 +r 2 p 2 .

(126)

Since r 1 > 0 and r 2 < 0, it follows from the last equation in (126) that det(J M ) = 0.

T 0 λ

 0 n (t)e Rt dt = 0 and ∃T 0 ∈ (0, T ) / λ n = 0 in (T 0 , T ), ∀n ∈ {1, . . . , N }.

Figure 1 :

 1 Figure 1: Behaviour of A 1 i : (a) Mean velocity vector (b) Unidirectional velocity field.

Figure 2 :Figure 3 :

 23 Figure 2: Flow of mean velocity: (c) Localized source in Ω 1 (d) Localized source in Ω 2 .

1 -r i+j- 1 1e r 1 2 e r 2

 11122 all (i, j) ∈ IN × IN * and (x, y) ∈ (0, ) 2 , it follows from (115)-(116) thatψ(y) ψ ( ) ζ n,m i,j (x, ) -ζ m,n j-1,i+1 (x, ) -ζ n,m i,j-1 (x, y) = βn,m i,jr 2 e r 1 (x+y)+r 2 + r 1 e r 2 (x+y)+r (x+y-Sn-Sm) + r i+j-1 (x+y-Sn-Sm)

1 D 2

 2 (r 1 -r 2 ) 2 ψ(y) ψ ( ) ζ n,m 0,2 (x, ) -ζ m,n 1,1 (x, ) -ζ n,m 0,1 (x, y) , * ϕ n,m SS (x, y) = -1 D 2 (r 1 -r 2 ) 2 ψ(y) ψ ( ) ζ n,m 1,2 (x, ) -ζ m,n 1,2 (x, ) -ζ n,m 1,1 (x, y) , * ϕ n,m λλ (x, y) = -1 D 2 (r 1 -r 2 ) 2 ψ(y) ψ ( )ζ n,m 0,1 (x, ) -ζ m,n 0,1 (x, ) -ζ n,m 0,0 (x, y) .

  there is no significant sources occurring in Ω i . Go to End do.• Otherwise, proceed as follows:1. Determine from (51)-(52) the total amount discharged in Ω i .2. Select m and n such that Theorem 3.5 yields injectivity of Ψ m,n i in Ω i . 3. Compute the coefficients in (54) and solve the last equation in (53): If NO solution in Ω i , according to Theorem 4.1, MULTIPLE sources occur in Ω i . If ONE solution, from Theorem 4.2, it's the position of the UNIQUE source in Ω i .

	End do

  * 1 λ n λ m ∂ Sn u(p i )∂ Sm u(P i ) -∂ Sm u(p i )∂ Sn u(P i ) = ϕ n,m SS (p i , P i ) -ψ (p i -S n )ψ (P i -S m ) 1 -H(P i -S n ) H(p i -S n ) -ϕ n,m SS (P i , p i ) -ψ (P i -S n )ψ (p i -S m ) 1 -H(p i -S n ) H(P i -S n ), (121) * ∂ λn u(p i )∂ λm u(P i ) -∂ λm u(p i )∂ λn u(P i ) = ϕ n,m λλ (p i , P i ) -ψ(p i -S n )ψ(P i -S m ) 1 -H(P i -S n ) H(p i -S n ) -ϕ n,m λλ (P i , p i ) -ψ(P i -S n )ψ(p i -S m ) 1 -H(p i -S n ) H(P i -S n ),

  (p 2 )∂ S 2 u(P 2 ) -∂ S 2 u(p 2 )∂ λ 2 u(P 2 ) ∂ λ 1 u(p 1 )∂ S 1 u(P 1 ) -∂ S 1 u(p 1 )∂ λ 1 u(P 1 ) -∂ S 1 u(p 2 )∂ S 2 u(P 2 ) -∂ S 2 u(p 2 )∂ S 1 u(P 2 ) ∂ λ 1 u(p 1 )∂ λ 2 u(P 1 ) -∂ λ 2 u(p 1 )∂ λ 1 u(P 1 ) -∂ λ 2 u(p 2 )∂ S 1 u(P 2 ) -∂ S 1 u(p 2 )∂ λ 2 u(P 2 ) ∂ λ 1 u(p 1 )∂ S 2 u(P 1 ) -∂ S 2 u(p 1 )∂ λ 1 u(P 1 ) -∂ λ 2 u(p 1 )∂ S 1 u(P 1 ) -∂ S 1 u(p 1 )∂ λ 2 u(P 1 ) ∂ λ 1 u(p 2 )∂ S 2 u(P 2 ) -∂ S 2 u(p 2 )∂ λ 1 u(P 2 ) -∂ S 1 u(p 1 )∂ S 2 u(P 1 ) -∂ S 2 u(p 1 )∂ S 1 u(P 1 ) ∂ λ 1 u(p 2 )∂ λ 2 u(P 2 ) -∂ λ 2 u(p 2 )∂ λ 1 u(P 2 ) + ∂ λ 2 u(p 1 )∂ S 2 u(P 1 ) -∂ S 2 u(p 1 )∂ λ 2 u(P 1 ) ∂ λ 1 u(p 2 )∂ S 1 u(P 2 ) -∂ S 1 u(p 2 )∂ λ 1 u(P 2 ) .

Acknowledgement: This work was supported by M2SINUM project co-financed by the European fund(ERDF,18P03390/18E01750/18P02733) and the Normandie Regional Council.