1 Introduction

The topic addressed in this article belongs to a class of inverse problems where the main objective is to identify the location and size of pollution sources. In this area of research, the rate of impurity emission from the pollution sources using the measurements of the concentration of the pollutants in some monitoring region is also an interesting issue to investigate. In general, the associated forward problems are governed by a parabolic advection-dispersion-reaction equation in one [15] or more spatial dimensions [3, 24, 34]. Such inverse problems are motivated from many physical phenomena and real-life events such as leakages in nuclear plants, refinery and chemical laboratories; and fire outbreaks in buildings and forests, for example. However, most of the practical applications reported in the literature concern the monitoring the water quality of rivers and groundwater [3, 15, 24] and the air quality in large cities [34].

In this paper, we reconstruct a set of pollution sources in a fluid medium by measuring the concentration of the pollutants within a subregion of the domain. To be more precise, we consider a two dimensional fluid domain where the impurities are getting injected with an unknown velocity through finitely many distributed unrecognized sources localized in a compactly supported region. The main idea is to measure the concentration of the pollutant away from their sources and use this information to completely characterize the sources and the intensity of the pollution, i.e. to determine the number of sources, their locations, their sizes and the velocity of impurity integration. It is not a difficult task to produce numerical results which show that if all these variables are unknown the uniqueness of the problem is lost and hence the well-posedness. This motivates us to consider two different inverse problems. In the first problem, we assume that the velocity field is given and reconstruct the topology of the pollution sources by recovering their locations and sizes. In the second one, the sizes of the pollution sources are assumed to be known and we determine their locations and the mean velocity of the respective leakages. In our simplified mathematical setting, we are dealing with an inverse reconstruction problem whose forward counterpart is governed by a steady-state convection-diffusion equation.

In this case, since the pollutant concentration is known in a small subregion of the fluid domain, the inverse problem can be written in the form of an over-determined boundary value problem. Following the classical well-known strategy, we rewrite our inverse problem as an optimization problem which consists in minimizing the misfit between the known measurement and the solution to the model problem with respect to the parameters related to the unknown sources of pollution. One can see [5, 6, 22] where the inverse problems associated with the reconstruction of geometrical subdomains are reformulated as topology optimization problems. A quite general approach to deal with such problems is based on the concept of topological derivatives [28], which can be seen as a particular case of the broader class of asymptotic methods. The bibliography is quite exhaustive, but for a general idea on the subject, readers may refer to [1, 7, 10,11,12].

The topological derivative is defined as the first term (correction) of the asymptotic expansion of a given shape functional with respect to a small parameter that measures the size of singular domain perturbations, such as holes, inclusions, defects, source-terms and cracks. This relatively new concept has applications in many different fields such as shape and topology optimization, inverse problems, image processing, multi-scale material design and mechanical modeling including damage and fracture evolution phenomena. The topological derivative was rigorously introduced by Sokołowski and Żochowski [31] and developed later by many authors [18, 27, 30, 33]. It can be seen as a particular case of the broader class of asymptotic methods fully developed in the book by Ammari and Kang [2], for instance. One can refer to the books by Maz’ya et al. [25] and Nazarov and Plamenevskij [26] for similar topics and discussions. In addition, one can also cite the paper by Sokołowski and Żochowski [32] where the domain decomposition method was proposed in the domain of topological derivatives. More recently, the concept of second-order topological derivative was introduced in [14], which has been successfully used for solving a class of inverse reconstruction problems. A variety of applications of the topological derivative method in topology optimization, inverse imaging problem and inverse reconstruction problems can be found in the complementary book by Novotny, Sokołowski and Żochowski [29]. Following the original ideas introduced in [16], in the current article, the second-order topological derivative is applied in the context of the proposed pollution sources reconstruction problem. The two main challenges in this article, in comparison to [16] are as follows: (i) the presence of the convective term in the governing forward equation, leading to a non-symmetric and non-coercive bilinear form; and (ii) the reconstruction of a vector quantity representing the mean velocity of a pollutant substance that escapes of a confined region characterizing a pollution source to be determined.

The outline of this paper is as follows. The inverse problem is reformulated as a topology optimization problem in Sect. 2. Such reformulation of the problem allows to solve it through a topological derivatives-based approach. Thus, we introduce some tools from the theory of topological derivatives in the beginning of the Sect. 3, namely, the functionals associated to the unperturbed and perturbed domains, the ansatz for a scalar field of interest and some auxiliary boundary value problems. In Sect. 4, we state the main result of this paper which consists of the asymptotic expansion of the shape functional. Estimates of the remainders, obtained in Sect. 4, are presented in Appendix A. Numerical results are presented in Sect. 5 where two reconstruction algorithms are devised from the asymptotic expansion of the shape functional. In addition, numerical experiments are conducted in order to test the capabilities of each of the proposed algorithms. Some concluding remarks are inferred in Sect. 6.

2 Problem Formulation

Let \(\varOmega \subset \mathbb {R}^2\) be an open and bounded domain with smooth boundary \(\partial \varOmega \). Measurements of the pollutant concentration field are collected in a subdomain \(\varOmega _o\subset \varOmega \). We assume that there may be an unknown number (\(N^*\in \mathbb {Z}^+\)) of isolated pollution sources \(\omega ^*_i\) within the domain \(\varOmega \). See Fig. 1a. Therefore, there is a set \(\omega ^*= \cup _{i=1}^{N^*} \omega ^*_i\), whose components \(\omega _i^*\) satisfy \(\overline{\omega _i^*}\cap \overline{\omega _j^*} = \emptyset \) for \(i\ne j\) and \(\overline{\omega _i^*}\cap \partial \varOmega = \emptyset \) for each \(i, j \in \lbrace 1, \ldots , N^*\rbrace \).

Fig. 1
figure 1

Domain \(\varOmega \) with (a) and without (b) the set of isolated pollution sources \(\omega ^*\)

We assume that the polluting substance is leaking in an incompressible flow with velocity \(V \in C(\varOmega ,\mathbb {R}^2)\), such that \(V \ne (0,0)\) at each \(\omega _i^*\subset \varOmega \), \(i=1,\ldots ,N^*\). For a given Dirichlet data g imposed on \(\partial \varOmega \), the resulting pollutant concentration z in \(\varOmega \) is observed in the subdomain \(\varOmega _o\). In this set up, the inverse problem consists in finding \(\chi _{\omega ^*}\) such that the pollutant concentration z satisfies the following boundary value problem

$$\begin{aligned} \left\{ \begin{array}{rcll} - \varDelta z + \chi _{\omega ^*} V \cdot \nabla z &{}=&{} 0 &{} \;\text {in}\;\varOmega ,\\ z &{}=&{} g &{} \;\text {on}\;\partial \varOmega ,\\ \end{array} \right. \end{aligned}$$
(1)

where the parameter \(\chi _{\omega ^*}\) is defined as

$$\begin{aligned} \chi _{\omega ^*} = \left\{ \begin{array}{cl} 0 &{} \;\text {in}\;\varOmega \setminus \overline{\omega ^*}, \\ 1 &{} \;\text {in}\;\omega ^*_i,\;i = 1,\ldots ,N^*\end{array} \right. \end{aligned}$$
(2)

and the velocity V is considered as a null divergence vector field, namely, \({{\,\mathrm{div}\,}}V = 0\) in \(\omega _i^*\), for \(i = 1,\ldots ,N^*\). Without loss of generality, we are considering only one boundary data g on \(\partial \varOmega \). The extension to several boundary data is trivial. See [16] for further detail. Here, function g is a purely synthetic data used to produce the observation of z in \(\varOmega _o\). The realistic problem of pollution source reconstruction is much more complicated and requires further investigation.

Now, for an initial guess \(\chi _\omega \) of \(\chi _{\omega ^*}\), we consider the pollutant concentration field u to be the solution to the boundary value problem

$$\begin{aligned} \left\{ \begin{array}{rclll} -\varDelta u + \chi _{\omega } V \cdot \nabla u &{}=&{} 0 &{} \text {in} &{} \varOmega , \\ u &{}=&{} g &{} \text {on} &{} \partial \varOmega , \end{array} \right. \end{aligned}$$
(3)

where

$$\begin{aligned} \chi _\omega = \left\{ \begin{array}{cl} 0 &{} \;\text {in}\;\varOmega \setminus {\overline{\omega }}, \\ 1 &{} \;\text {in}\;\omega _i,\;i = 1,\ldots ,N. \end{array} \right. \end{aligned}$$
(4)

The characteristic function \(\chi _{\omega ^*}\) is unknown and hence z but we assume that z can be measured in \(\varOmega _o\). Thus, if we look for an appropriate \(\chi _{\omega ^*}\), then we wish u to agree with z in \(\varOmega _o\), i.e., we want \(u = z_{|_{\varOmega _o}}\). In addition, in such a problem, we assume that the velocity field V is known and then we reconstruct the support of \(\omega ^*\) from partial measurements of the scalar field z taken on \(\varOmega _o\).

The inverse problem under investigation is written in the form of an ill-posed and over-determined boundary value problem. It is well known that such difficulty can be overcome by rewriting the inverse problem in the form of an optimization problem. In particular, if the velocity field V is a known data and we reconstruct the support of each \(\omega _i^*\), \(i=1,\ldots ,N^*\), then the inverse problem of finding \(\chi _{\omega ^*}\) can be rewritten as a topology optimization problem, namely

$$\begin{aligned} {\underset{\omega \subset \varOmega }{\text{ Minimize }} \;\; \mathcal {J}_\omega (u) = \int _{\varOmega _o} (u - z)^2,} \end{aligned}$$
(5)

where z and u are the solutions of the boundary value problems (1) and (3), respectively, corresponding to the Dirichlet data g.

3 Topological Derivatives-Based Approach

The inverse reconstruction problem of Sect. 2 has been written in the form of a topology optimization problem (5). A quite general approach for solution of such class of problems is based on the concept of topological derivative which requires the expansion of the shape functional \(\mathcal {J}_\omega (u)\) with respect to the parameters depend upon a set of small perturbations. Since the topological derivative does not depend on the initial guess of the unknown topology \(\omega ^*\), we start with the unperturbed domain by setting \(\omega =\emptyset \). See Fig. 1b. More precisely, we consider

$$\begin{aligned} {\mathcal {J}(u_0) = \int _{\varOmega _o} \left( u_0 - z \right) ^2,} \end{aligned}$$
(6)

where \(u_0\) be the solution of the unperturbed boundary value problem

$$\begin{aligned} \left\{ \begin{array}{rcll} -\varDelta u_0 &{}=&{} 0 &{} \;\text {in}\;\varOmega , \\ u_0 &{}=&{} g &{} \;\text {on}\;\partial \varOmega . \end{array} \right. \end{aligned}$$
(7)

In this article, we are considering the topology optimization problem (5) for the ball-shaped pollution sources and hence we define the topologically perturbed counter-part of (7) by introducing \(N \in \mathbb {Z}^+\) number of small circular pollution sources \(B_{\varepsilon _i} (x_i)\) with center at \(x_i\in \varOmega \) and radius \(\varepsilon _i\) for \(i=1,\ldots ,N\). The respective geometric set can be denoted as

$$\begin{aligned} B_\varepsilon (\xi ) = \bigcup _{i=1}^{N} B_{\varepsilon _i} (x_i), \end{aligned}$$
(8)

where \(\xi = (x_1,\ldots ,x_N)\) and \(\varepsilon = (\varepsilon _1,\ldots ,\varepsilon _N)\). Moreover, we assume that \(\overline{B_\varepsilon }\cap \partial \varOmega =\emptyset \), \(\overline{B_\varepsilon }\cap \varOmega _o = \emptyset \) and \(\overline{B_{\varepsilon _i}(x_i)}\cap \overline{B_{\varepsilon _j}(x_j)}=\emptyset \) for each \(i \ne j\) and \(i,j \in \lbrace 1,\ldots ,N \rbrace \). The shape functional associated with the topologically perturbed domain is written as

$$\begin{aligned} {\mathcal {J}(u_\varepsilon ) = \int _{\varOmega _o} \left( u_\varepsilon - z \right) ^2} \end{aligned}$$
(9)

with \(u_\varepsilon \) be the solution of the perturbed boundary value problem

$$\begin{aligned} \left\{ \begin{array}{rcll} -\varDelta u_\varepsilon + \chi _\varepsilon V \cdot \nabla u_\varepsilon &{}=&{} 0 &{} \;\text {in}\;\varOmega , \\ u_\varepsilon &{}=&{} g &{} \;\text {on}\;\partial \varOmega , \end{array} \right. \end{aligned}$$
(10)

where the parameter \(\chi _\varepsilon \) is defined as

$$\begin{aligned} \chi _\varepsilon = \left\{ \begin{array}{cl} 0 &{} \;\text {in}\;\varOmega \setminus \overline{B_\varepsilon (\xi )}, \\ 1 &{} \;\text {in}\;B_{\varepsilon }(\xi ). \end{array} \right. \end{aligned}$$
(11)

By construction, we consider that the velocity field V satisfies the additional condition \({{\,\mathrm{div}\,}}V = 0\) in \(B_\varepsilon (\xi )\). Since we are interested in reconstructing the mean velocity of the leakages, in order to simplify the analysis we assume that the velocity V is constant in the neighborhood of \(B_{\varepsilon _i}(x_i)\), so that \(V(x) \approx V_i:=V(x_i)\), with \(x \in B_{\varepsilon _i}(x_i) \), \(i=1, \ldots , N\).

As mentioned earlier, the topological derivatives measure the sensitivity of the shape functional with respect to the parameters \((\varepsilon ,\xi )\) depending upon a set of small circular pollution sources \(B_{\varepsilon }(\xi )\). Therefore, our idea is to obtain the number, location and radius of the components of the set \(B_{\varepsilon }(\xi )\) in order to produce the best approximation to the true set of pollution sources \(\omega ^*\) by using the concept of topological derivatives.

In this article, we are interested in expanding the shape functional \(\mathcal {J}(u_\varepsilon )\) defined in (9) in power of \(\varepsilon \). Therefore, we start by simplifying the difference between the perturbed shape functional \(\mathcal {J}(u_\varepsilon )\) and its unperturbed counter-part \(\mathcal {J}(u_0)\) defined in (9) and (6), respectively, as follows

$$\begin{aligned} \mathcal {J}(u_\varepsilon ) - \mathcal {J}(u_0) = \int _{\varOmega _o} \left[ 2(u_\varepsilon - u_0) (u_0 - z) + (u_\varepsilon - u_0)^2\right] . \end{aligned}$$
(12)

Since we are approximating the set of pollution sources to be reconstructed \(\omega ^*\) by a number of circular balls \(B_{\varepsilon }(\xi )\), the idea is to expand the perturbed shape functional \(\mathcal {J}(u_\varepsilon )\) with respect to the Lebesgue measure (volume) of the two-dimensional ball \(B_{\varepsilon _i}(x_i)\), namely, \(|B_{\varepsilon _i}(x_i)|= \pi \varepsilon _i^2\), \(i = 1, \ldots , N\). To simplify the notation, we introduce the vector

$$\begin{aligned} \alpha = (\alpha _1,\ldots ,\alpha _N)\quad \text { with }\quad \alpha _i = |B_{\varepsilon _i}(x_i)|, \end{aligned}$$
(13)

and the characteristic function

$$\begin{aligned} \chi _i = \left\{ \begin{array}{cl} 0 &{} \;\text {in}\;\varOmega \setminus \overline{B_{\varepsilon _i}(x_i)}, \\ 1 &{} \;\text {in}\;B_{\varepsilon _i}(x_i). \end{array} \right. \end{aligned}$$
(14)

Additionally, for \(i = 1,\ldots ,N\), let us define the linear operators \(f \mapsto \psi [f]\) satisfying

$$\begin{aligned} \left\{ \begin{array}{rcll} \varDelta \psi [f] &{}=&{} 0 &{} \text{ in } \varOmega , \\ \psi [f] &{}=&{} -f &{} \text{ on } \partial \varOmega , \end{array} \right. \end{aligned}$$
(15)

and \(f \mapsto \psi _i[f]\) such that

$$\begin{aligned} \left\{ \begin{array}{rcll} \varDelta \psi _i[f] &{}=&{} \alpha _i^{-1} f \chi _i &{} \text{ in } \varOmega , \\ \psi _i[f] &{}=&{} 0 &{} \text{ on } \partial \varOmega . \end{array} \right. \end{aligned}$$
(16)

For a positive real number \(R\gg \varepsilon \), we consider the ball \(B_R(x_i)\) such that \(B_{\varepsilon _i}(x_i) \subset \varOmega \subset B_R(x_i)\) with \(x_i\in \varOmega \) for \(i = 1,\ldots ,N\). Then, we define the linear operator \(f \mapsto \psi _i^R[f]\) as follows

$$\begin{aligned} \left\{ \begin{array}{rcll} \varDelta \psi _i^R[f] &{}=&{} \alpha _i^{-1} f \chi _i &{} \text{ in } B_R(x_i), \\ \psi _i^R[f] &{}=&{} 0 &{} \text{ on } \partial B_R(x_i). \end{array} \right. \end{aligned}$$
(17)

By adopting the Einstein summation convention and taking into account the notation introduced above, let us consider the following ansatz for the asymptotic expansion of \(u_\varepsilon \) with respect to the parameters describing the small circular pollution sources as presented in Sect. 2

$$\begin{aligned} u_\varepsilon = u_0 + \alpha _ih_i^\varepsilon + \alpha _i\alpha _jv_{ij}^\varepsilon + \alpha _i\alpha _j^2 w_{ij}^\varepsilon + {\tilde{u}}_\varepsilon , \end{aligned}$$
(18)

where

$$\begin{aligned} h_i^\varepsilon := \psi _i[V_i \cdot \nabla u_0],\quad v_{ij}^\varepsilon := \psi _i[V_i \cdot \nabla h_j^\varepsilon ] \quad \text { and }\quad w_{ij}^\varepsilon := \psi _i[V_i \cdot \nabla v_{ij}^\varepsilon ], \end{aligned}$$
(19)

for \(i,j = 1,\ldots ,N\). The last term of the expansion (18), namely \({\tilde{u}}_\varepsilon \), must satisfy the problem

$$\begin{aligned} \left\{ \begin{array}{rcll} -\varDelta {\tilde{u}}_\varepsilon + \chi _\varepsilon V \cdot \nabla {\tilde{u}}_\varepsilon &{}=&{} \varPhi _\varepsilon &{} \text{ in } \varOmega , \\ {\tilde{u}}_\varepsilon &{}=&{} 0 &{} \text{ on } \partial \varOmega , \end{array} \right. \end{aligned}$$
(20)

with

$$\begin{aligned} \varPhi _\varepsilon = \varPhi _{\varepsilon ,1} + \varPhi _{\varepsilon ,2} \end{aligned}$$
(21)

such that, for \(i,j,l = 1,\ldots ,N\), one has

$$\begin{aligned} \varPhi _{\varepsilon ,1} = - {\alpha _j\alpha _lV_i \cdot \nabla v_{jl}^\varepsilon \;{\chi _i}_{|_{l \ne j}}} \quad \text { and }\quad \varPhi _{\varepsilon ,2} = - {\alpha _j\alpha _l^2 V_i \cdot \nabla w_{jl}^\varepsilon \;\chi _i}. \end{aligned}$$
(22)

In order to simplify further analysis, we write \(h_i^\varepsilon \) as a sum of three functions \(p_i^\varepsilon \), \(q_i\) and \({\tilde{h}}_i^\varepsilon \) in the form

$$\begin{aligned} h_i^\varepsilon = (V_i \cdot \nabla u_0(x_i)) (p^\varepsilon _i + q_i) + {\tilde{h}}_i^\varepsilon . \end{aligned}$$
(23)

The function \(p_i^\varepsilon \) is defined as

$$\begin{aligned} p_i^\varepsilon:= & {} \psi _i^R[1] + \frac{1}{2\pi } \ln {R} \nonumber \\ {}= & {} \left\{ \begin{array}{ll} \displaystyle {\frac{1}{4\pi \varepsilon _i^2} \Vert x-x_i\Vert ^2 - \frac{1}{2\pi }\left( \frac{1}{2} - \ln \varepsilon _i\right) } &{} \text{ in } B_{\varepsilon _i}(x_i),\\ \displaystyle {\frac{1}{2\pi }\ln {\Vert x-x_i\Vert }} &{} \text{ in } B_R (x_i)\setminus \overline{B_{\varepsilon _i}(x_i)}. \end{array} \right. \end{aligned}$$
(24)

Building on the expression (24), we construct \(q_{i}\), \(i = 1,\ldots ,N\), in order to compensate for the discrepancies introduced by \(p_i^\varepsilon \) on \(\partial \varOmega \). Therefore, we define

$$\begin{aligned} {q_i := \frac{1}{2\pi }\psi [\ln \Vert x - x_i\Vert ]}. \end{aligned}$$
(25)

Additionally, the residual function \({\tilde{h}}_i^\varepsilon \) is given by

$$\begin{aligned} {{\tilde{h}}_i^\varepsilon := \psi _i[V_i \cdot (\nabla u_0 - \nabla u_0(x_i))].} \end{aligned}$$
(26)

Taking into account the decomposition (23) together with the expressions (24)–(26), we introduce an alternative representation for the function \(h_i^\varepsilon \) outside the ball \(B_{\varepsilon _i}(x_i)\), namely

$$\begin{aligned} {h_i^\varepsilon = (V_i \cdot \nabla u_0(x_i)) h_i + {\tilde{h}}_i^\varepsilon \quad \text { in } \varOmega \setminus \overline{B_{\varepsilon _i}(x_i)},} \end{aligned}$$
(27)

where

$$\begin{aligned} {h_i := p_i + q_i} \end{aligned}$$
(28)

with

$$\begin{aligned} {p_i := \frac{1}{2\pi }\ln \Vert x - x_i\Vert \quad \text { and }\quad q_i:= \psi [p_i],} \end{aligned}$$
(29)

that is, for \(i = 1,\ldots ,N\), \(q_i\) solves the homogeneous boundary value problem

$$\begin{aligned} {\left\{ \begin{array}{rcll} \varDelta q_i &{}=&{} 0 &{} \text{ in } \varOmega , \\ q_i &{}=&{} -p_i &{} \text{ on } \partial \varOmega . \end{array} \right. } \end{aligned}$$
(30)

Now, building on the expression for \({\tilde{h}}_i^\varepsilon \) given by (26), we first define the vector quantity \(W_i\) coming out from the Taylor expansion of the argument of \(\psi _i\) as

$$\begin{aligned} {W_i := (\nabla ^2 u_0(x_i))V_i } \end{aligned}$$
(31)

and then the function \({\tilde{h}}_i^\varepsilon \) can be split into the following sum

$$\begin{aligned} {{\tilde{h}}_i^\varepsilon = {\tilde{p}}_i^\varepsilon + \varepsilon _i^{2} {\tilde{q}}_i + \tilde{{\tilde{h}}}_i^\varepsilon ,} \end{aligned}$$
(32)

where \({\tilde{p}}_i^\varepsilon \) is defined as

$$\begin{aligned}&{\tilde{p}}_i^\varepsilon := \psi _i^R[W_i \cdot (x - x_i)] \nonumber \\&\quad =\left\{ \begin{array}{ll} \displaystyle {\frac{1}{8 \pi }\left( \frac{\Vert x-x_i\Vert ^2 - \varepsilon _i^2}{\varepsilon _i^2} - \frac{R^2 - \varepsilon _i^2}{R^2}\right) W_i \cdot (x- x_i)} &{} \text{ in } B_{\varepsilon _i}(x_i),\\ \displaystyle {\frac{\varepsilon _i^{2}}{8 \pi R^2}\frac{\Vert x-x_i\Vert ^2 - R^2}{\Vert x-x_i\Vert ^2}\; W_i \cdot (x - x_i)} &{} \text{ in } B_R (x_i)\setminus \overline{B_{\varepsilon _i}(x_i)}, \end{array} \right. \nonumber \\ \end{aligned}$$
(33)

the function \({\tilde{q}}_i\) is such that

$$\begin{aligned} {{\tilde{q}}_i := \frac{1}{8 \pi R^2} \psi \left[ \frac{\Vert x-x_i\Vert ^2 - R^2}{\Vert x-x_i\Vert ^2}\; W_i \cdot (x - x_i)\right] } \end{aligned}$$
(34)

and the residual term \(\tilde{{\tilde{h}}}_i^\varepsilon \) is given by

$$\begin{aligned} {\tilde{{\tilde{h}}}_{i}^\varepsilon := \psi _i[T_i(x)],} \end{aligned}$$
(35)

where \(T_i(x) = O(\Vert x-x_i\Vert ^2)\) in an appropriate norm.

Similarly to the decomposition proposed for the function \(h_i^\varepsilon \) in (23), we write the function \(v_{ii}^\varepsilon \) as a sum of the functions \(p_{ii}^\varepsilon \), \(q_{ii}\) and \({\tilde{v}}_{ii}^\varepsilon \) in the form

$$\begin{aligned} {v_{ii}^\varepsilon = (V_i \cdot \nabla u_0(x_i)) (p^\varepsilon _{ii} + q_{ii}) + {\tilde{v}}_{ii}^\varepsilon .} \end{aligned}$$
(36)

The function \(p_{ii}^\varepsilon \) is such that

$$\begin{aligned}&p_{ii}^\varepsilon := \psi _i^R[V_i \cdot \nabla p_i^\varepsilon ]\nonumber \\&\quad =\left\{ \begin{array}{ll} \displaystyle {\frac{1}{(4\pi \varepsilon _i)^2}\left( \frac{\Vert x-x_i\Vert ^2 - \varepsilon _i^2}{\varepsilon _i^2} - \frac{R^2 - \varepsilon _i^2}{R^2}\right) V_i \cdot (x- x_i)} &{} \text{ in } B_{\varepsilon _i}(x_i),\\ \displaystyle {\frac{1}{(4\pi R)^2}\frac{\Vert x-x_i\Vert ^2 - R^2}{\Vert x-x_i\Vert ^2}\;V_i \cdot (x - x_i)} &{} \text{ in } B_R (x_i)\setminus \overline{B_{\varepsilon _i}(x_i)}, \end{array} \right. \nonumber \\ \end{aligned}$$
(37)

where we have used expression (24) to obtain

$$\begin{aligned} {\nabla p_i^\varepsilon (x) = \frac{x - x_i}{2\pi \varepsilon _i^2} \quad \text { in } B_{\varepsilon _i}(x_i).} \end{aligned}$$
(38)

Additionally, the functions \(q_{ii}\) and \({\tilde{h}}_i^\varepsilon \) are given by

$$\begin{aligned} {q_{ii} := \frac{1}{(4\pi R)^2} \psi \left[ \frac{\Vert x-x_i\Vert ^2 - R^2}{\Vert x-x_i\Vert ^2}\;V_i \cdot (x - x_i)\right] } \end{aligned}$$
(39)

and

$$\begin{aligned} {{\tilde{v}}_{ii}^\varepsilon := \psi _i[(V_i \cdot \nabla u_0(x_i)) (V_i \cdot \nabla q_i) + V_i \cdot \nabla {\tilde{h}}_i^\varepsilon ],} \end{aligned}$$
(40)

respectively.

Finally, in order to simplify further calculations related to the asymptotic development of the topologically perturbed shape functional (9), we introduce an adjoint state v to be the solution to the problem

$$\begin{aligned} \left\{ \begin{array}{rcll} - \varDelta v &{}=&{} (u_0 - z) \chi _{\varOmega _o} &{} \text{ in } \varOmega , \\ v &{}=&{} 0 &{} \text{ on } \partial \varOmega . \end{array} \right. \end{aligned}$$
(41)

4 Asymptotic Expansion of the Shape Functional

Let us start by replacing (18) in (12) to obtain

$$\begin{aligned}&\mathcal {J}(u_\varepsilon ) - \mathcal {J}(u_0) = 2 \alpha _i\int _{\varOmega _o} h_i^\varepsilon (u_0 - z) + 2 \alpha _i\alpha _j\int _{\varOmega _o} v_{ij}^\varepsilon (u_0 - z) \nonumber \\&\quad + 2 \alpha _i\alpha _j^2 \int _{\varOmega _o} w_{ij}^\varepsilon (u_0 - z) + \alpha _i\alpha _j\int _{\varOmega _o} h_i^\varepsilon h_j^\varepsilon + \sum _{\ell =1}^{9} \mathcal {E}_\ell (\varepsilon ), \end{aligned}$$
(42)

where

$$\begin{aligned} \mathcal {E}_1(\varepsilon )= & {} \int _{\varOmega _o} ( 2 (u_0 - z) + {\tilde{u}}_\varepsilon ){\tilde{u}}_\varepsilon ,\quad \mathcal {E}_2(\varepsilon ) = 2 \alpha _i\alpha _j\alpha _l\int _{\varOmega _o} h_i^\varepsilon v_{jl}^\varepsilon ,\nonumber \\&\quad \mathcal {E}_3(\varepsilon ) = 2 \alpha _i\alpha _j\alpha _l^2 \int _{\varOmega _o} h_i^\varepsilon w_{jl}^\varepsilon , \end{aligned}$$
(43)
$$\begin{aligned} \mathcal {E}_4(\varepsilon )= & {} 2 \alpha _i\int _{\varOmega _o} h_i^\varepsilon {\tilde{u}}_\varepsilon ,\quad \mathcal {E}_5(\varepsilon ) = \alpha _i\alpha _j\alpha _l\alpha _p\int _{\varOmega _o} v_{ij}^\varepsilon v_{lp}^\varepsilon ,\nonumber \\&\quad \mathcal {E}_6(\varepsilon ) = 2 \alpha _i\alpha _j\alpha _l\alpha _p^2 \int _{\varOmega _o} v_{ij}^\varepsilon w_{lp}^\varepsilon , \end{aligned}$$
(44)
$$\begin{aligned} \mathcal {E}_7(\varepsilon )= & {} 2 \alpha _i\alpha _j\int _{\varOmega _o} v_{ij}^\varepsilon {\tilde{u}}_\varepsilon ,\quad \mathcal {E}_8(\varepsilon ) = \alpha _i\alpha _j^2\alpha _l\alpha _p^2 \int _{\varOmega _o} w_{ij}^\varepsilon w_{lp}^\varepsilon \quad \text { and }\nonumber \\&\quad \mathcal {E}_9(\varepsilon ) = 2 \alpha _i\alpha _j^2 \int _{\varOmega _o} w_{ij}^\varepsilon {\tilde{u}}_\varepsilon . \end{aligned}$$
(45)

On the basis of the boundary value problem for the adjoint state (41) together with the Green’s formula, we obtain

$$\begin{aligned} \int _{\varOmega _o} \psi _i[f] (u_0 - z) = - \int _{\varOmega } \psi _i[f] \varDelta v = - \frac{1}{\alpha _i} \int _{B_{\varepsilon _i} (x_i)} f v. \end{aligned}$$
(46)

The closed formula given by expression (46) allows to replace the first three integrals over \(\varOmega _o\) on the right-hand side of (42) by integrals over the ball \(B_{\varepsilon _i}(x_i)\). In fact, by taking into account the notation introduced in (19) and choosing f in the form of the scalar products \(V_i \cdot \nabla u_0\), \(V_i \cdot \nabla h_j^\varepsilon \) and \(V_i \cdot \nabla v_{ij}^\varepsilon \) in (46), expression (42) can be rewritten as

$$\begin{aligned} \mathcal {J}(u_\varepsilon ) - \mathcal {J}(u_0)= & {} - 2 \int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla u_0) v - 2 \alpha _j\int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla h_j^\varepsilon ) v \nonumber \\&- 2 \alpha _j^2 \int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla v_{jj}^\varepsilon ) v + \alpha _i\alpha _j\int _{\varOmega _o} h_i^\varepsilon h_j^\varepsilon + \sum _{\ell =1}^{9} \mathcal {E}_\ell (\varepsilon ),\nonumber \\ \end{aligned}$$
(47)

where the Einstein summation convention still holds. For the second and third integrals on the right-hand side of the expansion (47), we first split the sum over j in the cases when \(i = j\) and \(i \ne j\), then we replace \(h_j^\varepsilon \) and \(v_{jj}^\varepsilon \) by the decompositions (23) and (36), respectively. Since the fourth integral in (47) is evaluated over \(\varOmega _o\), then the functions \(h_i^\varepsilon \) (and also \(h_j^\varepsilon \)) are replaced by expression given in (27). Thus, the asymptotic expansion (47) can be expressed as

$$\begin{aligned}&\mathcal {J}(u_\varepsilon ) - \mathcal {J}(u_0) = - 2 \int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla u_0) v- 2 \alpha _i\int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla {\tilde{p}}_i^\varepsilon ) v \nonumber \\&\quad - 2 \alpha _i(V_i\cdot \nabla u_0(x_i)) \int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla p_i^\varepsilon ) v - 2 \alpha _i(V_i\cdot \nabla u_0(x_i)) \int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla q_i) v \nonumber \\&\quad - 2 \alpha _j(V_j\cdot \nabla u_0(x_j)) \int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla h_j) {v}_{|_{i \ne j}} - 2 \alpha _i^2 (V_i\cdot \nabla u_0(x_i)) \int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla p_{ii}^\varepsilon ) v \nonumber \\&\quad + \alpha _i\alpha _j(V_i\cdot \nabla u_0(x_i)) (V_j\cdot \nabla u_0(x_j)) \int _{\varOmega _o} h_i h_j + \sum _{\ell =1}^{15} \mathcal {E}_\ell (\varepsilon ). \end{aligned}$$
(48)

Here, the six new remainders are defined as

$$\begin{aligned} \mathcal {E}_{10}(\varepsilon )= & {} - 2 \alpha _i\int _{B_{\varepsilon _i}(x_i)} V_i \cdot (\varepsilon _i^{2} \nabla q_i + \nabla \tilde{{\tilde{h}}}_i^\varepsilon ) v,\nonumber \\&\quad \mathcal {E}_{11}(\varepsilon ) = - 2 \alpha _j\int _{B_{\varepsilon _i}(x_i)} (V_i \cdot \nabla {\tilde{h}}_j^\varepsilon ) {v}_{|_{i \ne j}}, \end{aligned}$$
(49)
$$\begin{aligned} \mathcal {E}_{12}(\varepsilon )= & {} - 2 \alpha _j^2 (V_j\cdot \nabla u_0(x_j)) \int _{B_{\varepsilon _i} (x_i)} (V_i \cdot \nabla q_{jj}) v,\nonumber \\&\quad \mathcal {E}_{13}(\varepsilon ) = - 2 \alpha _j^2 \int _{B_{\varepsilon _i} (x_i)} (V_i \cdot \nabla {\tilde{v}}_{jj}^\varepsilon )v, \end{aligned}$$
(50)
$$\begin{aligned} \mathcal {E}_{14}(\varepsilon )= & {} - 2 \alpha _j^2 (V_j\cdot \nabla u_0(x_j)) \int _{B_{\varepsilon _i} (x_i)} (V_i \cdot \nabla p_{jj}^\varepsilon ) {v}_{|_{i \ne j}} \end{aligned}$$
(51)

and

$$\begin{aligned} \mathcal {E}_{15}(\varepsilon ) = \alpha _i\alpha _j\int _{\varOmega _o} [(V_i\cdot \nabla u_0(x_i)) h_i{\tilde{h}}_j^\varepsilon + (V_j\cdot \nabla u_0(x_j)) h_j{\tilde{h}}_i^\varepsilon + {\tilde{h}}_i^\varepsilon {\tilde{h}}_j^\varepsilon ]. \end{aligned}$$
(52)

In order to obtain the final form of the asymptotic expansion of the shape functional \(\mathcal {J}(u_\varepsilon )\) from (48) we use the explicit representation of the functions \(p_i^\varepsilon \), \({\tilde{p}}_i^\varepsilon \) and \(p_{ii}^\varepsilon \), given by (24), (33) and (37), respectively, to compute their gradients in \(B_{\varepsilon _i}(x_i)\). Moreover, we also consider the Taylor’s expansions of the functions \(\nabla u_0\), \(\nabla q_i\), \(\nabla h_j\) and v around the point \(x_i\) up to produce a remainder of order \(o(|\varepsilon |^4)\) as done in [16]. Thus, we get

$$\begin{aligned}&\mathcal {J}(u_\varepsilon ) - \mathcal {J}(u_0) = - 2 \alpha _i(V_i\cdot \nabla u_0(x_i)) v(x_i) - \frac{1}{2\pi } \alpha _i^2\;V_i\cdot \nabla ^2 u_0(x_i) \nabla v(x_i)\nonumber \\&\quad +\, \frac{1}{4\pi } \alpha _i^2\;V_i \cdot (\nabla ^2 u_0(x_i)V_i) v(x_i)- \frac{1}{4\pi } \alpha _i^2 (V_i\cdot \nabla u_0(x_i))(V_i\cdot \nabla v(x_i)) \nonumber \\&\quad -\, 2 \alpha _i^2 (V_i\cdot \nabla u_0(x_i))(V_i \cdot \nabla q_i(x_i))v(x_i) +\, \frac{1}{8\pi } \alpha _i^2 (V_i\cdot \nabla u_0(x_i)) \Vert V_i\Vert ^{2} v(x_i)\nonumber \\&\quad - 2 \alpha _i\alpha _j(V_j\cdot \nabla u_0(x_j))(V_i\cdot \nabla h_j(x_i)){v(x_i)}_{|_{i \ne j}} \nonumber \\&\quad +\, \alpha _i\alpha _j(V_i\cdot \nabla u_0(x_i)) (V_j\cdot \nabla u_0(x_j)) \int _{\varOmega _o} h_i h_j + \mathcal {E}(\varepsilon ), \end{aligned}$$
(53)

provided that \(u_0\) and v are harmonic in \(\varOmega \setminus \varOmega _o\), with

$$\begin{aligned} \mathcal {E}(\varepsilon ) = \mathcal {E}_0 (\varepsilon ) + \sum _{\ell =1}^{15} \mathcal {E}_\ell (\varepsilon ). \end{aligned}$$
(54)

The new remainder \(\mathcal {E}_0(\varepsilon )\) contains the residual terms arising from different Taylor’s expansions arguing in a similar fashion to get the remainders \(\mathcal {E}_{10}(\varepsilon )\)\(\mathcal {E}_{19}(\varepsilon )\) in [16].

We are now in position to state the main result of this paper.

Theorem 1

Let \(u_0\), \(h_i\), \(q_i\) and v be the functions defined in (7), (28), (30) and (41), respectively, for \(i=1,\ldots ,N\). Additionally, let V be the null divergence velocity field introduced in problem (1). Then, the asymptotic expansion for the topologically perturbed shape functional \(\mathcal {J}(u_\varepsilon )\) defined in (9), with respect to \(\alpha _i = |B_{\varepsilon _i}(x_i)|\), \(i=1,\ldots ,N\), is given by (53), where the remaining terms, given by (54), are such that \(|\mathcal {E}(\varepsilon )| = o(|\varepsilon |^4)\), whose justification can be found in Appendix A.

5 Numerical Experiments

This section consists of two reconstruction algorithms and their respective numerical examples. Firstly, we propose the following contextualization of the inverse problem in order to recall some of its characteristic elements. Let us consider a pipeline network within a reference domain \(\varOmega \) which contains a \(N^*\) number of leakage points through which a polluting substance escapes with a velocity V. Each leakage point \(\omega ^*_i\) (\(i = 1,\ldots ,N^*\)) can be imagined as a pollution source which has to be reconstructed by a ball-shaped geometrical subdomain \(B_{\varepsilon _i}(x_i)\). Hence, it can be completely characterized by the center \(x_i\) and radius \(\varepsilon _i\) of \(B_{\varepsilon _i}(x_i)\), for \(i = 1,\ldots ,N^*\). Moreover, the vector field V at the point \(x_i\) provides us an additional information about each leakage.

Since the asymptotic expansion (53) depends on (a) the number, (b) the locations, (c) the sizes of the pollution sources and (d) the velocity of the pollution, we solve two different problems keeping in mind the physical aspects of the current scenario. In both problems to be considered, the number of pollution sources \(N^*\) is assumed to be given. In case it is not known, see Sect. 6.

In the first case, we assume that the velocity field V is given and we reconstruct the topology of the pollution sources by recovering their locations \(\xi \) and sizes \(\alpha \). We deal with this problem in Sect. 5.1 under the nomenclature support velocity reconstruction. In Sect. 5.2, we consider the second case with the name mean velocity reconstruction. Here, we assume that the sizes \(\alpha \) of the pollution sources are known and determine their locations \(\xi \) and respective velocities \(V_i\) of the leakages, \(i=1,\ldots ,N\).

Numerical experiments are conducted in the following scenario. The reference geometrical domain is given by a square \(\varOmega = (\)-\(\,0.5,0.5) \times (\)-\(\,0.5,0.5)\) which is discretized with three-node finite elements. The mesh is generated from a grid of \(160 \times 160\) squares. Each square is divided into four triangles which leads us to a resulting mesh comprising 102400 elements and 51521 nodes. Measurements of scalar fields of interest are taken in the subdomain \(\varOmega _o\) which is a union of four small regions in the neighborhood of the corners of the square domain \(\varOmega \). See Fig. 2a. The boundary \(\partial \varOmega \) of the reference domain is excited by imposing different Dirichlet data, namely, \(g_1 = x\), \(g_2 = y\) and \(g_3 = xy\).

The auxiliary boundary value problems are solved over the resulting mesh. From these solutions the topological derivatives can be numerically evaluated at any point of the mesh. However, due the high complexity of the algorithm [23], a sub-grid X consists of uniformly distributed points is extracted from the original mesh and then used as a set of admissible locations where a combinatorial search is performed which leads to the optimal solution defined in X. In the case of the examples presented below, we consider a fixed sub-grid X comprising 165 points, as illustrated in Fig. 2b.

Fig. 2
figure 2

Numerical scenario: a the set \(\varOmega _o\) (in gray color) as part of the reference domain \(\varOmega \) and b the set of admissible locations X

In order to obtain noisy synthetic data, the true target function z, corresponding to the solution of the problem (1) for a given Dirichlet data g, is corrupted with white Gaussian noise, where the resulting level of noise in the measurement of z is given by

$$\begin{aligned} \mu ^*= \frac{\Vert z_{\mu } - z\Vert _{L^2(\varOmega _o)}}{\Vert z\Vert _{L^2(\varOmega _o)}}\times 100 \%, \end{aligned}$$
(55)

where \(z_{\mu }\) denotes the corrupted target function used as synthetic data.

In all our numerical results, we represent the (true and reconstructed) pollution sources by black, the subdomain \(\varOmega _o\) by gray and the remaining domain by white colors.

5.1 Support Velocity Reconstruction

In this section, we present our first algorithm followed by two numerical examples concerning the reconstruction of the topology of the pollution sources. The solution algorithm is devised from the topological asymptotic expansion of the shape functional \(\mathcal {J}(u_\varepsilon )\) given by (53). After disregarding all terms of order \(o(|\alpha |^2)\) in (53), the resulting truncated expansion, in its matrix form, is written as

$$\begin{aligned} \delta J(\alpha ,\xi ,N):= \alpha \cdot d(\xi ) + \frac{1}{2} H(\xi ) \alpha \cdot \alpha , \end{aligned}$$
(56)

where the entries of the vector \(d\in \mathbb {R}^{N}\) and the matrix \(H \in \mathbb {R}^{N \times N}\) are given by

$$\begin{aligned} d_i := - 2 (V_i\cdot \nabla u_0(x_i)) v(x_i), \end{aligned}$$
(57)

and

$$\begin{aligned} H_{ii}:= & {} - \frac{1}{\pi } V_i \cdot \nabla ^2 u_0(x_i) \nabla v(x_i) + \frac{1}{2\pi } V_i \cdot (\nabla ^2 u_0(x_i)\;V_i) v(x_i)\nonumber \\&-\, \frac{1}{2\pi } (V_i \cdot \nabla u_0(x_i)) (V_i \cdot \nabla v(x_i))\nonumber \\&-\, 4 (V_i \cdot \nabla u_0(x_i)) (V_i \cdot \nabla q_i(x_i)) v(x_i) + \frac{1}{4\pi } (V_i \cdot \nabla u_0(x_i)) \Vert V_i\Vert ^2 v(x_i) \nonumber \\&+\, 2 (V_i \cdot \nabla u_0(x_i))^2 \int _{\varOmega _o} h_i^2, \end{aligned}$$
(58)
$$\begin{aligned} H_{ij}:= & {} - 4 (V_j \cdot \nabla u_0(x_j)) (V_i \cdot \nabla h_j(x_i)) v(x_i) \nonumber \\&+\, 2 (V_i \cdot \nabla u_0(x_i)) (V_j \cdot \nabla u_0(x_j)) \int _{\varOmega _o} h_i h_j, \end{aligned}$$
(59)

if \(i\ne j\), respectively, for \(i,j = 1,\ldots ,N\). Here there is no summation in the repeated indexes.

The truncated asymptotic expansion \(\delta J(\alpha ,\xi ,N)\), given by (56), is a quadratic form in \(\alpha \). Hence, in order to minimize (56), we proceed with the derivative of \(\delta J(\alpha ,\xi ,N)\) with respect to the variable \(\alpha \) to write the first-order optimality condition

$$\begin{aligned} \langle D_{\alpha } \delta J, \beta \rangle = [d(\xi ) + H^{s}(\xi ) \alpha ] \cdot \beta = 0, \quad \forall \beta \in \mathbb {R}^N, \end{aligned}$$
(60)

with \(H^{s} = (H + H^{\top })/2\), where \(H^\top \) denotes the transpose of the matrix H. From (60), we obtain the linear system

$$\begin{aligned} \alpha = - (H^{s}(\xi ))^{-1}d(\xi ). \end{aligned}$$
(61)

Note that, the quantity \(\alpha \), solution of (61), becomes a function of the locations \(\xi \), i.e., \(\alpha = \alpha (\xi )\). Let us replace \(\alpha = \alpha (\xi )\) in (56) to obtain

$$\begin{aligned} \delta J(\alpha (\xi ),\xi ,N) = \frac{1}{2} d(\xi ) \cdot \alpha (\xi ). \end{aligned}$$
(62)

Therefore, the optimal locations and sizes of the pollution sources, given by the pair of vectors (\(\xi ^\star \), \(\alpha ^\star \)), can be obtained by

$$\begin{aligned} \xi ^{\star } := \underset{\xi \in X}{\text {argmin}}\;\;\delta J(\alpha (\xi ),\xi ,N) \quad \text {and} \quad \alpha ^{\star } := \alpha (\xi ^{\star }), \end{aligned}$$
(63)

where X denotes the set of admissible locations for the pollution sources. The above procedure written in pseudo-code format can be found in [23].

5.1.1 Numerical Examples

Two numerical examples are presented below. The first one tests the robustness of the reconstruction in the presence of noisy data while the second one analyses the sensitivity of the reconstruction with respect to the number of measurements.

Example 1

reconstruction of a single leakage from noisy data. We consider a target domain containing a single leakage \(\omega ^*\) located at \(x^*= (0.2,0.1)\) whose radius is \(\varepsilon ^*= 0.03\), as shown in Fig. 3a. In addition, it is assumed that the pollution flows with constant velocity \(V = (2,2)\). The reconstruction of the topology of \(\omega ^*\) is performed from only one measurement by choosing \(g_1\) as Dirichlet data. In the absence of noise, the proposed algorithm is able to find the exact center \(x^*\) and size \(\varepsilon ^*\) of \(\omega ^*\). In addition, the pollution source is accurately reconstructed for levels of noise up to \(0.1\%\). On the other hand, the reconstruction degrades for a level of noise around \(0.3\%\). These results are reported in Fig. 3b–d.

Fig. 3
figure 3

Example 1: Target (a) and results (bd)

Example 2

reconstruction of multiple leakages. In this example, the pipeline network within the reference domain contains four leakages which are located at the points \(x_1^*= (0,0.3)\), \(x_2^*= (0.3,0)\), \(x_3^*= (0,-0.3)\) and \(x_4^*= (-0.3,0)\). The sizes of the leakages are equal, namely, \(\varepsilon _i^*= 0.03\), for \(i = 1,\ldots ,4\). We illustrate the target domain in Fig. 4a. The pollution substance escapes from the pipeline network with a velocity given by \(V = (1,10x^2)\). We start by considering only one measurement associated with \(g_1\). After comparing Fig. 4a, b, we observe that the algorithm fails in reconstructing the target domain. For a large number of pollution sources, the information obtained from only one measurement is not sufficient to accurately reconstruct them. Therefore, we increase the number of measurements by taking \(g_1\) and \(g_2\), whose result is shown in Fig. 4c. Since the reconstruction still fails, we add one more excitation, namely \(g_1\), \(g_2\) and \(g_3\). As we can observe in Fig. 4d, the pollution sources were successfully reconstructed after considering three measurements. In this case, we have obtained the exact location of the leakages while their sizes are approximately equal to the true values, namely, \(\varepsilon _i^\star = 0.0299\), for \(i = 1,\ldots ,4\).

Fig. 4
figure 4

Example 2: Target (a) and results (bd)

5.2 Mean Velocity Reconstruction

Here, we present our second algorithm concerning the mean velocity reconstruction for leakages whose sizes are known. We demonstrate the effectiveness of the algorithm with a numerical example at the end of this section. Analogous to the mathematical steps taken in Section 5.1, we start by considering the topological asymptotic expansion of the shape functional \(\mathcal {J}(u_\varepsilon )\) given by (53). Firstly, all the terms of order \(o(|\alpha |^2)\) together with the sixth term on the right-hand side of (53) are disregarded. This procedure allows us to obtain a linear system with respect to the variable of interest similar to (61). Next, we take the resulting truncated expansion and rewrite the velocity field in the form

$$\begin{aligned} V_i = \alpha _i \frac{V_i}{\alpha _i} = \alpha _i {\overline{V}}_i \end{aligned}$$
(64)

in order to work with the two-dimensional flow rate of the leakages at the points \(x_i\), \(i = 1,\ldots ,N\). Taking into account these two steps, we obtain the truncated expansion

$$\begin{aligned} \delta J({\overline{V}},\xi ,N) := {\overline{V}} \cdot d(\xi ) + \frac{1}{2} H(\xi ) {\overline{V}} \cdot {\overline{V}}, \end{aligned}$$
(65)

where \( d,{\overline{V}}\in \mathbb {R}^{2N}\) and \( H\in \mathbb {R}^{2N \times 2N}\) can be seen as the first-order topological derivative, the mean flow of the leakages and the higher-order topological derivative, respectively. More precisely, for a fixed number N of leakages to be reconstructed, the unknown flow vector \({\overline{V}}\in \mathbb {R}^{2N}\) is given by \({\overline{V}} := ({\overline{V}}_1, \ldots , {\overline{V}}_N)^\top \) where each entry \({\overline{V}}_i\) is a vector in \(\mathbb {R}^2\) written as \({\overline{V}}_i := ({\overline{V}}_i^1,{\overline{V}}_i^2)^\top \), for \(i = 1,\ldots ,N\). Thus, the components of \({\overline{V}}_i\) must be interpreted as follows: \({\overline{V}}_i^1\) and \({\overline{V}}_i^2\) are the flow rate of the polluting substance at the point \(x_i\) in the x- and y-directions, respectively. Similarly, the vector \( d\in \mathbb {R}^{2N}\) is written as \(d := (d_1, \ldots , d_N)^\top \) where each entry \(d_i\) is a vector in \(\mathbb {R}^2\) given by

$$\begin{aligned} d_i := - 2 \alpha _i^2 \nabla u_0(x_i) v(x_i) - \frac{1}{2\pi } \alpha _i^3 \nabla ^2 u_0(x_i) \nabla v(x_i), \end{aligned}$$
(66)

for \(i = 1,\ldots ,N\). Each entry \( H_{ij} \in \mathbb {R}^{2 \times 2}\) of the matrix \( H\in \mathbb {R}^{2N \times 2N}\)

$$\begin{aligned} H := \begin{pmatrix} H_{11} &{} \cdots &{} H_{1N} \\ \vdots &{} \ddots &{} \vdots \\ H_{N1} &{} \cdots &{} H_{NN} \end{pmatrix} \end{aligned}$$
(67)

is given by

$$\begin{aligned} H_{ii}:= & {} \frac{1}{2\pi } \alpha _i^4 \nabla ^2 u_0(x_i) v(x_i) - \frac{1}{2\pi } \alpha _i^4 \nabla u_0(x_i) \otimes \nabla v(x_i)\nonumber \\&-\, 4 \alpha _i^4 (\nabla u_0(x_i) \otimes \nabla q_i(x_i)) v(x_i) + 2 \alpha _i^4 \nabla u_0(x_i) \otimes \nabla u_0(x_i) \int _{\varOmega _o} h_i^2\quad \quad \end{aligned}$$
(68)

and

$$\begin{aligned} H_{ij} := - 4 \alpha _i^2 \alpha _j^2 (\nabla u_0(x_i) \otimes \nabla h_i(x_j)) v(x_j) + 2 \alpha _i^2 \alpha _j^2 \nabla u_0(x_i) \otimes \nabla u_0(x_j)\int _{\varOmega _o} h_i h_j,\nonumber \\ \end{aligned}$$
(69)

for \(i\ne j\) with \(i,j = 1,\ldots ,N\). Here there is no summation in the repeated indexes.

The truncated expansion \(\delta J({\overline{V}},\xi ,N)\), given by (65), is a quadratic form in \({\overline{V}}\). Hence, in order to minimize (65), we differentiate \(\delta J({\overline{V}},\xi ,N)\) with respect to \({\overline{V}}\) to get the first-order optimality condition

$$\begin{aligned} \langle D_{{\overline{V}}} \delta J, \beta \rangle = [d(\xi ) + H^{s}(\xi ) {\overline{V}}] \cdot \beta = 0, \quad \forall \beta \in \mathbb {R}^{2N}, \end{aligned}$$
(70)

with \( H^{s} = (H + H^{\top })/2\). From (70), we obtain the linear system

$$\begin{aligned} {\overline{V}} = - (H^{s}(\xi ))^{-1} d(\xi ). \end{aligned}$$
(71)

Note that, the quantity \({\overline{V}}\), solution of (71), becomes a function of the locations \(\xi \), i.e., \({\overline{V}} = {\overline{V}}(\xi )\). Let us replace \({\overline{V}} = {\overline{V}}(\xi )\) in (65) to obtain

$$\begin{aligned} \delta J({\overline{V}}(\xi ),\xi ,N) = \frac{1}{2} d(\xi ) \cdot {\overline{V}}(\xi ). \end{aligned}$$
(72)

Therefore, the optimal locations and the velocity field of each leakage, given by the pair (\(\xi ^\star \), \(V^\star \)), can be obtained by

$$\begin{aligned} \xi ^{\star } := \underset{\xi \in X}{\text {argmin}}\;\;\delta J({\overline{V}}(\xi ),\xi ,N) \quad \text {and} \quad {\overline{V}}^{\star } := {\overline{V}}(\xi ^{\star }), \end{aligned}$$
(73)

where X denotes the set of admissible locations for the pollution sources. Since the quantity \({\overline{V}}^\star \) is obtained, the velocity field \(V^{\star }\) at the point \(\xi ^\star \) can be calculated through the formula introduced in (64), for each leakage \(\omega _i^*\), \(i = 1,\ldots ,N\).

Fig. 5
figure 5

Example 3: Target (a) and result (b)

By construction, the resulting algorithm is analogous to the one developed in Sect.  5.1. The main difference between them is the reconstruction variable. In the first algorithm, we have reconstructed the size of each leakage which is a scalar quantity. In the second one, the unknown variable is a vector quantity, namely, the velocity at the leakage points.

5.2.1 Numerical Example

A numerical example is presented below. In this case, we demonstrate the capability of the algorithm in reconstructing the mean velocity field when multiple leakages are considered in the target domain.

Example 3

reconstruction of the velocity field for multiple leakages. In this example, the target domain contains two leakages of same size \(\varepsilon _1^*= \varepsilon _2^*= 0.03\) whose locations are given by the points \(x_1^*= (-0.1,0.3)\) and \(x_2^*= (0.1,-0.2)\). See Fig. 5a. The velocity of the leakage is given by the vector field \(V = (-5y^2 +5y,1)\). Therefore, the velocity evaluated at the center of each leakage is equal to \(V(x_1^*) = (1.05,1)\) and \(V(x_2^*) = (-1.2,1)\), respectively. The reconstruction is obtained from three measurements with the help of all Dirichlet data \(g_1\), \(g_2\) and \(g_3\). In this case, we successfully reconstructed the center of each pollution source, as can be seen in Fig. 5b. The velocities of the leakages were accurately obtained, namely, \(V(x_1^\star ) = (1.0418,0.9972)\) and \(V(x_2^\star ) = (-1.1946,0.9908)\).

6 Conclusions

As concluding remarks concerning the proposed reconstruction algorithms, we highlight some important points. Briefly, for a given velocity field V and number N of trial pollution sources, the reconstruction algorithm is able to find their locations \(\xi ^\star \) and sizes \(\alpha ^\star = \alpha (\xi ^\star )\) in one step. The proposed method approximates the unknown set of pollution sources by several balls which can be seen as a limitation of our approach. However, it can be used to get a good initial guess for more sophisticated iterative approaches based on level-sets methods [4, 8, 9, 19,20,21], for instance. Similarly, in the second algorithm, for a given number N of pollution sources whose sizes are assumed to be known, the reconstruction algorithm returns their locations \(\xi ^\star \) and velocities \(V^\star = V(\xi ^\star )\) in one step. In addition, such algorithms are independent of initial guess and they are capable of reconstructing the pollution sources in the presence of noise data. Finally, let us point out some interesting features/aspects of the type of algorithm proposed above: (a) when the number \(N^*\) of target subdomains is unknown, the algorithm can be started based on an assumption that there exists \(N>N^*\) subdomains and then we should find a number (\(N-N^*\)) of trial balls with negligible sizes [17]; (b) if the center of the target subdomain does not belong to the set of nodes of the sub-mesh, namely \(x^*\notin X\), the algorithm returns a location \(x^\star \) which is the closest to \(x^*\) [17]; and (c) since a combinatorial search over all the n-points of the set X has to be performed, this exhaustive search becomes rapidly infeasible for \(n\gg N\) as N increases [23].