Pollution Sources Reconstruction Based on the Topological Derivative Method

The topological derivative method is used to solve a pollution sources reconstruction problem governed by a steady-state convection-diffusion equation. To be more precise, we are dealing with a shape optimization problem which consists of reconstruction of a set of pollution sources in a fluid medium by measuring the concentration of the pollutants within some subregion of the reference domain. The shape functional measuring the misfit between the known data and solution of the state equation is minimized with respect to a set of ball-shaped geometrical subdomains representing the pollution sources. The necessary conditions for optimality are derived with the help of the topological derivative method which consists in expanding the shape functional asymptotically and then truncate it up to the second order term. The resulting expression is trivially minimized with respect to the parameters under consideration which leads to a noniterative second-order reconstruction algorithm. Two different cases are considered. Firstly, when the velocity of the leakages is given and we reconstruct the support of the unknown sources, including their locations and sizes. In the second case, we consider the size of the pollution sources to be known and find out the mean velocity of the leakages and their locations. Numerical examples are presented showing the capability of the proposed algorithm in reconstructing multiple pollution sources in both cases.

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 reallife 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 wellposedness. 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.

Problem Formulation
Let Ω ⊂ R 2 be an open and bounded domain with smooth boundary ∂Ω. Measurements of the pollutant concentration field are collected in a subdomain Ω o ⊂ Ω. We assume that there may be an unknown number (N * ∈ Z + ) of isolated pollution sources ω * i within the domain Ω. See Fig. 1a. Therefore, there is a set We assume that the polluting substance is leaking in an incompressible flow with velocity V ∈ C(Ω, R 2 ), such that V = (0, 0) at each ω * i ⊂ Ω, i = 1, . . . , N * . For a given Dirichlet data g imposed on ∂Ω, the resulting pollutant concentration z in Ω is observed in the subdomain Ω o . In this set up, the inverse problem consists in finding χ ω * such that the pollutant concentration z satisfies the following boundary value problem where the parameter χ ω * is defined as and the velocity V is considered as a null divergence vector field, namely, div V = 0 in ω * i , for i = 1, . . . , N * . Without loss of generality, we are considering only one boundary data g on ∂Ω. 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 Ω o . The realistic problem of pollution source reconstruction is much more complicated and requires further investigation. Now, for an initial guess χ ω of χ ω * , we consider the pollutant concentration field u to be the solution to the boundary value problem where The characteristic function χ ω * is unknown and hence z but we assume that z can be measured in Ω o . Thus, if we look for an appropriate χ ω * , then we wish u to agree with z in Ω o , i.e., we want u = z | Ωo . In addition, in such a problem, we assume that the velocity field V is known and then we reconstruct the support of ω * from partial measurements of the scalar field z taken on Ω 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 ω * i , i = 1, . . . , N * , then the inverse problem of finding χ ω * can be rewritten as a topology optimization problem, namely where z and u are the solutions of the boundary value problems (1) and (3), respectively, corresponding to the Dirichlet data g.

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 J ω (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 ω * , we start with the unperturbed domain by setting ω = ∅. See Fig. 1b. More precisely, we consider where u 0 be the solution of the unperturbed boundary value problem In this article, we are considering the topology optimization problem (5) for the ball-shaped pollution sources and hence we define the topologically perturbed counterpart of (7) by introducing N ∈ Z + number of small circular pollution sources B ε i (x i ) with center at x i ∈ Ω and radius ε i for i = 1, . . . , N . The respective geometric set can be denoted as where ξ = (x 1 , . . . , x N ) and ε = (ε 1 , . . . , ε N ). Moreover, we assume that The shape functional associated with the topologically perturbed domain is written as with u ε be the solution of the perturbed boundary value problem where the parameter χ ε is defined as By construction, we consider that the velocity field V satisfies the additional condition div V = 0 in B ε (ξ ). 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 As mentioned earlier, the topological derivatives measure the sensitivity of the shape functional with respect to the parameters (ε, ξ ) depending upon a set of small circular pollution sources B ε (ξ ). Therefore, our idea is to obtain the number, location and radius of the components of the set B ε (ξ ) in order to produce the best approximation to the true set of pollution sources ω * by using the concept of topological derivatives.
In this article, we are interested in expanding the shape functional J (u ε ) defined in (9) in power of ε. Therefore, we start by simplifying the difference between the perturbed shape functional J (u ε ) and its unperturbed counter-part J (u 0 ) defined in (9) and (6), respectively, as follows Since we are approximating the set of pollution sources to be reconstructed ω * by a number of circular balls B ε (ξ ), the idea is to expand the perturbed shape functional J (u ε ) with respect to the Lebesgue measure (volume) of the two-dimensional ball . . , N . To simplify the notation, we introduce the vector and the characteristic function Additionally, for i = 1, . . . , N , let us define the linear operators satisfying and For a positive real number R ε, we consider the ball 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 ε with respect to the parameters describing the small circular pollution sources as presented in Sect. 2 where for i, j = 1, . . . , N . The last term of the expansion (18), namelyũ ε , must satisfy the problem with such that, for i, j, l = 1, . . . , N , one has In order to simplify further analysis, we write h ε i as a sum of three functions The function p ε i is defined as Building on the expression (24), we construct q i , i = 1, . . . , N , in order to compensate for the discrepancies introduced by p ε i on ∂Ω. Therefore, we define Additionally, the residual functionh ε i is given bỹ Taking into account the decomposition (23) together with the expressions (24)-(26), we introduce an alternative representation for the function where with that is, for i = 1, . . . , N , q i solves the homogeneous boundary value problem Now, building on the expression forh ε i given by (26), we first define the vector quantity W i coming out from the Taylor expansion of the argument of ψ i as and then the functionh ε i can be split into the following sum wherep ε i is defined as the functionq i is such that and the residual termh ε i is given bỹ where The function p ε ii is such that where we have used expression (24) to obtain Additionally, the functions q ii andh ε i are given by 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

Asymptotic Expansion of the Shape Functional
Let us start by replacing (18) in (12) to obtain where On the basis of the boundary value problem for the adjoint state (41) together with the Green's formula, we obtain The closed formula given by expression (46) allows to replace the first three integrals over Ω o on the right-hand side of (42) by integrals over the ball B ε 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 ·∇u 0 , V i ·∇h ε j and V i ·∇v ε i j in (46), expression (42) can be rewritten as 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 = j, then we replace h ε j and v ε j j by the decompositions (23) and (36), respectively. Since the fourth integral in (47) is evaluated over Ω o , then the functions h ε i (and also h ε j ) are replaced by expression given in (27). Thus, the asymptotic expansion (47) can be expressed as Here, the six new remainders are defined as and In order to obtain the final form of the asymptotic expansion of the shape functional J (u ε ) from (48) we use the explicit representation of the functions p ε i ,p ε i and p ε ii , given by (24), (33) and (37), respectively, to compute their gradients in B ε i (x i ). Moreover, we also consider the Taylor's expansions of the functions ∇u 0 , ∇q i , ∇h j and v around the point x i up to produce a remainder of order o(|ε| 4 ) as done in [16]. Thus, we get The new remainder E 0 (ε) contains the residual terms arising from different Taylor's expansions arguing in a similar fashion to get the remainders E 10 (ε)-E 19 (ε) 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, . . . , N . Additionally, let V be the null divergence velocity field introduced in problem (1). Then, the asymptotic expansion for the topologically perturbed shape functional J (u ε ) defined in (9), with respect to where the remaining terms, given by (54), are such that |E(ε)| = o(|ε| 4 ), whose justification can be found in Appendix A.

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 Ω which contains a N * number of leakage points through 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 ξ and sizes α. 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 α of the pollution sources are known and determine their locations ξ and respective velocities V i of the leakages, i = 1, . . . , N .
Numerical experiments are conducted in the following scenario. The reference geometrical domain is given by a square Ω = (-0.5, 0.5) × (-0.5, 0.5) which is discretized with three-node finite elements. The mesh is generated from a grid of 160 × 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 Ω o which is a union of four small regions in the neighborhood of the corners of the square domain Ω. See Fig. 2a. The boundary ∂Ω of the reference domain is excited by imposing different Dirichlet data, namely, g 1 = x, g 2 = y and g 3 = x y.
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.
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 where z μ 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 Ω o by gray and the remaining domain by white colors.

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 J (u ε ) given by (53). After disregarding all terms of order o(|α| 2 ) in (53), the resulting truncated expansion, in its matrix form, is written as where the entries of the vector d ∈ R N and the matrix H ∈ R N ×N are given by and if i = j, respectively, for i, j = 1, . . . , N . Here there is no summation in the repeated indexes.
The truncated asymptotic expansion δ J (α, ξ, N ), given by (56), is a quadratic form in α. Hence, in order to minimize (56), we proceed with the derivative of δ J (α, ξ, N ) with respect to the variable α to write the first-order optimality condition with H s = (H + H )/2, where H denotes the transpose of the matrix H . From (60), we obtain the linear system Note that, the quantity α, solution of (61), becomes a function of the locations ξ , i.e., α = α(ξ ). Let us replace α = α(ξ ) in (56) to obtain Therefore, the optimal locations and sizes of the pollution sources, given by the pair of vectors (ξ , α ), can be obtained by 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].

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 ω * located at x * = (0.2, 0.1) whose radius is ε * = 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 ω * 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 ε * of ω * . 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. 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 0.3, 0). The sizes of the leakages are equal, namely, ε * i = 0.03, for i = 1, . . . , 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, ε i = 0.0299, for i = 1, . . . , 4.

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 J (u ε ) given by (53). Firstly, all the terms of order o(|α| 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 in order to work with the two-dimensional flow rate of the leakages at the points x i , i = 1, . . . , N . Taking into account these two steps, we obtain the truncated expansion where d, V ∈ R 2N and H ∈ R 2N ×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 V ∈ R 2N is given by V :  := (d 1 , . . . , d N ) where each entry d i is a vector in R 2 given by is given by and Here there is no summation in the repeated indexes. The truncated expansion δ J (V , ξ, N ), given by (65), is a quadratic form in V . Hence, in order to minimize (65), we differentiate δ J (V , ξ, N ) with respect to V to get the first-order optimality condition with H s = (H + H )/2. From (70), we obtain the linear system Note that, the quantity V , solution of (71), becomes a function of the locations ξ , i.e., Therefore, the optimal locations and the velocity field of each leakage, given by the pair (ξ , V ), can be obtained by where X denotes the set of admissible locations for the pollution sources. Since the quantity V is obtained, the velocity field V at the point ξ can be calculated through the formula introduced in (64), for each leakage ω * i , i = 1, . . . , N . 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.

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 ε * 1 = ε * 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 ) = (1.0418, 0.9972) and V (x 2 ) = (−1.1946, 0.9908).

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 ξ and sizes α = α(ξ ) 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 ξ and velocities V = V (ξ ) 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 * / ∈ X , the algorithm returns a location x 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 N as N increases [23].

Conflicts of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Justification for the Main Result
The justification for the result stated through Theorem 1 is presented in two steps.
Firstly, we prove a priori estimates related to the auxiliary states N . Then, in the final part of this section, the previously obtained results are used to estimate the remainders left in the asymptotic expansion of the topologically perturbed cost functional. These estimates justify our topological asymptotic expansion (53).

A.1 A Priori Estimates
In order to simplify the presentation, we denote all the constants independent of ε and i as C for i = 1, . . . , N , whose value changes according to the place it is used. Furthermore, we use the fact that α i ∼ ε 2 i , when necessary.

Lemma 1 Let Ω be an open and bounded domain in R 2 and let B ε be a ball of radius
ε, such that B ε ⊂ Ω. Then, for a function ϕ ∈ H 1 (Ω), the following estimate holds true with 0 < δ < 1 and the constant C independent of the small parameter ε.
Proof By takingh ε i as a test function in the weak formulation of (26), we get From the Cauchy-Schwarz inequality and the interior elliptic regularity of the function u 0 , there exists a positive constant C independent of ε and i such that where we have used Lemma 1, with ϕ =h ε i , to obtain the last inequality. From the Poincaré inequality, we have Combining the results obtained in (88) and (87), we get (77). By considering the decomposition (23) as well as the fact that p ε i and q i are independent of ε in Ω o , we can write From (77), we have (78). From the decomposition (23) and the triangular inequality, we have as well as By (24), simple calculation gives Notice that, we can use the estimate of Lemma 1 for ϕ =h ε i together with the result obtained in (77) to get Using the estimates (92) in (90), as well as the interior elliptic regularity of q i together with (94), we obtain (79) and (80), respectively. From (91), we have Using the estimates (93) in (95), as well as the interior elliptic regularity of q i together with (77), we obtain (81) and (82), respectively. Now, let us take h ε i as a test function in the weak form of (19) 1 to obtain From the Cauchy-Schwarz inequality and the interior elliptic regularity of the function u 0 , there exists a positive constant C independent of ε and i such that where we have used (79) to obtain the last inequality in (97). From the Poincaré inequality, we have Combining the results obtained in (97) and (98), we get (83).
By takingh ε i as a test function in the weak formulation of (35), we get From the Cauchy-Schwarz inequality, there exists a positive constant C independent of ε and i such that where we have used the fact that T i (x) = O( x − x i 2 ) and the estimate to obtain the last inequality. Now, we set ϕ =h ε i in Lemma 1 and we use the respective estimate to obtain, from (100), that From the Poincaré inequality, we have Combining the results obtained in (101) and (102), we get the Result (84). By considering the decomposition (32), we can write From (33), simple calculation gives The required Result (85) is obtained by taking into account the interior elliptic regularity of the function q j as well as the estimates (84) for any 0 < δ i , δ j < 1, with i, j = 1, . . . , N .
Proof By takingṽ ε ii as a test function in the weak formulation of (40), we have From the Cauchy-Schwarz inequality, there exists a positive constant C independent of ε and i such that where we have used the interior elliptic regularity of the function q i and Lemma 2, Result (77), to obtain the last inequality. By considering the estimate of Lemma 1 for From the Poincaré inequality, we have Combining the results obtained in (115) and (116), we get the Result (105). By taking v ε i j in the weak formulation of (19) 2 , we obtain From the Cauchy-Schwarz inequality, there exists a positive constant C independent of ε and i such that From (118), in the particular case where i = j, we have where we have used Lemma 2, Result (81), as well as the result of Lemma 1 with ϕ = v ε ii , to obtain the last inequality. In addition, we also have from (118), for i = j, where we have used Lemma 2, Result (82), as well as the result of Lemma 1 with ϕ = v ε i j . From the Poincaré inequality, we have Combining the results obtained in (119) and (121), for i = j, we obtain the Result (106). Analogously, we obtain the Result (107) from (120) and (121). By considering the decomposition (36) as well as the fact that p ε ii and q ii are inde- From (105), we have (108). From the decomposition (36) and the triangular inequality, we have By (37), simple calculation gives Using the estimates (124) in (123), as well as the interior elliptic regularity of q ii together with (105), we obtain (109) and (110), respectively. By taking w ε i j as a test function in the weak formulation of (19) 3 , we have From the Cauchy-Schwarz inequality, there exists a positive constant C independent of ε and i such that From (126), in the particular case where i = j, we have where we have used (109) and the estimate of Lemma 1 with ϕ = w ε ii to obtain the last inequality. In addition, we also have from (126), for i = j, that where we have used (110) and the estimate of Lemma 1 with ϕ = w ε i j . From the Poincaré inequality, we have Combining the results obtained in (127) and (129), for i = j, we obtain the Result (111). Analogously, we obtain the Result (112) from (128) and (129).

A.2 Estimates of the Remainders
We will successively prove that |E (ε)| = o(|ε| 4 ) for = 1, . . . , 15, where |ε| := ε 1 +· · ·+ε N . For simplicity, we use the symbol C to denote any constant independent of ε. Remainder E 0 (ε) can be trivially bounded by using similar arguments as in [16], leading to |E 0 (ε)| = o(|ε| 4 ). The estimates of the other remainders are more involved and will be presented in detail in what follows. In order to obtain the estimate of each remainder, we start by using the Cauchy-Schwarz inequality and then we use the appropriate lemmas of Sect. A.1 as well as additional results related to some functions introduced in Section 3. Proceeding in this way, we obtain for any 2/3 < δ < 1, where we have used Lemma 4; for any 0 < δ < 1, where we have used Lemma 2, result (78), and Lemma 3, Results (107) and (108); for any 0 < δ < 1, where we have used Lemma 2, Result (78), and Lemma 3, Results (111) and (112); for any 0 < δ < 1, where we have used Lemma 2, Result (78), and Lemma 4; for any 0 < δ < 1, where we have used Lemma 3, Results (107) and (108); for any 1/6 < δ < 1, where we have used Lemma 3, Results (107), (108), (111) and (112); for any 0 < δ < 1, where we have used Lemma 3, Results (107)-(108), and Lemma 4; for any 0 < δ < 1, where we have used Lemma 3, Results (111) and (112); for any 0 < δ < 1, where we have used Lemma 3, Results (111)-(112), and Lemma 4; for any 0 < δ < 1, where we have used the interior elliptic regularity of the functions v and q i as well as the estimate V i ·∇h ε i L 2 (B ε i ) ≤ C hε i H 1 (Ω) together with Lemma 2, Result (84); for any 0 < δ < 1, where we have used the interior elliptic regularity of the function v as well as Lemma 2, Result (85); where we have used the interior elliptic regularity of the functions q j j and v; for any 1/2 < δ < 1, where we have used the interior elliptic regularity of the function v as well as the estimate V i · ∇ṽ ε j j L 2 (B ε i ) ≤ C ṽ ε j j H 1 (Ω) together with Lemma 3, Result (105); where we have used the interior elliptic regularity of the function v and the explicit solution of p ε j j , given by (37), to estimate V i · ∇ p ε j j L 2 (B ε i ) = O(ε i ), for i = j; for any 0 < δ < 1, where we have used the fact that h i L 2 (Ω o ) = O(1) and Lemma 2, Result (77).