Source Identification of a Chemical Incident in an Urban Area

Abstract

This work deals aims to present a methodology for source identification of chemical incidents in urban areas. We propose an approximation of the problem within the framework of the optimal control theory and we provide an algorithm for its numerical resolution. Finally, we analyze the validity of the algorithm in several academic situations.

1. Introduction

The existence of weapons of mass destruction represents an increase in the potential threat to peace in different areas of the world. Some of these weapons, with capacity for killing and bringing significant harm to numerous humans, may not only be in the hands of great powers, but also of regional powers and even terrorist organizations.

Regarding the use of nuclear, biological, chemical (NBC) weapons, this possibility must be always considered in asymmetric armed conflicts, where their use is more likely than in conflicts between great powers. An adversary with NBC capacity may introduce agents, materials and weapons at any time in a more or less indirect way. The launch and dispersal of NBC agents or materials can be carried out with different means such as missiles, aircrafts of all kinds, field artillery, difficult-to-detect aerosols, etc. Once an NBC incident occurs, it is vitally important to be able to predict the danger area to evacuate the civilian population and alert the mobilized units in that area [1]. The methods used to estimate the hazard area can be classified into:

• Simplified methods: These are preliminary estimates based on the characteristics of the incident and meteorological data. They are usually carried out by hand by a trained person.
• Improved methods: These are automatic or manual estimates that are usually made taking into account the type of incident, the place where it occurred, and the weather conditions. They are more accurate than previous ones and update as weather conditions change.
• Methods based on mathematical simulation: These are fully automatic methods that estimate the hazard area by numerical simulation, from the type of incident, meteorological data, and space-time domain information.

Currently, the most widely used mathematical models to simulate the evolution of a chemical agent are based on the Gaussian models. They are closed-form analytical solutions of the classic advection-diffusion equation

$$\frac{\partial c}{\partial t}+\nabla ·\left(c\mathbf{u}\right)=\nabla ·(\mathbf{K}\nabla c)+S,\Omega \times ]0,T[,$$

(1)

derived from some suitable hypotheses [2,3]. Specifically, it is assumed that:

Hypothesis 1.

The concentration of the chemical agent c is not dependent on time (the solution is steady state) and, specifically, $\partial c/\partial t=0$.

Hypothesis 2.

The incident occurs at a fixed point $\mathbf{b}=(0,0,H)$, where the chemical agent is emitted at a constant rate $Q>0$, in such a way that the source term is given by ${\displaystyle S\left(\mathbf{x}\right)=Q{\delta }_{\mathbf{b}}\left(\mathbf{x}\right)}$, where ${\delta }_{\mathbf{b}}\left(\mathbf{x}\right)$ is the Dirac delta at point $\mathbf{b}$.

Hypothesis 3.

The wind velocity field $\mathbf{u}$ is constant and aligned with the positive $x_1$-axis, that is $\mathbf{u}=(u,0,0)$ for some constant $u\ge 0$.

Hypothesis 4.

The diffusion coefficients are the same in all directions and only depend on the downwind distance $x_1$, that is, the diffusion matrix $\mathbf{K}$ is given by $\mathbf{K}=K\left(x_1\right)I_3$, where $K\left(x_1\right)$ is a real function and $I_3$ is the identity matrix.

Hypothesis 5.

The effect of diffusion on the $x_1$-axis is neglected (dominant convection), that is, $K\left(x_1\right){\partial }^{2}c/\partial {x_1}^2=0$.

Hypothesis 6.

Topographic variations and obstacles (trees, buildings, etc.) are neglected, in such a way that the space domain is $\Omega =[0,\infty )\times (-\infty ,\infty )\times [0,\infty )$.

Hypothesis 7.

There is no chemical agent for $x_1<0$, and the chemical agent does not penetrate the soil. Thus, the following boundary conditions can be considered

$$\begin{array}{c}c(0,x_2,x_3)=c(+\infty ,x_2,x_3)=0,\forall x_2\in \mathbb{R},x_3\ge 0,\hfill \\ c(x_1,-\infty ,x_3)=c(x_1,+\infty ,x_3)=0,\forall x_1,x_3\ge 0,\hfill \\ K\left(x_1\right)\partial c/\partial x_3(x_1,x_2,0)=0,c(x_1,x_2,+\infty )=0,\forall x_1\ge 0,x_2\in \mathbb{R}.\hfill \end{array}$$

(2)

Under these hypotheses, the general Equation (1) is rewritten as

$$u\frac{\partial c}{\partial x_1}=K\left(\frac{{\partial }^{2}c}{\partial {x_2}^2}+\frac{{\partial }^{2}c}{\partial {x_3}^2}\right)+Q{\delta }_{\mathbf{b}}\left(\mathbf{x}\right),\Omega$$

(3)

and completed with condition (2). The explicit solution of systems (2) and (3) is obtained by using the Laplace transform [3], and it is known as the Gaussian-plume model:

$$c(r,x_2,x_3)=\frac{Q\exp \left(-\frac{{\left(x_2\right)}^2}{4r}\right)\left(\exp \left(-\frac{{(x_3-H)}^2}{4r}\right)+\exp \left(-\frac{{(x_3+H)}^2}{4r}\right)\right)}{4\pi ur},$$

(4)

where

$$r=\frac{1}{u}{\int }_0^{x_1}K\left(\xi \right)d\xi .$$

If the chemical incident is instantaneous and only occurs at the initial time $t=0$, the Gaussian-plume model is not suitable. In this situation, hypotheses (H${}_3$)–(H${}_7$) hold, but (H${}_1$) and (H${}_2$) must be replaced respectively by:

Hypothesis 8.

There is no chemical agent before the incident, that is,

$$c(x_1,x_2,x_3,0)=0,\forall x_1,x_3\ge 0,x_2\in \mathbb{R}.$$

(5)

Hypothesis 9.

The incident occurs at the initial time $t=0$, and at a fixed point $\mathbf{b}=(0,0,H)$, in such a way that the source term is given by ${\displaystyle S\left(\mathbf{x}\right)=Q_T{\delta }_{\mathbf{b}}\left(\mathbf{x}\right){\delta }_0\left(t\right)}$, where $Q_T$ is the total amount of chemical agent released.

Under these hypotheses, the general Equation (1) is rewritten as

$$\frac{\partial c}{\partial t}+u\frac{\partial c}{\partial x_1}=K\left(\frac{{\partial }^{2}c}{\partial {x_2}^2}+\frac{{\partial }^{2}c}{\partial {x_3}^2}\right)+Q_T{\delta }_{\mathbf{b}}\left(\mathbf{x}\right){\delta }_0\left(t\right),\Omega \times ]0,T[,$$

and completed with the initial condition (5), and boundary conditions (2) formulated as $\forall t>0$. The system is solved again by using the Laplace transform, and the explicit solution

$${\displaystyle c(r,x_2,x_3,t)=\frac{Q_T\exp \left(-\frac{{(x_1-ut)}^2+{\left(x_2\right)}^2}{4r}\right)\left(\exp \left(-\frac{{(x_3-H)}^2}{4r}\right)+\exp \left(-\frac{{(x_3+H)}^2}{4r}\right)\right)}{8{\left(\pi r\right)}^{3/2}}}$$

(6)

is referred as the Gaussian-puff model.

Previous considerations (hypotheses (H${}_3$)–(H${}_7$)) clearly reveal the limitations of Gaussian models (4) and (6) to simulate the evolution of a chemical agent in an urban region, and consequently, to determine the hazard area if the chemical incident occurs in an urban domain. The main objective of this paper is just to develop a novel method to deal with chemical incidents in urban areas. To avoid the limitations of the Gaussian models, we will deal with the general equation, Equation (1), to characterize the source of the chemical agent, from measurements made at atmospheric monitoring stations located at different points of the city. The scientific literature on this subject is very rich, and there are many papers dealing not only with the mathematical study of inverse source problems [4,5,6], but also with interesting environmental applications in surface water [7,8,9,10], in groundwater [11,12] and in the atmosphere [2,8]. In this paper, the problem will be studied within the framework of optimal control problem of partial differential equations (PDEs). Taking advantage of previous works of the authors on the control of the urban heat island [13,14], and thinking about a 3D urban domain, the main novelty of the model proposed in this paper is that the classic advection-diffusion equation, Equation (1), will be completed with a reaction term depending on the air temperature, and combined with a 3D microclimatic model to simulate the wind velocity between buildings and the heat transfer between air, soil and buildings. Additionally, taking into account that the admissible set may be nonconnected, the inverse problems will be formulated and solved within the framework of mixed integer nonlinear programming (MINLP).

This paper is organized as follows. The mathematical model proposed to simulate the chemical agent evolution is presented in Section 2.1, and completed in Appendix A, where the 3D microclimatic model is detailed. From this model, in Section 2.2 MINLP is used to formulate the inverse problems within the framework of optimal control of PDEs. A completed numerical method to solve these problems is detailed in Section 2.3, and numerical results are presented and discussed in Section 3. Finally, some conclusions are summarized in Section 4.

2. Materials and Methods

2.1. Numerical Simulation: The State Model

In this section, we present the 3D mathematical model that we will use in the numerical resolution of the problem. We will consider a three-dimensional bounded domain ${\Omega }_A^{3D}\subset {\mathbb{R}}^3$ corresponding to an urban area, where $\partial {\Omega }_A^{3D}$ is the boundary of said domain, this is walls and ceilings of buildings, floor and, additionally, fictitious borders that delimit our domain (see Figure 1). We will denote by ${\Gamma }_A^{IN}$ the boundary corresponding to an incoming air flow. We will assume that the air temperature can affect the concentration of the chemical agent and that, eventually, there may be sedimentation effects. Therefore, the evolution of the concentration of the chemical agent $c_A$ (gr/m${}^3$) will be given by the solution of the following equation:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial c_A}{\partial t}+{\mathbf{u}}_A·\nabla c_A+w_C\frac{\partial c_A}{\partial x_3}-\nabla ·(K_C\nabla c_A)=F_C+G({\theta }_A,c_A),{\Omega }_A^{3D}\times ]0,T[,}\hfill \\ {\displaystyle c_A=c_A^{IN},{\Gamma }_A^{IN}\times ]0,T[,}\hfill \\ {\displaystyle K_C\frac{\partial c_A}{\partial {\mathbf{n}}_A}=0,(\partial {\Omega }_A^{3D}\setminus {\Gamma }_A^{IN})\times ]0,T[,}\hfill \\ {\displaystyle c_A\left(0\right)=c_A^0,{\Omega }_A^{3D}.}\hfill \end{array}\right.$$

(7)

where $w_C$ (m/s) is the sedimentation velocity (constant), ${\mathbf{u}}_A$ (m/s) is the air velocity, ${\theta }_A$ (K) is the air temperature, $F_C$ (gr/m${}^3$ s) is the source term, G (gr/m${}^3$ s) represents the influence of air temperature on the chemical agent (in order to simplify the model we will assume that $G({\theta }_A,c_A)=G\left({\theta }_A\right)c_A$, with $G(·)$ a negative function), $K_A$ (m${}^2$/s) is the diffusion constant, $c_A^{IN}$ (gr/m${}^3$) is the concentration of the chemical agent on the inlet boundary ${\Gamma }_A^{IN}$, and $c_A^0$ (gr/m${}^3$) is the initial concentration of the chemical agent. We must mention that $\nabla ·\left({\mathbf{u}}_{\mathbf{A}}c\right)=\nabla ·{\mathbf{u}}_A{c}_A+{\mathbf{u}}_A·\nabla c_A=\mathbf{u}·\nabla c_A$ since we are assuming that $\nabla ·{\mathbf{u}}_A=0$ (when we are considering air layers close to the ground, it is usually considered that the air behaves like an incompressible fluid). Regarding the source term, it is frequently considered [15] to be of the form:

$$F_C(\mathbf{x},t)=Q_C\left(t\right){\delta }_{{\mathbf{b}}_C}\left(\mathbf{x}\right),$$

(8)

where $Q_C\left(t\right)$ (gr/s) is the release rate and ${\mathbf{b}}_C$ is the point in which the source term is located. In the case of an instantaneous release, we will consider the following term:

$$F_C(\mathbf{x},t)=Q_C{\delta }_0\left(t\right){\delta }_{{\mathbf{b}}_C}\left(\mathbf{x}\right),$$

(9)

where $Q_C$ (gr) is the total amount of chemical agent released at time $t=0$. In the case we have an instantaneous release, we will rewrite (9) in terms of an initial condition:

$$c_A\left(0\right)=Q_C{\delta }_{{\mathbf{b}}_C}\left(\mathbf{x}\right).$$

The temperature ${\theta }_A$ (K) and air velocity ${\mathbf{u}}_A$ (m/s) will be obtained by solving a microclimatic model in which we will take into account the temperature of the soil ${\theta }_S$ (K) and buildings ${\theta }_B$ [K] (see Appendix A).

2.2. Optimal Control: The Inverse Problem

In this section, we will formulate the inverse problem consisting of the characterization of the source term associated with the chemical incident from a set of measurements taken in the urban area. For this, we will formulate the inverse problem by means of an optimal control problem [16]. Let us start by establishing the following notations:

• We will assume that the possible release point locations (${\mathbf{b}}_C$) are in a bounded region ${\mathcal{U}}_{ad}$ given by the union of M convex closed and bounded subsets (admissible release zones):

  $${\mathcal{U}}_{ad}=\bigcup _{k=1}^M{Z}_{ad}^k,$$

  where $\mathrm{int}\left(Z_{ad}^i\right)\cap \mathrm{int}\left(Z_{ad}^j\right)=\varnothing$, $\forall i\ne j$, and

  $$Z_{ad}^k=[l_1^k,u_1^k]\times [l_2^k,u_2^k]\times [l_3^k,u_3^k],k=1,\dots ,M,$$

  where $l_i^k$ and $u_i^k$ are, respectively, the lower and upper bounds for the $x_i$ coordinate of $Z_{ad}^k$, $i=1,2,3$. Let us empathize that the set ${\mathcal{U}}_{ad}$ can be nonconnected.
• $Q_C$ is the release rate or the total amount of chemical agent released. We will asume that $Q_C\in {\mathcal{V}}_{ad}$, where ${\mathcal{V}}_{ad}=\{Q\in X:0\le Q\left(t\right)\le Q_{max},\forall t\in [0,T]\}$ if the source term is given by (8) or ${\mathcal{V}}_{ad}=\{Q\in X:0\le Q_C\le Q_{max}\}$ if the source term is given by (9), where $X=L^2(0,T)$ in the first case and $X=\mathbb{R}$ in the second one.
• $\{{\tilde{c}}_i^n:i=1,\dots ,N_P,n=1,\dots ,N_T\}$ is the set of measurements taken in the urban area at points $({\mathbf{x}}_i,t_n)\in (\overline{{\Omega }_A^{3D}}\setminus {\Gamma }_A^{IN})\times [0,T]$, $i=1,\dots ,N_P$, $n=1,\dots ,N_T$.
• We consider the following objective function:

  $$J({\mathbf{b}}_C,Q_C)=\sum _{i=1}^{N_P}\sum _{n=1}^{N_T}{\left(c_A({\mathbf{x}}_i,t_n)-{\tilde{c}}_i^n\right)}^2.$$
  We must observe that to evaluate the previous function at each element $({\mathbf{b}}_C,Q_C)$ it is necessary to solve the state equation, Equation (7).

We will study the following optimal control problems [16].

• Problem 1. Estimation of the release point: We will estimate the release point assuming that the emission rate ${\tilde{Q}}_C$ of the chemical agent is known:

  $$\underset{{\mathbf{b}}_C\in {\mathcal{U}}_{ad}}{min}J({\mathbf{b}}_C,{\tilde{Q}}_C).$$

  (10)
• Problem 2. Estimation of the release rate: We will estimate the release rate (or the total amount) $Q_C$ assuming that the release point ${\tilde{\mathbf{b}}}_C$ is known:

  $$\underset{Q_C\in {\mathcal{V}}_{ad}}{min}J({\tilde{\mathbf{b}}}_C,Q_C).$$

  (11)
• Problem 3. Estimation of the release point and rate:

  $$\underset{({\mathbf{b}}_C,Q_C)\in {\mathcal{U}}_{ad}\times {\mathcal{V}}_{ad}}{min}J({\mathbf{b}}_C,Q_C).$$

  (12)

Let us observe that problems 1 and 3, Equations (10) and (12), respectively, can be formulated as Nonlinear Mixed Integer Programming Problems (MINLPs); if we introduce an integer variable ${\mathbf{y}}_C\in {\{0,1\}}^M$ such that $y_C^k=1$, if the release point is in the zone $Z_{ad}^k$, and $y_C^k=0$ in other cases. Taking into account the previous variable, we can reformulate problem 3, Equation (12), in the following classical framework of MINLPs (problem 1, Equation (10), is analogous):

$$\begin{array}{cc}\hfill {\displaystyle \underset{({\mathbf{b}}_C,Q_C,{\mathbf{y}}_C)\in {\mathbb{R}}^3\times {\mathcal{V}}_{ad}\times {\{0,1\}}^M}{min}}& \hfill {\displaystyle J({\mathbf{b}}_C,Q_C)}\\ \hfill {\displaystyle \mathrm{s}.\mathrm{t}.}& \hfill {\displaystyle \mathbf{h}({\mathbf{b}}_C,{\mathbf{y}}_C)}& \le & \mathbf{0}\hfill \\ \hfill & \hfill {\displaystyle A{\mathbf{y}}_C}& =& \mathbf{c},\hfill \end{array}$$

(13)

where:

• $A\in {\mathcal{M}}_{m_2\times M}$ with $m_2=1$ is such that $A{\mathbf{y}}_C={\sum }_{k=1}^M{y}_C^k$, therefore:

  $$A=\left(\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right),$$

  and $\mathbf{c}=1\in {\mathbb{R}}^{m_2}$. Observe that the constraint $A{\mathbf{y}}_C=\mathbf{c}$ implies that there is only one 1 in the control vector ${\mathbf{y}}_C$ with the rest of its components being equal to 0. We will denote by

  $${\mathcal{Y}}_{ad}=\{{\mathbf{y}}_C\in {\{0,1\}}^M:A{\mathbf{y}}_C=\mathbf{c}\}$$

  the admissible control set for the discrete variable ${\mathbf{y}}_C$.
• Given ${\mathbf{b}}_C=(p_C^1,p_C^2,p_C^3)\in {\mathbb{R}}^3$ and ${\mathbf{y}}_C\in {\mathcal{Y}}_{ad}$,

  $$\mathbf{h}({\mathbf{b}}_C,{\mathbf{y}}_C)=\left(\begin{array}{c}{\displaystyle {\mathbf{h}}^1({\mathbf{b}}_C,{\mathbf{y}}_C)}\\ {\displaystyle {\mathbf{h}}^2({\mathbf{b}}_C,{\mathbf{y}}_C)}\\ {\displaystyle {\mathbf{h}}^3({\mathbf{b}}_C,{\mathbf{y}}_C)}\end{array}\right)\in {\mathbb{R}}^6,$$

  being

  $${\mathbf{h}}^j({\mathbf{b}}_C,{\mathbf{y}}_C)=\left(\begin{array}{c}{\displaystyle \sum _{k=1}^M{y}_C^k(l_j^k-p_C^j)}\\ {\displaystyle \sum _{k=1}^M{y}_C^k(p_C^j-u_j^k)}\end{array}\right)\in {\mathbb{R}}^2,j=1,2,3.$$
  Observe that if ${\mathbf{y}}_C\in {\mathcal{Y}}_{ad}$ and ${\mathbf{b}}_C\in {\mathbb{R}}^3$ satisfies $\mathbf{h}({\mathbf{b}}_C,{\mathbf{y}}_C)\le \mathbf{0}$, then ${\mathbf{b}}_C\in {\mathcal{U}}_{ad}$.

If we fix the variable ${\mathbf{y}}_C^{*}\in {\mathcal{Y}}_{ad}$ we obtain a classical NLP problem:

$$\begin{array}{cc}\hfill {\displaystyle \underset{({\mathbf{b}}_C,Q_C)\in {\mathbb{R}}^3\times {\mathcal{V}}_{ad}}{min}}& \hfill {\displaystyle J({\mathbf{b}}_C,Q_C)}\\ \hfill {\displaystyle \mathrm{s}.\mathrm{t}.}& \hfill {\displaystyle \mathbf{h}({\mathbf{b}}_C,{\mathbf{y}}_C^{*})}& \le & \mathbf{0}.\hfill \end{array}$$

(14)

We denote by $({\mathbf{b}}_C^{*},Q_C^{*})$ a solution of the optimal control problem (14) associated with the discrete variable ${\mathbf{y}}_C^{*}$. Thus, if we consider the following set:

$$\begin{array}{c}{\displaystyle \mathcal{F}\left({\mathcal{Y}}_{ad}\right)=\{({\mathbf{b}}_C^{*},Q_C^{*},{\mathbf{y}}_C^{*})\in {\mathcal{U}}_{ad}\times {\mathcal{V}}_{ad}\times {\mathcal{Y}}_{ad}:}\\ {\displaystyle ({\mathbf{b}}_{\mathrm{C}}^{*},{\mathbf{Q}}_{\mathrm{C}}^{*})\mathrm{solution}\mathrm{of}\left(14\right)\mathrm{associated}\mathrm{with}{\mathbf{y}}_{\mathrm{C}}^{*}\},}\end{array}$$

we can solve problem (13) by taking $({\mathbf{b}}_C,{\mathbf{Q}}_C,{\mathbf{y}}_C)\in \mathcal{F}\left({\mathcal{Y}}_{ad}\right)$ such that:

$$J({\mathbf{b}}_C,{\mathbf{Q}}_C)=min\left\{J({\mathbf{b}}_C^{*},{\mathbf{Q}}_C^{*}):({\mathbf{b}}_C^{*},{\mathbf{Q}}_C^{*},{\mathbf{y}}_C^{*})\in \mathcal{F}\left({\mathcal{Y}}_{ad}\right)\right\}.$$

(15)

For low values of M, the size of the set $F\left(Y_{ad}\right)$ is small, and the MINLP problem, Equation (15), can be easily solved by an exhaustive search. For large-size problems, more appropriate methods, such as Branch and Bound, Generalized Benders Decomposition or external approximation (see, for instance, [17,18,19]) must be used. In any case, a quick method for solving the NLP problem, Equation (14), is the key for solving the MINLP problem, Equation (15). Consequently, the next section is devoted to present a numerical method for solving the NLP problem, Equation (14).

2.3. Numerical Resolution

To solve the NLP problem, Equation (14), we will use an algorithm of interior points, more specifically, IPOPT [20]. The use of this class of interior point algorithm requires, at least, the evaluation of the cost functional and the constraints, the evaluation of the gradient of the cost functional and the Jacobian matrix associated with the constraints. Since the constraints are linear, the problem lies in calculating the cost function and its gradient. In this section, we will detail how to carry out these computations.

The first step is the numerical resolution of the state equation, Equation (7). In order to achieve this, we will propose a space-time discretization based on the method of the characteristics and the finite element method [21,22]. Spatial and time discretizations have been performed in the scientific software FreeFem ++ [23]. For simplicity in the notations, and without loss of generality, we will assume that that the sedimentation rate is incorporated in the term ${\mathbf{u}}_A·\nabla c_A$ and we will use the notation $\mathbf{u}$ instead of ${\mathbf{u}}_A$.

Let us consider $N+1$ points ${\left\{t^n\right\}}_{n=0}^N$ in the interval $[0,T]$ such that:

• ${\displaystyle t^0=0}$,
• ${\displaystyle t^N=T}$,
• ${\displaystyle t^{n+1}-t^n=\Delta t,\forall n=0,\dots ,N-1}$.

We define ${\displaystyle \alpha =\frac{1}{\Delta t}}$ and we consider the material derivative for a scalar field c:

$$\frac{Dc}{Dt}(\mathbf{x},t)=\frac{\partial }{\partial t}c(X(\mathbf{x},t),t)=\frac{\partial c}{\partial t}(\mathbf{x},t)+\mathbf{u}(\mathbf{x},t)·\nabla c(\mathbf{x},t),$$

where ${\displaystyle \frac{\partial X}{\partial t}(\mathbf{x},t)=\mathbf{u}(\mathbf{x},t)}$. We can consider the following approximation for the material derivative in a time $t^{n+1}$:

$$\frac{Dc}{Dt}\left(t^{n+1}\right)\simeq \alpha (c^{n+1}-c^n\circ X_h^n),$$

with $X_h^n\left(\mathbf{x}\right)=X(\mathbf{x},t_{n+1},t_n)$ being the solution of the following initial value problem:

$$\left\{\begin{array}{c}{\displaystyle \frac{dX}{d\tau }=\mathbf{u}(X(\mathbf{x},t,\tau ),\tau ).}\hfill \\ {\displaystyle X(\mathbf{x},t,t)=\mathbf{x}.}\hfill \end{array}\right.$$

Thus, $c^n\circ X_h^n\simeq c^n(\mathbf{x}-{\mathbf{u}}_n\left(\mathbf{x}\right)\Delta t)$. This approximation, as we will see later, is very important for obtaining the gradient of the objective function.

So, given $c_A^0$, we compute ${\left\{c_A^{n+1}\right\}}_{n=0}^{N-1}$ solving the following equation:

$$\left\{\begin{array}{c}{\displaystyle \alpha c_A^{n+1}-\nabla ·(K_C\nabla c_A^{n+1})=F_C^{n+1}+G\left({\theta }_A^{n+1}\right)c_A^{n+1}+\alpha (c_A^n\circ X_h^n),{\Omega }_A^{3D},}\hfill \\ {\displaystyle c_A^{n+1}=0,{\Gamma }_A^{IN},}\hfill \\ {\displaystyle K_C\frac{\partial c_A^{n+1}}{\partial {\mathbf{n}}_A}=0,\partial {\Omega }_A^{3D}\setminus {\Gamma }_A^{IN}.}\hfill \end{array}\right.$$

(16)

For spatial discretization, we will assume that the domain ${\Omega }_A^{3D}$ is polyhedral and we consider ${\left\{{\tau }_h^A\right\}}_{h>0}$ as a family of regular meshes of the domain ${\Omega }_A^{3D}$. We define the following finite element space:

$$X_h^A=\{z\in \mathcal{C}\left(\overline{{\Omega }_A^{3D}}\right):z_{{|}_{\mathcal{T}}}\in {\mathcal{P}}_1\left(\mathcal{T}\right),\forall \mathcal{T}\in {\tau }_h^A,z_{{|}_{{\Gamma }_A^{IN}}}=0\}.$$

Thus, the fully discretized problem consists of ${\left\{c_A^{n+1}\right\}}_{n=0}^{N-1}\subset X_h^A$ solving the following variational formulation:

$$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}{c}_A^{n+1}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla c_A^{n+1}·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)c_A^{n+1}zd\mathbf{x}}\\ {\displaystyle ={\int }_{{\Omega }_A^{3D}}{F}_C^{n+1}zd\mathbf{x}+\alpha {\int }_{{\Omega }_A^{3D}}(c_A^n\circ X_h^n)zd\mathbf{x},\forall z\in W_h^A.}\end{array}$$

(17)

Remark 1.

Taking into account the definition of the Dirac Delta as a distribution, the first addition in the second term of Equation (17) can be computed using the following formula:

$${\int }_{{\Omega }_A^{3D}}{F}_C^{n+1}zd\mathbf{x}={\int }_{\Omega }{Q}_C\left(t^{n+1}\right){\delta }_{{\mathbf{b}}_C}\left(\mathbf{x}\right)zd\mathbf{x}=Q_C\left(t^{n+1}\right)z\left({\mathbf{b}}_C\right).$$

(18)

However, if we want to model a case in which the emission is not strictly punctual, the following approximation of the Dirac Delta can be considered:

$$\rho \left(\mathbf{x}\right)=\left\{\begin{array}{cc}Cexp\left(\frac{1}{{\parallel \mathbf{x}\parallel }^2-1}\right)\hfill & \parallel x\parallel <1\hfill \\ 0\hfill & \parallel x\parallel \ge 1\hfill \end{array}\right.$$

where $C>$ is such that:

$${\int }_{{\mathbb{R}}^n}\rho \left(\mathbf{x}\right)dx=1.$$

From the previous function, we can define

$${\phi }_{\mathbf{q},ϵ}\left(\mathbf{x}\right)=\frac{1}{{ϵ}^n}\rho \left(\frac{\mathbf{x}-\mathbf{q}}{ϵ}\right).$$

which converges to ${\delta }_{\mathbf{q}}\left(\mathbf{x}\right)=\delta (\mathbf{x}-\mathbf{q})$ when $ϵ\to 0$. Taking into account the previous sequence,

$$F_C(\mathbf{x},t)\simeq Q_C\left(t\right){\phi }_{{\mathbf{b}}_C,ϵ}\left(\mathbf{x}\right),$$

thus

$${\int }_{{\Omega }_A^{3D}}{F}_C^{n+1}zd\mathbf{x}\simeq Q_C^{n+1}{\int }_{{\Omega }_A^{3D}}{\phi }_{{\mathbf{b}}_C,ϵ}\left(\mathbf{x}\right)zdx.$$

(19)

If $F_C$ is given by (8), let us consider the following discretization for the control variable $Q_C\left(t\right)$:

$$Q_C\left(t\right)=Q_C^1{\chi }_{[t_0,t_1]}\left(t\right)+\sum _{n=2}^N{Q}_C^n{\chi }_{(t_{n-1},t_n]}\left(t\right),$$

where ${\mathbf{Q}}_C=(Q_C^1,Q_C^2,\dots ,Q_C^n)\in {\mathbb{R}}^n$. Thus,

$${\mathbf{s}}_C=(\underset{{\mathbf{Q}}_C}{\underbrace{Q_C^1,Q_C^2,\dots ,Q_C^n}},\underset{{\mathbf{b}}_C}{\underbrace{p_C^1,p_C^2,p_C^3}},\underset{{\mathbf{y}}_C}{\underbrace{y_C^1,y_C^2,\dots ,y_C^M}})\in {\mathbb{R}}^n\times {\mathbb{R}}^3\times {\{0,1\}}^M$$

denotes the discrete global control.

On the contrary, if $F_C$ is given by (9), we obtain

$${\mathbf{s}}_C=(\underset{{\mathbf{Q}}_C}{\underbrace{Q_C}},\underset{{\mathbf{b}}_C}{\underbrace{p_C^1,p_C^2,p_C^3}},\underset{{\mathbf{y}}_C}{\underbrace{y_C^1,y_C^2,\dots ,y_C^M}})\in {\mathbb{R}}^n\times {\mathbb{R}}^3\times {\{0,1\}}^M$$

as the global control variable.

In order to simplify the notations, we will assume that the measurements are made at the times associated with the time discretization. Thus,

$$J({\mathbf{b}}_C,{\mathbf{Q}}_C)=\sum _{i=1}^{N_P}\sum _{n=1}^N{\left(c_A^n\left({\mathbf{x}}_i\right)-{\tilde{c}}_i^n\right)}^2,$$

is the discretized objective function, where ${\left\{c_A^n\right\}}_{n=1}^N$ are the solutions of the fully discretized state equation, Equation (17).

To calculate the gradient of the cost functional, we can use the linearized equations or the adjoint state equations. In the computations that we present below, we will assume that $F_C$ is given by (8) and the modifications for treating case (9) are straightforward. We will denote by ${\delta }_{\mathbf{Q}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)\left(\delta \mathbf{Q}\right)$ the directional derivative of J with respect to ${\mathbf{Q}}_C$ in the direction $\delta \mathbf{Q}$ and by ${\delta }_{\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)\left(\delta \mathbf{b}\right)$ the directional derivative of J with respect to $\mathbf{b}$ in the direction $\delta \mathbf{b}$. Next, we present the expressions for the previous directional derivatives considering the linearized equations and the adjoint state equations, considering the two approaches for the computation of Dirac delta.

• Directional derivative of J with respect to $\mathbf{Q}$ using the linearized equations:

  $${\delta }_{\mathbf{Q}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)\left(\delta \mathbf{Q}\right)=2\sum _{i=1}^{N_p}\sum _{n=1}^N(c_A^n\left({\mathbf{x}}_i\right)-{\tilde{c}}_i^n){\delta }_{\mathbf{Q}}{c}_A^n\left({\mathbf{x}}_i\right),$$

  (20)

  where, given ${\delta }_{\mathbf{Q}}{c}_A^0=0$, ${\delta }_{\mathbf{Q}}{c}_A^n\in W_h^A$, $n=0,\dots ,N-1$, is the solution to:

  $$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}{\delta }_{\mathbf{Q}}{c}_A^{n+1}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla {\delta }_{\mathbf{Q}}{c}_A^{n+1}·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right){\delta }_{\mathbf{Q}}{c}_A^{n+1}zd\mathbf{x}}\\ {\displaystyle =\delta Q^{n+1}{\int }_{{\Omega }_A^{3D}}{\phi }_{{\mathbf{b}}_C,ϵ}zd\mathbf{x}+\alpha {\int }_{{\Omega }_A^{3D}}({\delta }_{\mathbf{Q}}{c}_A^n\circ X_h^n)zd\mathbf{x},\forall z\in W_h^A,}\end{array}$$

  (21)

  in case (19). In case (18) ${\delta }_{\mathbf{Q}}{c}_A^n\in W_h^A$ is the solution to:

  $$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}{\delta }_{\mathbf{Q}}{c}_A^{n+1}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla {\delta }_{\mathbf{Q}}{c}_A^{n+1}·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right){\delta }_{\mathbf{Q}}{c}_A^{n+1}zd\mathbf{x}}\\ {\displaystyle =\delta Q^{n+1}z\left({\mathbf{b}}_C\right)+\alpha {\int }_{{\Omega }_A^{3D}}({\delta }_{\mathbf{Q}}{c}_A^n\circ X_h^n)zd\mathbf{x},\forall z\in W_h^A.}\end{array}$$

  (22)
• Directional derivative of J with respect to $\mathbf{b}$ using the linearized equations:

  $${\delta }_{\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)\left(\delta \mathbf{b}\right)=2\sum _{i=1}^{N_p}\sum _{n=1}^N(c_A^n\left({\mathbf{x}}_i\right)-{\tilde{c}}_i^n){\delta }_{\mathbf{b}}{c}_A^n\left({\mathbf{x}}_i\right),$$

  (23)

  where, given ${\delta }_{\mathbf{b}}{c}_A^0=0$, ${\delta }_{\mathbf{b}}{c}_A^{n+1}\in W_h^A$, $n=0,\dots ,N-1$, is the solution to

  $$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}{\delta }_{\mathbf{b}}{c}_A^{n+1}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla {\delta }_{\mathbf{b}}{c}_A^{n+1}·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right){\delta }_{\mathbf{b}}{c}_A^{n+1}zd\mathbf{x}}\\ {\displaystyle =Q_C^{n+1}{\int }_{{\Omega }_A^{3D}}{\delta }_{\mathbf{b}}{\phi }_{\mathbf{b},ϵ}\left(\delta \mathbf{b}\right)zd\mathbf{x}+\alpha {\int }_{{\Omega }_A^{3D}}({\delta }_{\mathbf{b}}{c}_A^n\circ X_h^n)zd\mathbf{x},\forall z\in W_h^A,}\end{array}$$

  (24)

  in case (19). In case (18) ${\delta }_{\mathbf{b}}{c}_A^{n+1}\in W_h^A$ is the solution to:

  $$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}{\delta }_{\mathbf{b}}{c}_A^{n+1}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla {\delta }_{\mathbf{b}}{c}_A^{n+1}·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right){\delta }_{\mathbf{b}}{c}_A^{n+1}zd\mathbf{x}}\\ {\displaystyle =Q_C^{n+1}\nabla z\left({\mathbf{b}}_C\right)·\delta \mathbf{b}+\alpha {\int }_{{\Omega }_A^{3D}}({\delta }_{\mathbf{b}}{c}_A^n\circ X_h^n)zd\mathbf{x},\forall z\in W_h^A.}\end{array}$$

  (25)

Remark 2.

It should be noted that in Equation (25) the term $Q_C^{n+1}\nabla z\left({\mathbf{b}}_C\right)·\delta \mathbf{b}$ appears. This term is the result of the approximation of the derivative of the Dirac delta [24,25]. The basic idea is to consider a polynomial approximation ${\delta }_h(·,\mathbf{b})$ of the Dirac delta $\delta (\mathbf{x}-\mathbf{b})$. Indeed, let us assume that $\mathbf{b}\in T\in {\tau }_h^A$ and let $\widehat{\varphi }\left(\widehat{\mathbf{x}}\right)\in {\mathcal{P}}_1\left(\widehat{T}\right)$ such that

$${\int }_{\widehat{T}}\widehat{p}\left(\widehat{\mathbf{x}}\right)\widehat{\varphi }\left(\widehat{\mathbf{x}}\right)=\widehat{p}\left(F_T\left(\mathbf{b}\right)\right),\forall \widehat{p}\in {\mathcal{P}}_1\left(\widehat{T}\right),$$

(26)

where $\widehat{T}$ is the reference element and $F_T:T\to \widehat{T}$, with $F_T\left(\mathbf{x}\right)\in {\left[{\mathcal{P}}_1\left(T\right)\right]}^3$. To obtain the above polynomial, let us consider a basis $\widehat{B}=\{{\widehat{p}}_1,{\widehat{p}}_2,\dots ,{\widehat{p}}_L\}$ of the vector space ${\mathcal{P}}_1\left(\widehat{T}\right)$. We know that (26) is equivalent to:

$${\int }_{\widehat{T}}{\widehat{p}}_k\left(\widehat{\mathbf{x}}\right)\widehat{\varphi }\left(\widehat{\mathbf{x}}\right)={\widehat{p}}_k\left(F_T\left(\mathbf{b}\right)\right),\forall k=1,\dots ,L.$$

We denote by ${\widehat{\varphi }}_B={({\widehat{\alpha }}_1,\dots ,{\widehat{\alpha }}_L)}^t$ the coordinates of $\widehat{\phi }$ on the basis of $\widehat{B}$ ($\widehat{\varphi }={\sum }_{l=1}^L{\widehat{\alpha }}_l{\widehat{p}}_l$). We know that ${\widehat{\varphi }}_B$ is the solution to the following linear system:

$$G_{\widehat{B}}{\widehat{\varphi }}_{\widehat{B}}={\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),$$

where ${\widehat{p}}_{\widehat{B}}={({\widehat{p}}_1,\dots ,{\widehat{p}}_L)}^t$ and $G_{\widehat{B}}\in Sym\left(M_{L\times L}\left(\mathbb{R}\right)\right)$ is the Gram matrix associated with $\widehat{B}$:

$${\left[G_{\widehat{B}}\right]}_{i,j}={\int }_{\widehat{T}}{\widehat{p}}_i\left(\widehat{\mathbf{x}}\right){\widehat{p}}_j\left(\widehat{\mathbf{x}}\right)d\widehat{\mathbf{x}},i,j=1,\dots ,L.$$

Thus,

$$\begin{array}{ccc}\hfill {\displaystyle \widehat{\varphi }\left(\widehat{\mathbf{x}}\right)}& =& {\displaystyle <{\widehat{\varphi }}_B,{\widehat{p}}_{\widehat{B}}\left(\widehat{\mathbf{x}}\right)>}\hfill \\ & =& {\displaystyle <G_{\widehat{B}}^{-1}{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),{\widehat{p}}_{\widehat{B}}\left(\widehat{\mathbf{x}}\right)>}\hfill \\ \hfill & =& {\displaystyle <{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{\widehat{p}}_{\widehat{B}}\left(\widehat{\mathbf{x}}\right)>.}\hfill \end{array}$$

(27)

Now, we define

$${\delta }_h(\mathbf{x},\mathbf{b})=\left\{\begin{array}{ccc}|J_{\mathbf{x}}{F}_T\left(\mathbf{x}\right)|\widehat{\varphi }\left(F_T\left(\mathbf{x}\right)\right)\hfill & if\hfill & \mathbf{x}\in T,\hfill \\ 0\hfill & if\hfill & \mathbf{x}\notin T,\hfill \end{array}\right.$$

where $J_{\mathbf{x}}{F}_T\left(\mathbf{x}\right)$ is the Jacobian matrix of $F_T$. If we use expression (27):

$${\delta }_h(\mathbf{x},\mathbf{b})=\left\{\begin{array}{ccc}|J_{\mathbf{x}}{F}_T\left(\mathbf{x}\right)|<{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{x}\right)\right)>\hfill & if\hfill & \mathbf{x}\in T,\hfill \\ 0\hfill & if\hfill & \mathbf{x}\notin T\hfill \end{array}\right.$$

Taking into account the above definition, given an element $z_h\in W_h^A$:

$$\begin{array}{ccc}\hfill {\displaystyle {\int }_{{\Omega }_A^{3D}}{\delta }_h(\mathbf{x},\mathbf{b})z_h\left(\mathbf{x}\right)d\mathbf{x}}& =& {\displaystyle {\int }_T\left|J_{\mathbf{x}}{F}_T\left(\mathbf{x}\right)\right|<{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{x}\right)\right)>z_h\left(\mathbf{x}\right)d\mathbf{x}}\hfill \\ & =& {\displaystyle <{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{\int }_T\left|J_{\mathbf{x}}{F}_T\left(\mathbf{x}\right)\right|{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{x}\right)\right)z_h\left(\mathbf{x}\right)d\mathbf{x}>}\hfill \\ & =& {\displaystyle <{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{\int }_{\widehat{T}}{\widehat{p}}_{\widehat{B}}\left(\widehat{\mathbf{x}}\right){\widehat{z}}_h\left(\widehat{\mathbf{x}}\right)d\widehat{\mathbf{x}}>,}\hfill \end{array}$$

where ${\widehat{z}}_h\left(\widehat{\mathbf{x}}\right)=z_h\left(F_T^{-1}\left(\widehat{\mathbf{x}}\right)\right)$. Now, ${\widehat{z}}_h={{\widehat{z}}_{h_B}}^t{\widehat{p}}_{\widehat{B}}$; therefore:

$$\begin{array}{ccc}\hfill {\displaystyle <{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{\int }_{\widehat{T}}{\widehat{p}}_{\widehat{B}}\left(\widehat{\mathbf{x}}\right){\widehat{z}}_h\left(\widehat{\mathbf{x}}\right)d\widehat{\mathbf{x}}>}& =& {\displaystyle <{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{\int }_{\widehat{T}}{\widehat{p}}_{\widehat{B}}\left(\widehat{\mathbf{x}}\right){{\widehat{z}}_{h_B}}^t{\widehat{p}}_{\widehat{B}}\left(\widehat{\mathbf{x}}\right)d\widehat{\mathbf{x}}>}\hfill \\ & =& {\displaystyle <{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{G}_{\widehat{B}}{\widehat{z}}_{h_B}>}\hfill \\ & =& {\displaystyle {\widehat{z}}_{h_B}^t{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right)}\hfill \\ & =& {\widehat{z}}_h\left(F_T\left(\mathbf{b}\right)\right)\hfill \\ & =& {\displaystyle z_h\left(\mathbf{b}\right).}\hfill \end{array}$$

Thus:

$${\int }_{{\Omega }_A^{3D}}{\delta }_h(\mathbf{x},\mathbf{b})z_h\left(\mathbf{x}\right)d\mathbf{x}=z_h\left(\mathbf{b}\right).$$

Finally:

$$\begin{array}{ccc}\hfill {\displaystyle {\int }_{{\Omega }_A^{3D}}{\nabla }_{\mathbf{b}}{\delta }_h(\mathbf{x},\mathbf{b})z_h\left(\mathbf{x}\right)d\mathbf{x}}& =& {\nabla }_{\mathbf{b}}<{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right),G_{\widehat{B}}^{-1}{G}_{\widehat{B}}{\widehat{z}}_{h_B}>\hfill \\ & =& {\displaystyle {\widehat{z}}_{h,B}^t{J}_{\widehat{\mathbf{x}}}{\widehat{p}}_{\widehat{B}}\left(F_T\left(\mathbf{b}\right)\right)J_{\mathbf{x}}{F}_T\left(\mathbf{b}\right)}\hfill \\ & =& {\displaystyle \nabla z_h\left(\mathbf{b}\right).}\hfill \end{array}$$

In view of the expressions (20) and (23), we observe that to calculate the gradient of the objective function using the linearized equations, we have to solve N times (21) or (22) and 3 times (24) or (25). Indeed,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial J}{\partial Q^1}({\mathbf{Q}}_C,{\mathbf{b}}_C)}& =& {\displaystyle {\delta }_{\mathbf{Q}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(1,0,\dots ,0),}\hfill \\ \hfill {\displaystyle \frac{\partial J}{\partial Q^2}({\mathbf{Q}}_C,{\mathbf{b}}_C)}& =& {\displaystyle {\delta }_{\mathbf{Q}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(0,1,\dots ,0),}\hfill \\ \hfill ⋮& & ⋮\hfill \\ \hfill {\displaystyle \frac{\partial J}{\partial Q^N}({\mathbf{Q}}_C,{\mathbf{b}}_C)}& =& {\displaystyle {\delta }_{\mathbf{Q}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(0,0,\dots ,1),}\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial J}{\partial b^1}({\mathbf{Q}}_C,{\mathbf{b}}_C)}& =& {\displaystyle {\delta }_{\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(1,0,0),}\hfill \\ \hfill {\displaystyle \frac{\partial J}{\partial b^2}({\mathbf{Q}}_C,{\mathbf{b}}_C)}& =& {\displaystyle {\delta }_{\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(0,1,0),}\hfill \\ \hfill {\displaystyle \frac{\partial J}{\partial b^3}({\mathbf{Q}}_C,{\mathbf{b}}_C)}& =& {\displaystyle {\delta }_{\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(0,0,1).}\hfill \end{array}$$

Therefore, to calculate, for example, $\frac{\partial J}{\partial Q^1}({\mathbf{Q}}_C,{\mathbf{b}}_C)$, we have to solve (21) or (22) taking ${\delta }_{\mathbf{Q}}=(1,0,\dots ,0)$, for computing $\frac{\partial J}{\partial Q^2}({\mathbf{Q}}_C,{\mathbf{b}}_C)$ and we have to solve (21) or (22) taking ${\delta }_{\mathbf{Q}}=(0,1,\dots ,0)$, and so on.

Next, we will see that if we use the adjoint state equations it will only be necessary to solve one equation to calculate the gradient of the cost functional.

• Directional derivative of J with respect to $\mathbf{Q}$ and $\mathbf{b}$ using the adjoint state equation. On the one hand,

  $${\delta }_{\mathbf{Q},\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(\delta \mathbf{Q},\delta {\mathbf{b}}_C)=2\sum _{i=1}^{N_p}\sum _{n=1}^N(c_A^n\left({\mathbf{x}}_i\right)-{\tilde{c}}_i^n){\delta }_{\mathbf{Q},\mathbf{b}}{c}_A^n\left({\mathbf{x}}_i\right),$$

  (28)

  where, given $\delta c_A^0={\delta }_{\mathbf{Q},\mathbf{b}}{c}_A^0=0$, $\delta c_A^{n+1}={\delta }_{\mathbf{Q},\mathbf{b}}{c}_A^{n+1}\in W_h^A$, $n=0,\dots ,N-1$ is the solution to:

  $$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}\delta c_A^{n+1}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla \delta c_A^{n+1}·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)\delta c_A^{n+1}zd\mathbf{x}}\\ {\displaystyle =\delta Q^{n+1}{\int }_{{\Omega }_A^{3D}}{\phi }_{{\mathbf{b}}_C,ϵ}zd\mathbf{x}+Q_C^{n+1}{\int }_{{\Omega }_A^{3D}}{\delta }_{{\mathbf{b}}_C}{\phi }_{\mathbf{b},ϵ}\left(\delta \mathbf{b}\right)zd\mathbf{x}}\\ {\displaystyle +\alpha {\int }_{{\Omega }_A^{3D}}(\delta c_A^n\circ X_h^n)zd\mathbf{x},\forall z\in W_h^A,}\end{array}$$

  (29)

  in case (19). In case (18), $\delta c_A^{n+1}={\delta }_{\mathbf{Q},\mathbf{b}}{c}_A^{n+1}\in W_h^A$ is such that:

  $$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}\delta c_A^{n+1}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla \delta c_A^{n+1}·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)\delta c_A^{n+1}zd\mathbf{x}}\\ {\displaystyle =\delta Q^{n+1}z\left({\mathbf{b}}_C\right)+Q_C^{n+1}\nabla z\left({\mathbf{b}}_C\right)·\delta \mathbf{b}+\alpha {\int }_{{\Omega }_A^{3D}}(\delta c_A^n\circ X_h^n)zd\mathbf{x},\forall z\in W_h^A.}\end{array}$$

  (30)
  Now, we will see how we can obtain the equations for the adjoint state. Let us consider ${\left\{r_A^n\right\}}_{n=0}^N\subset W_h^A$ such that $r_A^N=0$ and let us take, for each $n=1,\dots ,N-1$, $r_A^n$ as a test function in (29) or (30):

  $$\begin{array}{c}{\displaystyle \sum _{n=0}^{N-1}\left\{\alpha {\int }_{{\Omega }_A^{3D}}\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla c_A^{n+1}·\nabla r_A^{n}d\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}\right\}}\\ {\displaystyle =\sum _{n=0}^{N-1}\{\alpha {\int }_{{\Omega }_A^{3D}}(\delta c_A^n\circ X_h^n)r_A^{n}d\mathbf{x}+\delta Q^{n+1}{\int }_{{\Omega }_A^{3D}}{\phi }_{{\mathbf{b}}_C,ϵ}{r}_A^{n}d\mathbf{x}}\\ {\displaystyle +Q_C^{n+1}{\int }_{{\Omega }_A^{3D}}{\delta }_{{\mathbf{b}}_C}{\phi }_{\mathbf{b},ϵ}\left(\delta \mathbf{b}\right)r_A^{n}d\mathbf{x}\},}\end{array}$$

  in case (19) and for (18):

  $$\begin{array}{c}{\displaystyle \sum _{n=0}^{N-1}\left\{\alpha {\int }_{{\Omega }_A^{3D}}\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla c_A^{n+1}·\nabla r_A^{n}d\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}\right\}}\\ {\displaystyle =\sum _{n=0}^{N-1}\left\{\alpha {\int }_{{\Omega }_A^{3D}}(\delta c_A^n\circ X_h^n)r_A^{n}d\mathbf{x}+\delta Q^{n+1}{r}_A^n\left({\mathbf{b}}_C\right)+Q_C^{n+1}\nabla r_A^n\left({\mathbf{b}}_C\right)·\delta \mathbf{b}\right\}.}\end{array}$$
  Taking into account that $\delta c_A^0=r_A^N=0$,

  $$\begin{array}{c}{\displaystyle \sum _{n=0}^{N-1}\left\{\alpha {\int }_{{\Omega }_A^{3D}}\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla c_A^{n+1}·\nabla r_A^{n}d\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}\right\}}\\ {\displaystyle =\sum _{n=0}^{N-1}\{\alpha {\int }_{{\Omega }_A^{3D}}(\delta c_A^{n+1}\circ X_h^{n+1})r_A^{n+1}d\mathbf{x}+\delta Q^{n+1}{\int }_{{\Omega }_A^{3D}}{\phi }_{{\mathbf{b}}_C,ϵ}{r}_A^{n}d\mathbf{x}}\\ {\displaystyle +Q_C^{n+1}{\int }_{{\Omega }_A^{3D}}{\delta }_{{\mathbf{b}}_C}{\phi }_{\mathbf{b},ϵ}\left(\delta \mathbf{b}\right)r_A^{n}d\mathbf{x}\}}\end{array}$$

  (31)

  in case (19) and for (18): -5.0cm0cm

  $$\begin{array}{c}{\displaystyle \sum _{n=0}^{N-1}\left\{\alpha {\int }_{{\Omega }_A^{3D}}\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla c_A^{n+1}·\nabla r_A^{n}d\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)\delta c_A^{n+1}{r}_A^{n}d\mathbf{x}\right\}}\\ {\displaystyle =\sum _{n=0}^{N-1}\left\{\alpha {\int }_{{\Omega }_A^{3D}}(\delta c_A^{n+1}\circ X_h^{n+1})r_A^{n+1}d\mathbf{x}+\delta Q^{n+1}{r}_A^n\left({\mathbf{b}}_C\right)+Q_C^{n+1}\nabla r_A^n\left({\mathbf{b}}_C\right)·\delta \mathbf{b}\right\}.}\end{array}$$

  (32)
  Therefore, if we define $r_A^n\in W_h^A$, $n=N-1,\dots ,0$, as the solution to:

  $$\begin{array}{c}{\displaystyle \alpha {\int }_{{\Omega }_A^{3D}}{r}_A^{n}zd\mathbf{x}+{\int }_{{\Omega }_A^{3D}}{K}_C\nabla r_A^n·\nabla zd\mathbf{x}-{\int }_{{\Omega }_A^{3D}}G\left({\theta }_A^{n+1}\right)r_A^{n}zd\mathbf{x}}\\ {\displaystyle =\alpha {\int }_{{\Omega }_A^{3D}}(z\circ X_h^{n+1})r_A^{n+1}d\mathbf{x}+2\sum _{i=1}^{N_P}(c_A^{n+1}\left({\mathbf{x}}_i\right)-{\tilde{c}}_i^{n+1})z\left({\mathbf{x}}_i\right),}\end{array}$$

  (33)

  we know that:

  $$\begin{array}{ccc}\hfill {\displaystyle {\delta }_{\mathbf{Q},\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(\delta \mathbf{Q},\delta {\mathbf{b}}_C)}& =& {\displaystyle \sum _{n=0}^{N-1}\{\delta Q^{n+1}{\int }_{{\Omega }_A^{3D}}{\phi }_{{\mathbf{b}}_C,ϵ}{r}_A^{n}d\mathbf{x}}\hfill \\ & & {\displaystyle +Q_C^{n+1}{\int }_{{\Omega }_A^{3D}}{\delta }_{{\mathbf{b}}_C}{\phi }_{\mathbf{b},ϵ}\left(\delta \mathbf{b}\right)r_A^{n}d\mathbf{x}\},}\hfill \end{array}$$

  (34)

  in case we take an approximation of the Dirac delta (19) and if we consider (18):

  $${\delta }_{\mathbf{Q},\mathbf{b}}J({\mathbf{Q}}_C,{\mathbf{b}}_C)(\delta \mathbf{Q},\delta {\mathbf{b}}_C)=\sum _{n=0}^{N-1}\left\{\delta Q^{n+1}{r}_A^n\left({\mathbf{b}}_C\right)+Q_C^{n+1}\nabla r_A^n\left({\mathbf{b}}_C\right)·\delta \mathbf{b}\right\}.$$

  (35)

It should be noted that the equation for the discrete adjoint state (33) is valid for any choice of approximation of the Dirac delta. The considered approximation of the Dirac delta appears in the expression of the gradient of cost functionals (34) or (35). Furthermore, to calculate the gradient of the cost functional it is only necessary to solve the equation for the adjoint state once. Indeed, given ${\left\{r_A^n\right\}}_{n=0}^N\subset W_h^A$, with $r_A^N=0$, the solution to (33) is

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial Q^k}({\mathbf{Q}}_C,{\mathbf{b}}_C)}=& {\displaystyle {\int }_{{\Omega }_A^{3D}}{\phi }_{{\mathbf{b}}_C,ϵ}{r}_A^{k-1}d\mathbf{x},k=1,\dots ,N,}\hfill \end{array}$$

(36a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial b^k}({\mathbf{Q}}_C,{\mathbf{b}}_C)=}& {\displaystyle \sum _{n=0}^{N-1}{Q}_C^{n+1}{\int }_{{\Omega }_A^{3D}}\frac{\partial {\phi }_{\mathbf{b},ϵ}}{\partial b^k}\left({\mathbf{b}}_C\right)r_A^{n}d\mathbf{x},k=1,2,3,}\hfill \end{array}$$

(36b)

in this case we consider (19); for (18):

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial Q^k}({\mathbf{Q}}_C,{\mathbf{b}}_C)=}& {\displaystyle r_A^{k-1}\left({\mathbf{b}}_C\right),k=1,\dots ,N,}\hfill \end{array}$$

(37a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial J}{\partial b^k}({\mathbf{Q}}_C,{\mathbf{b}}_C)=}& {\displaystyle \sum _{n=0}^{N-1}{Q}_C^{n+1}\frac{\partial r_A^n}{\partial x_i}\left({\mathbf{b}}_C\right),k=1,2,3.}\hfill \end{array}$$

(37b)

3. Results and Discussion

In this section, we will present the results that we obtained in the numerical simulations. All the numerical simulations were carried out with scientific software FreeFem++ [23] interfaced with IPOPT [20] on a 2019 MacBook Pro (2.5 GHz Intel Core i5 with four kernels).

We considered a scale three-dimensional mesh composed of nine buildings with heights of, respectively, 8, 5, 4, 5, 6, 8, 8, 5 and 4 m (back to front and left to right), with the geometrical configuration presented in Figure 2 (the depth of the soil considered is 3 m).

Associated with the previous geometrical configuration, we have considered the domain occupied by the air (effective computational domain for the control problem); in our case, we considered the upper boundary 10 m from the ground, see Figure 3.

In order to obtain the solution of the microclimate model (see Appendix A), we have considered the parameters listed in

Table 1. Microclimate model parameters.

Coefficient	Air	Asphalt	Buildings
Density	${\rho }_A=1.16\times {10}^3$	${\rho }_C=2.11\times {10}^6$	${\rho }_B=2.3\times {10}^6$
Specific heat	$c{p}_A=1.007$	$c{p}_C=0.92$	$c{p}_B=0.88$
Conductivity	${\alpha }_A=0.05$	${\alpha }_C=0.06$	${\alpha }_B=0.06$
Emissivity		${ϵ}_C=0.95$	${ϵ}_B=0.95$
Albedo		$a_C=0.05$	$a_B=0.08$

.

The convective heat transfer coefficients for the corresponding interfaces were $h_C^A=100$ on ${\Gamma }_C^A$, $h_G^A=100$ on ${\Gamma }_G^A$, $h_C^B=100$ on ${\Gamma }_C^B$, $h_W^A=100$ on ${\Gamma }_B^W$, and $h_R^A=1$ on ${\Gamma }_B^R$.

To compute the radiation temperatures appearing in the heat equations for soil and buildings, we assumed that $R_{sw,net}(\mathbf{x},t)=(R{M}_{sw,dir}+R{M}_{sw,diff})$$\sigma (\mathbf{x},t)$, and $R_{lw,dow}(\mathbf{x},t)=R{M}_{lw,dow}\sigma (\mathbf{x},t)$, where, for our particular problem, we considered $R{M}_{sw,dir}=650$ Wm${}^{-2}$, $R{M}_{sw,diff}=350$ Wm${}^{-2}$, and $R{M}_{lw,dow}=450$ Wm${}^{-2}$. The function $\sigma (\mathbf{x},t)\in [0,1]$ models the attenuation of the previous maximum values, taking into account the movement of the sun: effect of shadows, night and day, and so on. In our simplified case, we have assumed that $\sigma (\mathbf{x},t)=max\{sin(2\pi t/86,400),0\},\forall (\mathbf{x},t)\in ({\Gamma }_C^A\cup {\Gamma }_G^A\cup {\Gamma }_B^R)\times [0,T]$ (that is, over soil and roofs we considered no attenuation due to shadows and radiation only depending on time). Finally, we considered the initial values ${\mathbf{u}}_A^0=(0,0)$ m s${}^{-1},{\theta }_A^0={\theta }_S^0={\theta }_B^0=300K$, and the boundary conditions ${\mathbf{u}}_A^{IN}=(0,{10}^{-4},0)$ m s${}^{-1},{\theta }_A^{IN}=300K$.

We also considered that $G=0$, $c_A^0=0$ gr m${}^{-3}$ and $K_C={10}^{-1}$ m${}^2$ s${}^{-1}$ in Equation (16). For the time discretization, we have considered a final time T = 10,800 s and $\Delta t=900$ s ($N=12$ time steps). In Figure 4, we can see the air velocity ${\mathbf{u}}_A$ and the air temperature ${\theta }_A$ on the boundary $\partial {\Omega }_A^{3D}$ at time T = 10,800 s.

We will assume that in each building there are three sensors placed at 1, 2 and 3 m from the ground (see Figure 5). Therefore, we have a total of $N_p=27$ measurement points.

We also assumed that we have measurements in each time step ($N_T=12$). To validate the methodology proposed in this work, we have generated artificial measurements taking as the release point ${\mathbf{b}}_R=(2,12,6)\in Z_{ad}^2=[1.2,31.8]\times [10.2,13.8]\times [4.2,11.8]$ and ${\mathbf{Q}}_R=\{c_R({\mathbf{x}}_i,t_n):i=1,\dots ,N_p,n=1,\dots ,N_T\}$, where $c_R$ is the solution of the discretized state equation, Equation (17), associated with the release point ${\mathbf{b}}_R$ and the release rate $Q_R\left(t\right)=3(1+sin(2\pi t/86,400-2\pi \Delta t/86,400))$ gr m${}^{-3}$ s${}^{-1}$, $t\in [0,T]$, such that $Q_R\left(t\right)\in {\mathcal{V}}_{ad}=[0,10]$. For instance, Figure 6 shows the concentrations in the measurement points placed at 3 m from the ground, if the delta approximation (19) is used. In this case, the distribution of the chemical agent at T = 10,800 s can be seen in Figure 7.

Therefore, the main objective of this section is to show how the methodology we propose allows one to recover the release point ${\mathbf{b}}_R$ and the discharge rate ${\mathbf{Q}}_R$ using artificial measures (see Figure 6).

We have obtained the following results taking a convergence tolerance for the IPOPT algorithm equal to ${10}^{-10}$:

• Problem 1, Equation (10). Estimation of the release point: In this case, we fix the release rate to ${\mathbf{Q}}_R$ and we try to recover the release point ${\mathbf{b}}_R$ using the artificial measurements. To this end, we start from the initial control ${\mathbf{b}}_0=[1.2,10.2,4.2]\in Z_{ad}^2$ and run the optimization algorithm in the following cases:

  (a)

  Delta approximation (19) and linearized Equation (23): 22 iterations, optimal cost $6.3054131\times {10}^{-22}$, CPU time 6502 s;

  (b)

  Delta approximation (19) and adjoint state Equation (36b): 32 iterations, optimal cost $3.2106931\times {10}^{-21}$, CPU time 4766 s;

  (c)

  Delta definition (18) and linearized Equation (23): 37 iterations, optimal cost $1.3769812\times {10}^{-20}$, CPU time 1779 s;

  (d)

  Delta definition (18) and adjoint state Equation (37b): 35 iterations, optimal cost $8.3016154\times {10}^{-21}$, CPU time 1958 s.

  In all the above cases, the release point used for the generation of artificial data is recovered by the algorithm.
• Problem 2, Equation (11). Estimation of the release rate: In this case, we fix the release point to ${\mathbf{b}}_R$ and we try to recover the release rate ${\mathbf{Q}}_R$ using the artificial measurements. We start from the initial control ${\mathbf{Q}}_0=\mathbf{0}\in {\mathcal{V}}_{ad}$ and run the optimization algorithm in the following cases:

  (a)

  Delta approximation (19) and linearized Equation (20): 26 iterations, optimal cost $1.8307661\times {10}^{-11}$, CPU time $11,177$ s.

  (b)

  Delta approximation (19) and adjoint state Equation (36a): 29 iterations, optimal cost $1.8307668\times {10}^{-11}$, CPU time 4696 s.

  (c)

  Delta definition (18) and linearized Equation (20): 25 iterations, optimal cost $1.8564142\times {10}^{-19}$, CPU time 2798 s.

  (d)

  Delta definition (18) and adjoint state Equation (37a): 27 iterations, optimal cost $1.0237655\times {10}^{-19}$, CPU time 709 s.

  In all the above cases, the release rate used for the generation of artificial data is recovered by the algorithm.
• Problem 3, Equation (12). Estimation of the release point and rate. We try to recover the release point ${\mathbf{b}}_R$ and the release rate ${\mathbf{Q}}_R$ using the artificial measurements. We start from the initial control ${\mathbf{b}}_0=[1.2,10.2,4.2]\in Z_{ad}^2$ and ${\mathbf{Q}}_0=\mathbf{0}\in {\mathcal{V}}_{ad}$ and run the optimization algorithm in the following cases:

  (a)

  Delta approximation (19) and linearized Equation (28): 252 iterations, optimal cost $1.2557307\times {10}^{-15}$, CPU time $180,895$ s.

  (b)

  Delta approximation (19) and adjoint state Equations (36a) and (36b): 263 iterations, optimal cost $1.1373912\times {10}^{-18}$, CPU time $90,977$ s.

  (c)

  Delta definition (18) and linearized Equation (28): 209 iterations, optimal cost $1.1029007\times {10}^{-12}$, CPU time $40,522$ s.

  (d)

  Delta definition (18) and adjoint state Equations (37a) and (37b): 211 iterations, optimal cost $4.4262276\times {10}^{-14}$, CPU time $14,785$ s.

  In all the above cases, the release rate and point used for the generation of artificial data are recovered by the algorithm.

4. Conclusions

In this paper, we propose a methodology for source identification of chemical incidents in urban areas. To achieve this, the classic advection-diffusion equation was completed with a reaction term depending on the air temperature, combined with a 3D microclimatic model, and numerically solved within the framework of mixed integer nonlinear programming (MINLP). In view of the results obtained, we observe that the proposed methodology is effective in identifying a source of contamination produced by a chemical agent in the academic cases studied. Both the proposed Dirac delta approximation and its own definition provide us with comparable results. The numerical resolution of the optimization problem using the adjoint state equations is more effective from the point of view of CPU time. Combining the search for the release point and rate requires a very high number of algorithm iterations, which penalizes CPU time.