Jilin Wang , Juanjuan Liu ,∗ , Bin Wang , Wei Cheng , Jiping Zhang

a State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, Institute of Atmospheric Physics, Chinese Academy of Sciences,Beijing, China

b College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China

c Beijing Institute of Applied Meteorology, Beijing, China

d College of Environmental Sciences and Engineering, Peking University, Beijing, China

Keywords:Pollution source identification Multiple sources Synthetic experiment Inverse modeling

ABSTRACT Rapid and accurate identification of the characteristics (source location, number, and intensity) of pollution sources is essential for emergency assessment of contamination events. Compared with single-point source identification, the reconstruction of multiple sources is more challenging. In this study, a two-step inversion method is proposed for multi-point pollution source reconstruction from limited measurements with the number of sources unknown. The applicability of the proposed method is validated with a set of synthetic experiments corresponding to one-, two-, and three-point pollution sources. The results show that the number and locations of pollution sources are retrieved exactly the same as prescribed, and the source intensities are estimated with negligible errors. The algorithm exhibits good performance in single- and multi-point pollution source identification, and its accuracy and efficiency of identification do not deteriorate with the increase in the number of sources. Some limitations of the algorithm, together with its capabilities, are also discussed in this paper.

1. Introduction

In the last decade, emissions from factories, sudden releases of hazardous contaminants or terrorist attacks have begun to attract increasing attention. Such incidents can pose a serious threat to human health or homeland security. Therefore, it is crucial for emergency assessment to rapidly and accurately identify the characteristics (source location,number, and intensity) of contaminant sources, which is a key factor in improving the forecasting skill for contaminant dispersion.

With increasing concerns about pollutant events, several methods have been proposed to identify emissions from point sources. Most of them have focused on the reconstruction of a single source ( Allen et al.,2007 ; Hutchinson et al., 2017 ; Keats et al., 2007 ; Ma et al., 2018 ;Matthes et al., 2005 ; Pudykiewicz, 1998 ; Qiu et al., 2018 ; Seibert, 2001 ;Singh et al., 2015a ; Zheng and Chen, 2011 ). However, releases from multi-point pollution sources are more common in real-world applications. For example, urban pollution may originate from different factory emissions, and toxic gases may be deliberately released by terrorists in different places. However, due to the overlapping concentrations of multiple sources, it is difficult to distinguish the contribution of each individual source and to accurately estimate its characteristics. In addition,different release scenarios may lead to similar concentration measurements. The identification of multiple sources is therefore assumed to be more challenging than that of a single source. Although limited attention has been paid to reconstructing multiple sources, some progress has nevertheless been made. Based on Bayesian theory and assuming that the number of sources is known, ( Yee, 2007 ) successfully estimated the characteristics (location and intensity) of multiple sources.Later, the Reversible Jump Markov Chain Monte Carlo algorithm was introduced to determine the number of sources ( Yee, 2008 ). Wade and Senocak (2013) proposed another method, incorporating Bayesian inference with a ranking system to reconstruct contaminants with an unknown number of sources. Bayesian methods have advantages in terms of providing uncertainties or confidence levels of source parameters.However, the applicability of this approach in an emergency response is limited due to the low efficiency. In addition, a linear least-squares approach was proposed by Lushi and Stockie (2010) to calculate source intensities, provided that the locations of multiple releases are known.Annunzio et al. (2012) used a Multi-Entity Field Approximation technique to determine the source locations of multiple contaminant releases. Later, a simple and weighted least-squares two-step approach was proposed by Sharan et al. (2012) and Singh et al. (2013) , respectively, to retrieve multiple point releases from limited concentration measurements. This algorithm does not require the initial guess of unknown parameters, but it assumes that the number of sources is known in advance.

It is worth noting that most multi-point pollution source identification methods have difficulty in determining the number of sources.Moreover, their computational costs rises significantly with an increasing number of sources. In view of that, a new algorithm free from any prior information, including the location and number of sources, is introduced in this study, which can automatically determine the number of sources and efficiently reconstruct multiple releases. Accordingly, the paper is organized as follows: The methodological theory is illustrated in Section 2 . Section 3 validates the applicability of the proposed method to single- and multi-point sources through a set of synthetic experiments.In Section 4 , the algorithm is summarized and discussed.

2. Methodology

In this subsection, a new two-step inversion method is introduced to reconstruct multi-point pollution sources. The new algorithm achieves the identification of source parameters in two steps: first, by identifying the number and locations of sources through an iterative process;and second, by minimizing the cost function to determine the respective source intensities. To improve computational efficiency, the source—receptor relationship (SRR) is introduced. As a vital concept in the field of inverse modeling, the SRR describes the potential sensitivity of the measurements with respect to the sources. Considering that the number of potential sources (namely, the number of grid points) is much larger than the number of receptors, the SRR is thus calculated from adjoint runs of the dispersion model. In this study, the SRR plays two major roles in the process of iteration. On the one hand, it serves as a linear operator converting the emission field into the simulated concentrations, which can be written as

whereHis them×nmatrix of SRRs derived from the dispersion model,in whichmis the number of measurements, andnis the number of computational grid points in the dispersion model.xis then× 1 vector representing the emission flux of sources, andysis them× 1 vector of simulated concentrations at the receptors. The utilization of SRRs (adjoint functions) avoids the forward run of the dispersion model for each computational grid point. Simulated concentrations can be directly obtained by multiplying the emission flux by the adjoint functions, which considerably reduces the computational cost. On the other hand, the SRR is weighted by observations to obtain the weighted influence function covering all possible source areas. Assuming that its spatial distribution is consistent with the weighted influence function, the field of emission is numerically proportional to that of the weighted influence function, and can be expressed in the form of

whereφ∗is then× 1 vector of the weighted influence function calculated from the SRRs and observations, andαis a constant corresponding to the scale factor between two fields. As with other inversion methods,we define the cost function as the squared difference between observed and simulated concentrations,

whereαbecomes the only unknown constant parameter in the cost function. It is evident that the minimization of Eq. (4) is much simpler than that of Eq. (3) . In this way, the estimation of the emission fieldxhas been transformed into the estimation of the constant parameterα, which significantly reduces the number of unknowns and further improves the computational efficiency. The optimal estimation of the parameterαcan be directly obtained with the least-squares method. Onceαis determined, the vector of emission fieldxcan be subsequently attained according to Eq. (2) .

The algorithm for identifying multiple sources is described in detail below. In the first iteration, the parameterα1is retrieved through the minimization of Eq. (4) , and the initial guess of the emission fieldx1 is accordingly obtained by the calculation of Eq. (2) . The residuals between observed and simulated concentrations related tox1are subsequently calculated in terms of Eq. (5) ,

Afterwards, it comes to the second iteration. It is known that if the initial guess of the emission field is not estimated accurately enough,there will be some information about sources remaining in the residuals. In the second iteration, the increment of the emission field (x2)is retrieved from the signal of residuals and utilized to further correct the initial guess (x1 ). The process of the second retrieval is exactly the same as that of the first one, except that the residuals are used to replace the observations when calculating the weighted influence function and minimizing Eq. (4). Once we get the emission increment, the estimated emission in the second iteration (S2) can be accordingly obtained by summing up the emission field of the initial guess and increment;namely,

Fig. 1. Maps of the relative locations of sources and receptors and simulation domains. (a) Sources are denoted by the black “×”symbols, and receptors by red triangles. (b) The black outer frames show the domains of the WRF simulation on a Lambert projection, and the inner red frame shows the domain of the FLEXPART model on a latitude—longitude grid.

whereHeis them×npmatrix of the SRRs at estimated locations, in whichnprepresents the number of identified point sources, andqis thenp× 1 vector of source intensities to be estimated. It is found that the form of Eq. (8) is extremely similar to Eq. (3) . However, the minimization of Eq. (8) is considered much simpler, since the number and locations of sources have already been identified and the dimensionality of unknowns is greatly reduced.

In short, the computational process can be summarized as follows:

Step 1. Generate the source—receptor relationship matrix from an atmospheric dispersion model and set the number of iterationp= 1.

Step 2. Calculate the weighted influence function with SRRs matrix and observations.

Step 3. Estimate the scale parameter by minimizing Eq. (4) .

Step 4. Calculate the initial guess or increment of the emission field with Eq. (2) .

Step 5. Calculate the estimated emission field and apply the nonlinear constraints of Eq. (7) to it.

Step 6. Calculate the residuals between observed and simulated concentrations.

Step 7. If the termination criterion is reached, go to step 8; otherwise,go to step 2 and replace the vector of observations with residuals.

Step 8. Determine the number and locations of sources and calculate the intensities of identified sources by minimizing Eq. (8) .

Step 9. Output the characteristics of identified sources, and the algorithm ends.

3. Synthetic experiments

3.1. Case description

In this section, the mathematical consistency of the source reconstruction algorithm is preliminarily examined with some synthetic experiments, involving the identification of one-, two-, and three-point sources. Ten cases are implemented for evaluation, where synthetic concentrations are generated from the simulation of dispersion model, without considering measurement errors or model uncertainties. The twostep inversion method described above is then applied to reconstruct the unknown source parameters with the obtained measurements. Thesecases mainly differ in their source locations and intensities. The respective source locations are shown in Fig. 1 (a). For the reconstruction of single-point sources, there are four cases corresponding to releases from S1, S2, S3, and S4, respectively. The source intensities are all set to 1.0 g s−1 . In addition, reconstructions of two- and three-point sources are carried out to demonstrate the practicability of the proposed method in multi-point source identification. Releases with equal (S2—S3, S3—S4)and unequal (S2—S3 ′ , S3—S4 ′ ) intensities at different locations are investigated in the two-point source cases. For cases S2—S3 and S3—S4, tracers are released with an intensity of 1.0 g s−1 at two source locations; while in cases of S2—S3 ′ and S3—S4 ′ , two point sources have different emission rates of 1.0 g s−1 and 0.5 g s−1 , respectively. Moreover, the remaining two cases (S1—S2—S3, S2—S3—S4) are related to three-point release scenarios, in which the intensities for all three sources are set to 1.0 g s−1.Detailed source term information is illustrated in Tables 1—3 .

Table 1 Reconstruction of one-point source cases.

Table 2 Reconstruction of two-point source cases.

For each case, the inert gas tracers are released continuously throughout the simulation period, from 1200 UTC 11 March to 0000 UTC 13 March 2011. Twenty receptors ( Fig. 1 (a)) are arranged in the downwind direction to sample tracers hourly from 0000 UTC 12 Marchto 0000 UTC 13 March 2011. A total of 480 synthetic measurements are thus obtained and can be utilized in the inversion. These receptors are distributed within a rectangular zone, with a spacing of 0.1°. The nearest distance between sources and receptors is approximately 20 km.

Table 3 Reconstruction of three-point source cases.

3.2. Numerical setup

The Lagrangian particle dispersion model FLEXPART-WRF( Brioude et al., 2013 ) is utilized to produce artificial observations and to retrieve the SRR matrix (adjoint functions) for source identification. It has been widely used to simulate the atmospheric dispersion over complex terrain and to identify unknown releases with limited measurements ( Seibert, 2001 ; Stohl et al., 2011 , 2012 ; Evangeliou et al.,2017 ; Tichý et al., 2017 ). For the simulation of meteorological fields,the WRF model is configured with two nested domains of 12 km and 4 km horizontal resolution ( Fig. 1 (b)). The physics options used are the Kain—Fritsch convection scheme ( Kain and Fritsch, 1993 ), WSM 3-class simple ice microphysics scheme ( Hong et al., 2004 ), Yonsei University planetary boundary layer scheme ( Hong et al., 2006 ), Noah land surface scheme ( Chen et al., 2001 ), RRTM-G (Rapid Radiative Transfer Model) longwave radiation scheme ( Mlawer et al., 1997 ), and Dudhia shortwave radiation scheme ( Dudhia, 1989 ). Both the initial and boundary conditions are provided by ERA-Interim ( Dee et al.,2011 ). Simulations are restarted every day from 11 to 13 March 2011,each of which lasts for 30 h. The first 6 h are discarded as spin up, and outputs of the last 24 h obtained from the inner domain (4 km) are finally combined to provide the meteorological input for dispersion.

For the backward run of the FLEXPART model, 500 000 particles are released at each receptor during the sampling intervals and followed backwards in time to calculate the SRR without considering any chemical or removal process. Its computational domain is illustrated in Fig. 1 (b). The outputs are sampled with 1 h intervals on a 351 ×251 grid with a 0.05°×0.05°resolution. As emissions are usually emitted at the ground level, we use the SRR values at a footprint layer 0—100 m above ground as the input of the inversion algorithm. As shown in Section 2 ,the SRRs will play important roles in the process of source reconstruction.

3.3. Results

In order to evaluate the retrieval performance, two indicators are introduced where the error in localization (EL ) is defined in the form of Euclidean distance between the estimated and actual source position,and the relative error of source intensity (EI) is formulated as

whereqe is the estimation andqtrepresents the true value of source intensity. For both indicators, a lower value corresponds to a higher accuracy of source reconstruction.

The retrieval results for one-point source identification are displayed in Table 1 , where source locations and intensities are expressed in grid units and g s−1 , respectively. It is shown that all four cases can be identified as single point sources, which is consistent with the actual situation. All sources are well identified with no location errors, and their intensities are accurately estimated with a maximum relative error of 0.00001%. In addition, the retrieval results for the first and last iteration are presented in Fig. 2 . The estimated locations of the top 10 highest emissions and their corresponding intensities are indicated. Note that the first iteration has already found the correct location of sources for most cases (except for the retrieval of source S4). However, the top 10 sources identified in the first iteration exhibit comparable intensities,making it difficult to determine the number of sources. Also the accuracy of convergence is not achieved; therefore, more iterations are required.An accurate estimation of the number and locations of sources can be gradually achieved with an increasing number of iterations. It is shown in Fig. 2 that the strongest emission is highlighted in the last iteration,and the source intensities at other grid points have been rapidly reduced to zero due to the nonlinear constraint of Eq. (7) .

Retrievals for two-point sources are displayed in Table 2 . It is observed that all four cases are well reconstructed by the new algorithm (Fig. S1). The number and locations are retrieved exactly the same as prescribed. Also, the maximum error for intensity estimation is 0.00004%, which is negligible. Checking the results in different cases, it can be found that the estimated errors with regards to identical (S2—S3 and S3—S4) and different (S2—S3 ′ and S3—S4 ′ ) intensities are comparable, which implies that the accuracy of retrieval is independent of the relative magnitude of source intensities. Likewise, the reconstruction of three-point sources also shows good performance ( Table 3 ). Locations are correctly identified, and intensities are estimated within an acceptable range with relative errors of no more than 0.00004% (Fig. S2).Comparing the estimated results of one-, two-, and three-point source cases, it is shown that the accuracy of reconstruction does not deteriorate with the increase in the number of sources. Moreover, the computational cost of the algorithm will not increase as the number of sources increases, which is a vital feature in multi-point source identification.For example, the average time for three-point source identification is about 47 min with a single CPU (Intel(R) Xeon(R) CPU E5-4640 v2 @2.20GHz machine). In fact, the efficiency of this algorithm can be further greatly improved by parallel computations.

4. Discussion and conclusions

In this paper, we have presented a two-step algorithm for single- and multi-point source identification. It first determines the number and locations of sources through an iterative process, and then minimizes a cost function to estimate the source intensities within a least-squares framework. The superiority of this algorithm lies in its ability to determine not only the source locations and intensities, but also to automatically identify the number of sources. Also, there is no need to provide prior information about unknown parameters.

The applicability of the proposed method is verified by a set of synthetic experiments. It is shown that the algorithm performs well in the reconstruction of one-, two-, and three-point sources, where the number of sources and respective parameters (source locations and intensities)of each source are retrieved exactly the same as predefined. Its accuracy and efficiency in source identification does not deteriorate with an increasing number of sources. Moreover, the ability of the proposed method to reconstruct five- and seven-point sources has also been verified with another group of experiments, the results of which are not reported in this paper for brevity.

It is worth noting that these synthetic tests conducted in this paper do not take any errors in account. Studies investigating the stability of the inversion algorithm with respect to measurement errors and meteorological field uncertainties, as well as the application of the algorithm to the real European Tracer Experiment, will be given in a separate subsequent paper. In addition, it can be seen from the algorithm that it can also be utilized to identify releases from area sources. Further attempts to invert the emission inventories of greenhouse gases could be made in the future. Also, we do not consider the release height of sources at present. Although ground-level approximation is commonly used in source term estimation, it was proved in Singh et al. (2015b) that neglecting the height of release can cause under-prediction in the retrieval of intensity. Further modifications or extensions of the algorithm are therefore required for more complex scenarios of contaminant dispersion.

Fig. 2. Retrieval results of the first iterations, (a) S1, (c) S2, (e) S3, and (g) S4, and last iterations, (b) S1, (d) S2, (f) S3, and (h) S4, for single-point source identification.The numbers 1—10 represent the top 10 points with high source intensities.

Funding

This study was supported by the National Key R&D Program of China[Grant Nos. 2017YFC1501803 and 2017YFC1502102].

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.aosl.2021.100098 .