Source Identiﬁcation of a Chemical Incident in an Urban Area

: This work deals aims to present a methodology for source identiﬁcation 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.


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: (i) 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. (ii) 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. (iii) 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 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: where r = 1 u 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, ∀x 1 , x 3 ≥ 0, x 2 ∈ R. (5) Hypothesis 9. The incident occurs at the initial time t = 0, and at a fixed point b = (0, 0, H), in such a way that the source term is given by S(x) = Q T δ b (x)δ 0 (t), where Q T is the total amount of chemical agent released.
Under these hypotheses, the general Equation (1) is rewritten as , Ω×]0, T[, and completed with the initial condition (5), and boundary conditions (2) formulated as ∀t > 0. The system is solved again by using the Laplace transform, and the explicit solution 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.

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 Ω 3D A ⊂ R 3 corresponding to an urban area, where ∂Ω 3D A 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 Γ I N A 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: where w C (m/s) is the sedimentation velocity (constant), u A (m/s) is the air velocity, θ 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 is the concentration of the chemical agent on the inlet boundary Γ I N A , and c 0 A (gr/m 3 ) is the initial concentration of the chemical agent. We must mention that ∇ · (u A c) = ∇ · u A c A + u A · ∇c A = u · ∇c A since we are assuming that ∇ · 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: where Q C (t) (gr/s) is the release rate and 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: 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: Version May 19, 2021 submitted to J 3 and completed with the initial condition (5), and boundary conditions (2) formulated 8t > 0. The system is solved again by using the Laplace transform, and the explicit solution is referred as the Gaussian-puff model. will be studied within the framework of optimal control problem of partial differential 67 equations (PDF). Thinking about a 3D urban domain like the one outlined in Figure   68 1, the classic advection-diffusion equation (1) will be completed with a reaction term 69 depending on the air temperature, and combined with a 3D micro-climatic model to    The temperature θ A (K) and air velocity u A (m/s) will be obtained by solving a microclimatic model in which we will take into account the temperature of the soil θ S (K) and buildings θ B [K] (see Appendix A).

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 (b C ) are in a bounded region U ad given by the union of M convex closed and bounded subsets (admissible release zones): where l k i and u k i are, respectively, the lower and upper bounds for the x i coordinate of Z k ad , i = 1, 2, 3. Let us empathize that the set U ad can be nonconnected. • Q C is the release rate or the total amount of chemical agent released. We will asume that the source term is given by (9), where X = L 2 (0, T) in the first case and X = R in the second one.
We consider the following objective function: We must observe that to evaluate the previous function at each element (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 Q C of the chemical agent is known: • Problem 2. Estimation of the release rate: We will estimate the release rate (or the total amount) Q C assuming that the release point b C is known: • Problem 3. Estimation of the release point and rate: 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 y C ∈ {0, 1} M such that y k C = 1, if the release point is in the zone Z k ad , and y k C = 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): where: • A ∈ M m 2 ×M with m 2 = 1 is such that Ay C = ∑ M k=1 y k C , therefore: and c = 1 ∈ R m 2 . Observe that the constraint Ay C = c implies that there is only one 1 in the control vector y C with the rest of its components being equal to 0. We will denote by Y ad = {y C ∈ {0, 1} M : Ay C = c} the admissible control set for the discrete variable y C .
If we fix the variable y * C ∈ Y ad we obtain a classical NLP problem: We denote by (b * C , Q * C ) a solution of the optimal control problem (14) associated with the discrete variable y * C . Thus, if we consider the following set: we can solve problem (13) by taking (b C , Q C , y C ) ∈ F (Y ad ) such that: For low values of M, the size of the set F(Y ad ) 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).

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 u A · ∇c A and we will use the notation u instead of u A .
Let us consider N + 1 points {t n } N n=0 in the interval [0, T] such that: We define α = 1 ∆t and we consider the material derivative for a scalar field c: We can consider the following approximation for the material derivative in a time t n+1 : being the solution of the following initial value problem: Thus, c n • X n h c n (x − u n (x)∆t). This approximation, as we will see later, is very important for obtaining the gradient of the objective function.
So, given c 0 A , we compute {c n+1 A } N−1 n=0 solving the following equation: For spatial discretization, we will assume that the domain Ω 3D A is polyhedral and we consider {τ A h } h>0 as a family of regular meshes of the domain Ω 3D A . We define the following finite element space: Thus, the fully discretized problem consists of {c n+1 A } N−1 n=0 ⊂ X A h solving the following variational formulation: 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: 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: From the previous function, we can define which converges to δ q (x) = δ(x − q) when → 0. Taking into account the previous sequence, If F C is given by (8), let us consider the following discretization for the control variable Q C (t): On the contrary, if F C is given by (9), we obtain 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, is the discretized objective function, where {c n A } N n=1 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 δ Q J(Q C , b C )(δQ) the directional derivative of J with respect to Q C in the direction δQ and by δ b J(Q C , b C )(δb) the directional derivative of J with respect to b in the direction δ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 Q using the linearized equations: where, given δ Q c 0 in case (19). In case (18) • Directional derivative of J with respect to b using the linearized equations: where, given δ b c 0 in case (19). In case (18) Remark 2. It should be noted that in Equation (25) the term Q n+1 C ∇z(b C ) · δ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 δ h (·, b) of the Dirac delta δ(x − b). Indeed, let us assume that b ∈ T ∈ τ A h and let φ( x) ∈ P 1 ( T) such that where T is the reference element and F T : T → T, with F T (x) ∈ [P 1 (T)] 3 . To obtain the above polynomial, let us consider a basis B = { p 1 , p 2 , . . . , p L } of the vector space P 1 ( T). We know that (26) is equivalent to: We denote by φ B = ( α 1 , . . . , α L ) t the coordinates of ϕ on the basis of B ( φ = ∑ L l=1 α l p l ). We know that φ B is the solution to the following linear system: where p B = ( p 1 , . . . , p L ) t and G B ∈ Sym(M L×L (R)) is the Gram matrix associated with B: Now, we define where J x F T (x) is the Jacobian matrix of F T . If we use expression (27): Taking into account the above definition, given an element z h ∈ W A h : Thus: 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, Therefore, to calculate, for example, ∂J ∂Q 1 (Q C , b C ), we have to solve (21) or (22) taking δ Q = (1, 0, . . . , 0), for computing ∂J ∂Q 2 (Q C , b C ) and we have to solve (21) or (22) taking δ Q = (0, 1, . . . , 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 Q and b using the adjoint state equation. On the one hand, where, given δc 0 in case (19). In case (18), δc n+1 Now, we will see how we can obtain the equations for the adjoint state. Let us consider {r n A } N n=0 ⊂ W A h such that r N A = 0 and let us take, for each n = 1, . . . , N − 1, r n A as a test function in (29) or (30): in case (19) and for (18): Taking into account that δc 0 in case (19) and for (18): Therefore, if we define r n A ∈ W A h , n = N − 1, . . . , 0, as the solution to: we know that: in case we take an approximation of the Dirac delta (19) and if we consider (18): 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 {r n in this case we consider (19); for (18):

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. The convective heat transfer coefficients for the corresponding interfaces were h A C = 100 To compute the radiation temperatures appearing in the heat equations for soil and buildings, we assumed that R sw,net (x, t) = (RM sw,dir + RM sw,di f f ) σ(x, t), and R lw,dow (x, t) = RM lw,dow σ(x, t), where, for our particular problem, we considered RM sw,dir = 650 Wm −2 , RM sw,di f f = 350 Wm −2 , and RM lw,dow = 450 Wm −2 . The function σ(x, t) ∈ [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 σ(x, t) = max{sin(2πt/86,400), 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 u 0 A = (0, 0) m s −1 , θ 0 A = θ 0 S = θ 0 B = 300 K, and the boundary conditions u I N A = (0, 10 −4 , 0) m s −1 , θ I N A = 300 K.
We also considered that G = 0, c 0 A = 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 ∆t = 900 s (N = 12 time steps). In Figure 4, we can see the air velocity u A and the air temperature θ A on the boundary ∂Ω 3D A 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.

Consideraremos un dominio acotado tridimensional ⌦ 3D
A ⇢ R 3 correspondiente a un núcleo urbano, siendo @⌦ 3D A la frontera de dicho dominio, esto es paredes y techos de edificios, suelo y, adicionalmente, fronteras ficticias que delimitan nuestro dominio. Dentro de estasúltimas, destacamos la frontera que llamaremos IN A por la que tenemos un flujo de aire entrante. Supondremos que la temperatura del aire puede afectar a la concentración del contaminante y que, eventualmente, puede haber efectos de sedimentación. Por lo tanto, la evolución de la concentración del agente químico c A vendrá dada por la solución de la siguiente ecuación en derivadas parciales: 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 b R = (2, 12, 6) ∈ Z 2 ad = [1.  (17), associated with the release point b R and the release rate Q R (t) = 3 (1 + sin(2 π t /86,400 − 2 π ∆t /86,400)) gr m −3 s −1 , t ∈ [0, T], such that Q R (t) ∈ 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.

Measurement points
Therefore, the main objective of this section is to show how the methodology we propose allows one to recover the release point b R and the discharge rate 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 :

1.
Problem 1, Equation (10). Estimation of the release point: In this case, we fix the release rate to Q R and we try to recover the release point b R using the artificial measurements. To this end, we start from the initial control b 0 = [1.2, 10.2, 4.2] ∈ Z 2 ad and run the optimization algorithm in the following cases: (a) Delta approximation (19) and linearized Equation (23) Problem 2, Equation (11). Estimation of the release rate: In this case, we fix the release point to b R and we try to recover the release rate Q R using the artificial mea-surements. We start from the initial control Q 0 = 0 ∈ V ad and run the optimization algorithm in the following cases: (a) Delta approximation (19) and linearized Equation (20) Problem 3, Equation (12). Estimation of the release point and rate. We try to recover the release point b R and the release rate Q R using the artificial measurements. We start from the initial control b 0 = [1.2, 10.2, 4.2] ∈ Z 2 ad and Q 0 = 0 ∈ V ad and run the optimization algorithm in the following cases: Delta approximation (19) and linearized Equation (28)

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.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results'.

Appendix A
In this appendix, we include the microclimate mathematical model that we have used in the numerical simulations. The microclimate model is similar to that used by authors in [14]. We summarize it here for convenience of the reader. So, we consider a 2D domain Ω 2D = {(x, y) ∈ R 2 : 0 < x < a, 0 < y < b} and two positive functions H S , H A : Ω 2D → [0, ∞) that represent the height of the layers of soil and air. We also consider a subdomain Ω 2D B ⊂ Ω 2D corresponding to the buildings and a function H B : Ω 2D B → [0, ∞) representing the height of these buildings. Then, we define the following 3D domains: which correspond, respectively, to the domain occupied by soil, buildings and the air. We can see that ∂Ω 3D We consider the following equations that model the behavior of the air temperature θ A (K), density ρ A (m 2 s −2 ) and velocity u A (m s −1 ): where ν A is the kinematic viscosity coefficient, θ REF A is a reference temperature, g is the gravity acceleration, n A is the unit outward normal vector to the boundary ∂Ω 3D A , and u I N A , u OUT A and u 0 A are given boundary and initial conditions. Observe that we are considering a directional do-nothing condition for the Navier-Stokes equations [26] where are the convection coefficients, F A is a heat source term, and θ I N A and θ 0 A are given boundary and initial conditions. The soil temperature θ S (K): where K S is the diffusion coefficient, b A,S 1 and b B,S 1 are the convection coefficients, b A,S 2 is the radiation coefficient, T A,S r is the radiation temperature induced by solar radiation, F S is a source term, and θ SUB S and θ 0 S are given boundary and initial conditions. The buildings temperature θ S (K): where K B is the diffusion coefficient, b A,R 1 , b A,W are the radiation temperatures, F B is a source term, and θ 0 B is a given initial condition. The characteristic parameters that define the thermal behavior of the materials involved in the problem are the following: • ρ A , ρ S and ρ B (g m −3 ) are the densities of air, soil and buildings, respectively; • cp A , cp S and cp B (Ws g −1 K −1 ) are the specific heat capacities of air, soil and buildings; • α A , α S and α B (Ws g −1 K −1 ) are the thermal conductivities of air, soil and buildings; • S , W and R (dimensionless constants) are the emissivities of the surfaces corresponding to soil, walls and roofs, respectivel; • a S , a W and a R (dimensionless constants) are the albedos of soil, walls and roofs, representing the ratio of reflected radiation from the surface to incident radiation upon it; • h A S , h B S , h A W and h A R (W m −2 K −1 ) are the convective heat transfer coefficients between soil/air, soil/buildings, walls/air and roofs/air, respectively.
From the above coefficients, we can define the coefficients associated with Equations (A1)-(A4): • K A , K S and K B (m 2 s −1 ) are the thermal diffusivities of air, soil and buildings, defined from the above data in the following way: and b A,R 1 (m s −1 ) are the coefficients related to convective heat transfer, obtained from the following relations: for the temperature of air: for the temperature of soil: for the temperature of buildings: • b A,S 2 , b A,W 2 and b A,R 2 (m s −1 K −3 ) are the coefficients related to radiative heat transfer for soil, walls and roofs, respectively, obtained from following relations: with σ B (W m −2 K −4 ) the Stefan-Boltzmann constant. • Finally, in order to compute the radiation temperatures T A,S r , T A,W r and T A,R r (K) on the different solid boundaries (soil, walls and roofs), we use the following expressions (involving the corresponding solar radiations, albedos and emissivities): where R sw,net (x, t) denotes the net incident shortwave radiation on the surface, and R lw,dow (x, t) denotes the downwelling longwave radiation, both measured in W m −2 .