Inverse Dispersion for an Unknown Number of Sources: Model Selection and Uncertainty Analysis

A simple recursive method is presented for performing the inverse dispersion modeling of an unknown number of localized sources, given a finite number of noisy concentration data acquired by an array of detectors. Bayesian probability theory is used to address the problem of selecting the source model which is most plausible in view of the given concentration dataset and all the available prior information. The recursive algorithm involves subtracting a predicted concentration signal arising from a source model consisting of N localized sources from the measured concentration data for increasing values ofN and examining the resulting residual data to determine if the residuals are consistent with the estimated noise level in the concentration data. The method is illustrated by application to a real concentration dataset obtained from an atmospheric dispersion experiment involving the simultaneous release of a tracer from four sources.


Introduction
Significant advances in sensing technology for concentration measurements of contaminants e.g., toxic industrial materials; chemical, biological, and radiological agents released into the atmosphere, either accidentally or deliberately, have fostered interest in exploiting this information for the reconstruction of the contaminant source distribution responsible for the observed concentration.Indeed, innovative methods for measuring and collecting concentration data in formal observation or monitoring networks of ever more compact nanotechnology-enabled sensors, including the adoption of new wireless sensor network technologies for acquisition of data in situ and for rapid delivery of this information through wireless transmission 1 , have placed the emphasis on data assimilation and inverse dispersion modeling.
where n is the number of parameters unknowns used to parameterize the source model, f O θ is an objective functional that measures the total mismatch usually taken as the sum-of-squared differences between the measured concentration data and the synthetic predicted concentration data associated with the current source model θ, and f R θ is a regularization functional that is used either to impose an additional constraint on the solution or to incorporate a priori information which is required to produce a mathematically unique solution .Most efforts in application of the optimization approach to inverse dispersion modeling have been focused on the development of numerical methods for solution of the optimization problem or on the prescription of the regularization functional for incorporation of various forms of a priori information.The optimization method for inverse dispersion modeling has been used by Robertson  While the optimization approach for inverse dispersion modeling seeks to provide a single optimal solution for the problem, the stochastic Bayesian approach seeks to generate multiple solutions to the problem with the evaluation of the degree of plausibility for each solution.In contrast to the optimization approach, the stochastic Bayesian approach allows the quantification of the uncertainty in the source reconstruction arising as such from the use of incomplete and noisy concentration data and imperfect models of atmospheric dispersion for the inverse dispersion modeling .This approach has been developed and used by a number of researchers e.g., 7-10 .In all the previous studies cited above, the inverse dispersion modeling involved the problem of identification of source parameters for a single localized source.The determination of the characteristics of multiple localized sources was briefly addressed by Yee 11 and by Sharan et al. 12 using, respectively, a stochastic Bayesian approach and an optimization least-squares approach that were similar to those applied previously for the recovery of a single source.In both these cases, this was possible because it was assumed in these investigations that the number of sources was known a priori.The problem of reconstruction of an unknown number of localized sources using a finite number of noisy concentration measurements is a significantly more difficult problem.A solution for this problem was proposed by Yee 13,14 who approached the problem as a generalized parameter estimation problem in which the number of localized sources, N, in the unknown source distribution was included explicitly in the parameter vector θ, in addition to the usual parameters that characterize each localized source e.g., location, emission rate, source-on time, source-off time .
Solving the reconstruction of an unknown number of localized sources as a generalized parameter estimation problem posed a number of conceptual and technical difficulties.The primary conceptual difficulty resided in the fact that when the number of sources is unknown a priori, the dimension of the parameter space or, equivalently, the dimension or length of the associated parameter vector θ is an unknown that needs to be estimated using the noisy concentration data.Yee 13,14 overcame this conceptual difficulty by using Bayesian probability theory to formulate the full joint probability density function PDF for the number of sources N and the parameters of the N localized sources and then demonstrated how a reversible-jump Markov chain Monte Carlo RJMCMC algorithm can be designed to draw samples of candidate source models from this joint PDF, allowing for changes in the dimensionality of the source model or equivalently, changes in the number of localized sources in the source distribution .Furthermore, it was found that the RJMCMC algorithm sampled the parameter space of the unknown source distribution rather inefficiently.To overcome this technical difficulty, Yee 13,14 showed how the RJMCMC algorithm can be combined either with a form of parallel tempering based on a Metropolis-coupled MCMC algorithm 13 or with a simpler and computationally more efficient simulated annealing scheme 14 to improve significantly the sampling efficiency or, mixing rate of the Markov chain in the variable dimension parameter space.
In this paper, we address the inverse dispersion modeling problem for the difficult case of multiple sources when the number of sources is unknown a priori as a model selection problem, rather than as a generalized parameter estimation problem as described by Yee 13,14 .The model selection problem here is formulated in the Bayesian framework which involves the evaluation of model selection statistics that gauges rigorously the tradeoff between the goodness-of-fit of the source model structure to the concentration data and the complexity of the source model structure.The model selection approach proposed herein for reconstruction of an unknown number of localized sources is conceptually and algorithmically simpler than that based on a generalized parameter estimation approach 13, 14 , leading as such to a simpler and more efficient computer implementation.More specifically, the model selection approach is simpler than the generalized parameter estimation approach in that it does not need to deal with the complexity of the variable dimensionality of the parameter space, which requires generally the development of sophisticated and computationally intensive algorithms.

Bayesian Model Selection
Let m N denote a source distribution model which consists of a known number, N, of localized sources.The model m N is characterized by a set of parameters θ.Bayesian inference for θ, when given the concentration data D, is based on Bayes' theorem rule : where I denotes the background contextual information that defines the multiple source reconstruction problem e.g., background meteorology, atmospheric dispersion model used to determine the source-receptor relationship, etc. and the vertical bar "|" denotes "conditional upon" or "given."Furthermore, in 2.
However, in our problem, the number of sources N is unknown a priori, so the relevant question that needs to be addressed is as follows: given a set of possible source models S ≡ {m 0 , m 1 , . . ., m N max } m N denotes the source model consisting of N localized sources , which source model is the most plausible probable given the concentration data D? This question is addressed rigorously in the Bayesian framework by using Bayes' theorem to compute the posterior probability for source model m N N 0, 1, . .Given any two source models m N 0 and m N 1 from the set S, the "odds" in favor of model m N 1 relative to model m N 0 are given by the ratio of the posterior probabilities of the two models, which on using 2.3 and 2.2 , reduces to Note that the posterior odds ratio K 10 is the product of the evidence ratio Z m N 1 /Z m N 0 and the prior odds ratio p m N 1 | I /p m N 0 | I .If every source model in the set S is equally probable namely, p m N | I 1/ N max 1 for N 0, 1, . . ., N max , then the prior odds ratio is unity, and the posterior odds ratio is identically equal to the evidence ratio.When viewed as a model selection problem, the inverse dispersion modeling of an unknown number of sources is conceptually simple: compute the evidence Z m N for each source model m N in S given the input concentration data D and the background information I, select the model with the largest value of the evidence as the most probable model, and recover the parameters θ corresponding to this most probable model.
The key to the model selection problem reduces to the computation of the evidence Z m N for an arbitrary value of N. A perusal of 2.2 shows that the evidence for the source model structure m N namely, a source model consisting of exactly N localized sources involves a computation of the overlap integral or inner product of the likelihood L θ and the prior distribution π θ over the entire parameter space for model m N .This inner product can be interpreted as a metric statistic for model selection and, as such, is an intrinsic element of the source model structure m N in the sense that it depends on both the set of parameters to be varied and the prior ranges for those parameters but is independent of the most probable preferred values for the parameters.
Note that the model selection metric selects preferentially a source model structure with the largest overlap of the likelihood and prior distribution, implying that in order for a more complex source model with a necessarily higher-dimensional parameter space to be selected, this requires the prior information to be smeared over a larger hypervolume in the parameter space in order to overlap significantly the likelihood.This mechanism embodies automatically an Occam's razor in that a simpler source model structure is preferentially selected unless the likelihood of the data for a more complex source model increases significantly more than what would be expected from simply fitting the noise residual between the measured concentration and the predicted concentration provided by the model .
The solution of the model selection problem is conceptually simple involving, as such, the computation of Z m N m N ∈ S .Unfortunately, this computation is a technically difficult problem for two reasons.Firstly, as N or, equivalently, as the number of localized sources increases, the overlap integral in 2.2 needs to be evaluated over a parameter space of increasing dimensionality n.More specifically as N increases, the dimensionality of the parameter space n increases linearly with N. Secondly, the integrand of the overlap integral of 2.2 becomes ever more complex with increasing N, owing to a counting degeneracy in the problem.In particular, it is noted that L θ is invariant under a reordering relabelling of the identifiers used to label the individual localized sources in the source model.This degeneracy simply corresponds to different but physically equivalent identifications of what is meant by a particular localized source in the source distribution, implying that the number of modes in L θ increases with m N as N!.The factor N! here corresponds simply to the number of possible permutations of the labels for the localized sources in the source model m N .In summary, the increase in the dimensionality of the parameter space and in the complexity of L θ with increasing N makes the problem of computation of Z m N technically difficult.
In view of the difficulty arising from the numerical evaluation of Z m N for increasing values of N, we consider an alternative but closely related methodology for addressing the model selection problem in the context of inverse dispersion modeling of an a priori unknown number of localized sources.To this purpose, the question that needs to be posed is after the removal of a predicted concentration for a source model m N from the measured concentration data, are the residuals that remain consistent with the estimated noise level or is there further evidence for the existence of an additional localized source in the residuals?This question is easily addressed using Bayes' theorem as follows.
Towards this purpose, the following model is assumed for the measured concentration data d i i 1, 2, . . ., N m : where C i is the true albeit unknown mean concentration, e is the measurement error, C i is the predicted mean concentration obtained from an atmospheric dispersion model, is the model error incurred by using the necessarily imperfect atmospheric model to predict C i , and e i is the composite error that includes both the measurement and model error.Furthermore, it is assumed that the expectation value e i of the composite error is zero, and the variance e 2 i is assumed for now to be known and given by σ 2 i .Hence, it is assumed that we are given the concentration data D ≡ {d 1 , d 2 , . . ., d N m } and corresponding uncertainties {σ 1 , σ 2 , . . ., σ N m } where σ i is the standard deviation square root of the variance for the i-th concentration datum d i .Normally, exact values for the uncertainties σ i are unknown so that all we would have available would be some estimated values for these uncertainties s i , rather than the true values σ i .
If we assume a priori that there are N localized sources in the domain, we can analyse the concentration data D given that m N is the correct model structure and draw samples of source models with exactly N localized sources encoded by θ from the posterior distribution given by 2.1 .Assume that N s samples of source models with N localized sources have been drawn from the posterior distribution; namely, we have available the samples θ k m N k 1, 2, . . ., N s where θ k m N is kth source model sample consisting of exactly N localized sources.To determine if N is the correct number of localized sources, we can compute the residual data E ≡ { 1 , 2 , . . ., N m } where i is the mean residual for the ith concentration datum with realizations of the residuals obtained by subtracting the predicted concentration C i m N for a sample of the source model m N encoded in θ k m N from the measured concentration d i .The mean residual datum i is estimated from the ensemble of source model samples θ k m N k 1, 2, . . ., N s as follows: Similarly, the uncertainty s i error in the residual datum i is estimated as s i , where

2.8
Now, given the residual data E, we evaluate the evidence for the presence of an additional source by considering two models hypotheses for E: 1 H 0 : there is no additional source and the residual data i C i m 0 e i e i where C i m 0 is the predicted concentration for a source model with no localized sources and hence has zero value , and 2 H 1 : there is an additional source and the residual data i C i m 1 e i where C i m 1 is the predicted concentration for a source model consisting of one localized source .
Under hypothesis H 0 null hypothesis , the residuals are consistent with the estimated noise level e i .To decide which hypothesis is favoured by the data, we calculate the posterior odds ratio for H 1 relative to H 0 .From 2.5 and assuming that hypotheses H 0 and H 1 are equally probable so p H 0 | I p H 1 | I 1/2, we get If K 10 exceeds some preassigned threshold, then there is evidence for an additional source namely, hypothesis H 1 is favoured with respect to hypothesis H 0 .If this is the case, then the number of sources N is increased by one, and we repeat the analysis of the concentration data D given that m N 1 is the correct model structure.This analysis would involve drawing samples from the posterior distribution of 2.1 with m N replaced by m N 1 , computing the residual data E obtained now by subtracting the predicted concentration C i m N 1 from the measured concentration d i , and then determining once more the evidence in favour of H 1 relative to that of H 0 .When this test fails namely, K 10 is less than or equal to the preassigned threshold value , the algorithm terminates.At this point, the number of sources N * has been determined, and the samples drawn from the posterior distribution of 2.1 for m N * can be used to estimate any posterior statistic of interest for the source parameters θ.
To summarize, the algorithm for the inverse dispersion modeling of an unknown number of sources consists of the following steps.The algorithm summarized previously determines the minimum number of localized sources needed to represent the concentration data D down to the estimated noise level, without the requirement to compute the evidence Z m N for N ≥ 2 a technically difficult problem computationally .Steps 2 and 3 in the algorithm correspond to the signal detection phase which addresses the specific question: is there evidence in the concentration data D for the existence of one or more localized sources in the domain?On the other hand, Steps 4 to 6 in the algorithm address the estimation phase and answer the question given that there is evidence for one or more localized sources in the domain detection has been confirmed , how many localized sources are present, and what are the values of the source parameters that best characterize each of these localized sources?

Source-Receptor Relationship (Dispersion Modeling)
In this paper, we are interested in the reconstruction of an unknown number N of localized sources.In view of this, the model for the source density function has the following explicit form for source model m N : 3.1 where δ • and U • are the Dirac delta and Heaviside unit step functions, respectively, and, Q j , x s,j , T b,j , and T e,j are the emission rate, vector position location , source-on time, and source-off time, respectively, of the jth source j 1, 2, . . ., N .For a source model or molecule composed of N localized sources or source atoms , it is convenient to define the parameter vector as follows: To apply the Bayesian probability theory, it is necessary to relate the hypotheses of interest about the unknown source model S to the modeled predicted concentration C cf. 2.6 .The predicted concentration C can be compared directly to the measured concentration d "seen" by a sensor at receptor location x r and time t r , averaged over the sensor volume and sampling time.The comparison can be effected by averaging the mean concentration C x, t over this sensor volume and sampling time to give where h x, t | x r , t r is the spatial-temporal filtering function for the sensor concentration measurement at x r , t r and R × 0, T corresponds to the space-time volume that contains the source distribution and the receptors sensors .Note that the mean concentration C x r , t r "seen" by a sensor can be expressed succinctly as the inner or scalar product C | h of the mean concentration C and the sensor response function h.
A source-oriented forward-time Lagrangian stochastic LS model can be used to predict C x, t .In this approach, C is estimated from the statistical characteristics of "marked" particle trajectories described by the following stochastic differential equation 15 : is the Lagrangian velocity of a "marked" fluid element at time t marked by the source as the fluid element passes through it at some earlier time t .The state of the fluid particle at any time t after its initial release from the source density function S is defined by X, V .In 3.3 , C 0 is the Kolmogorov "universal" constant that is associated with the Kolmogorov similarity hypothesis for the form of the second-order Lagrangian velocity structure function in the inertial subrange; is the mean dissipation rate of turbulence kinetic energy; dW t ≡ dW i t dW 1 t , dW 2 t , dW 3 t are the increments of a vector-valued three-dimensional Wiener process; a ≡ a i a 1 , a 2 , a 3 is the drift coefficient vector or, more precisely, the conditional mean acceleration vector .
Unfortunately, for the current application, the computational demands of the sourceoriented approach make it impractical as a direct tool for sampling from the posterior distribution because this involves necessarily exploring a large number of source parameter hypotheses.This is highly computer intensive, as the simulation-based Bayesian inference procedure requires a large number of forward calculations of the mean concentration to be performed, each of which involves the numerical solution of 3.2 and 3.3 .In this particular case, it may be useful to construct an emulator for the simulation model and use this emulator as an inexpensive surrogate to replace the forward model 16-18 .Fortunately, applying this type of approximation method for the forward model is not required for the current problem.It was shown by Keats et al. 7 and Yee et al. 10 that an exact computationally efficient procedure appropriate for use in a Bayesian inference scheme exists in the form of a receptor-oriented strategy for the representation of the source-receptor relationship.
Towards this purpose, the modeled sensor concentration C x r , t r can be computed alternatively using the following dual representation:

3.4
where C * x 0 , t 0 | x r , t r is the adjunct or dual concentration at space-time point x 0 , t 0 associated with the concentration datum at location x r at time t r with t 0 ≤ t r .A comparison of 3.2 and 3.4 implies the duality relationship C | h C * | S between C and C * through the source functions h and S.Moreover, C * is uniquely defined in the sense that it is explicitly constructed so that it verifies this duality relationship.
The adjunct field C * in the dual representation can be estimated from the statistical characteristics of "marked" particle trajectories corresponding to a receptor-oriented backward-time LS trajectory model, defined as the solution to the following stochastic differential equation with dt > 0 : where at any given time t , X b , V b is a point in the phase space along the backward trajectory of the "marked" fluid element here assumed to be marked or tagged as a fluid particle which at time t r passed through the spatial volume of the detector at location x r .Here, X b , V b defines the state of a fluid particle at any time t < t r before its "final" release from the sensor space-time volume h at time t r .It can be shown 15, 19 that C * obtained from 3.5 for a detector with the filtering function h and C obtained from 3.2 for a release from the source density S, are exactly consistent with the duality relationship V t , and 2 a b in 3.5 is related to a in 3.3 as where P E u is the background Eulerian velocity PDF.Thomson 15 provides one particular solution for the drift coefficient vector a that is consistent with the well-mixed criterion for the case where P E u has a Gaussian functional form namely, for Gaussian turbulence .
Finally, if we substitute 3.1 into 3.4 , the predicted concentration C θ ≡ C x r , t r ; θ "seen" by the sensor at x r , t r is given explicitly by the following expression: T b,j C * x s,j , t s | x r , t r dt s . 3.8

Prior
Assuming that each parameter in θ for source model m N is logically independent, the prior π θ can be factored as follows:

4.1
In this paper, the component prior distributions are assigned uniform distributions over an appropriate range for each source parameter.Furthermore, it is noted that the parameter ranges specified in the prior define the region in R n n ≡ 6N over which the evidence integral is carried out cf. 2.2 .More specifically, for j 1, 2, . . ., N, Q j ∼ U 0, Q max , where Q max is an a priori upper bound on the emission rate; x s,j ∼ U R where R ⊂ R 3 is the a priori specified spatial domain that is assumed to contain the source S; and, T b,j ∼ U 0, T max and T e,j ∼ U T b,j , T max , where T max is the upper bound on the time at which the source was turned on or turned off.Here, the symbol ∼ denotes "is distributed as" and U a, b is the uniform distribution defined on the closed interval a, b .Finally, note in the specification of the prior for T e,j that the lower bound of the parameter range is T b,j , encoding the fact that the time the source is turned off must necessarily occur after it has been turned on.

Likelihood
In the absence of a detailed knowledge of the noise distribution e i in 2.6 , other than the fact that it has a finite variance σ 2 i , the application of the principle of maximum entropy 20 informs us that a Gaussian distribution is the most conservative choice for the direct probability likelihood of the concentration data D. In consequence, the likelihood function L θ has the following form: where

4.3
In 4.2 , the notation for the likelihood function was adjusted to include the standard deviation for the noise ς ≡ σ 1 , σ 2 , . . ., σ N m in the conditioning to emphasize the fact that {σ i } N m i 1 are assumed to be known quantities.
Unfortunately, as mentioned previously, the noise term e i in 2.6 is rather complicated, arising as such from a superposition of both model and measurement errors.In consequence, reliable estimates for σ i i 1, 2, . . ., N m are difficult to obtain in practical applications.In view of this, let us denote by s i the quoted estimate of the standard deviation for the noise term e i for which the true but unknown standard deviation is σ i .Now, let us characterize the uncertainty in the specification of the standard deviation of e i with an inverse gamma distribution of the following form or, equivalently, prior distribution for σ i : where Γ x denotes the gamma function and α and β are scale and shape parameters, respectively, that define the inverse gamma distribution.The inverse gamma distribution for σ i has mean Again, the parameters α and β have been added to the PDF of the noise uncertainty in 4.4 to indicate that the values for these parameters are assumed to be known.
In view of 4.4 , the true but unknown standard deviations σ i of e i that appear in 4.2 can be treated as nuisance parameters and eliminated by considering the integrated likelihood The integrated likelihood of 4.6 can be interpreted simply as an average over all conditional likelihoods given ς, weighted by their prior probabilities.In so doing, the integrated likelihood incorporates the uncertainties regarding the standard deviations for e i i 1, 2, . . ., N m .The integrated likelihood function given by 4.6 depends explicitly on the hyperparameters α and β for which values need to be assigned.In this paper, the values of α and β are assigned as α π −1 and β 1.The assignment α π −1 ensures that σ i s i encoding our belief that our estimates s i of the standard deviation of e i are unbiased.Furthermore, the assignment β 1 results in a very heavy-tailed distribution for π σ i | s i , α, β which allows significant deviations of the noise uncertainty from the quoted value of s i provided by the user .Indeed, with the choice β 1, the variance associated with p σ i | s i , α, β in 4.4 becomes infinite.The heavy tail of the distribution is chosen to account for possibly significant under-estimations of the actual uncertainty namely, the quoted uncertainty s i σ i .This could arise from inconsistencies in the model concentration predictions owing to structural model error or to "outliers" in the measured concentration data owing to either measurement error or perhaps distortion of the measured concentration data due to some unrecognized spurious background source.
In order to compute Z m 1 required in the evaluation of K 10 in 2.9 , we are required to specify a functional form for the likelihood L θ m 1 | E, s of the estimated mean residual data E ≡ { 1 , 2 , . . ., N s } and the associated uncertainty in this data s ≡ {s 1 , s 2 , . . ., s N m }.To this end, the integrated likelihood function of 4.6 is also used for the specification of

4.7
Here, θ m 1 are the parameters of a source model consisting of a single localized source, used in the test to determine if the mean residual data E contains evidence for the presence of an additional source.Towards this objective, the likelihood function of 4.7 is used to calculate the evidence Z m 1 for this test as where n 6 corresponding to the dimension of the parameter space for one source .

Computational Methods
There are two major computational problems that need to be addressed in order to fully specify the algorithm summarized in Section 2, namely, 1 specification of a method to draw samples from the posterior distribution p θ | m N , D, I Steps 2 and 6 and 2 specification of a method for computation of the evidence Z m 1 Step 3, cf.The are a number of options that are available for the computation of the evidences Z m 1 and Z m 1 .Owing to the fact that the determination of Z m 1 and Z m 1 involves only an evaluation of a multidimensional integral in a low-dimensional space R 6 in this case , it is in principle possible to apply a brute force numerical integration 25 to address this problem.An alternative would be to apply a methodology that is applicable to the evaluation of the evidence in the general case, involving the evaluation of an overlap integral in a high-dimensional parameter space R n n 1 .These methodologies include importance sampling estimators 26 and thermodynamic integration 27 .
Although there are various alternatives leading to many different combinations that could be used potentially to draw samples from p θ | m N , D, I and to evaluate Z m 1 and Z m 1 , we have used instead a single methodology to do both.The methodology that is used in this paper to address both these problems is nested sampling developed by Skilling 28 for the efficient evaluation of the evidence Z for the general case.The nested sampling algorithm transforms the multidimensional evidence integral of 2.2 into a simple one-dimensional representation: where is the prior mass in the parameter hypothesis space enclosed with likelihood greater than L * and L χ is the inverse which labels the likelihood contour that encloses a prior mass χ.
If we evaluate the likelihood L χ at a sequence of m points , the evidence Z can be approximated from the likelihood-ordered samples as the following weighted sum: If the prior mass points χ i are sampled in a logarithmic manner as χ i i j 1 t j where t j ∈ 0, 1 is a shrinkage ratio, then the nested algorithm consists of the following steps.The reader is referred to Skilling 28 for further details of the algorithm.
Randomly draw M samples θ from the prior π θ to give an ensemble of samples.Evaluate the likelihood L θ for each of the samples in the ensemble.in Step 5 on termination of the algorithm.The key part of the algorithm is Step 3 involving drawing a sample from the prior π θ within a prescribed hard likelihood constraint L θ > L * .To this purpose, Feroz et al. 29 developed a very efficient algorithm which the authors refer to as MULTINEST for sampling from a prior within a hard likelihood constraint involving a very sophisticated procedure for decomposition of the support of the likelihood above a given bound L * into a set of overlapping ellipsoids and then sampling from the resulting ellipsoids.The algorithm is appropriate for sampling from posterior distributions with multiple modes and with pronounced curving degeneracies in a high-dimensional parameter space.
Finally, the posterior odds ratio K 10 in 2.9 is a summary of the evidence for H 1 additional source is present in the residuals against H 0 no source is present in the residuals .To interpret how strong or weak is the evidence for H 1 against H 0 provided by the residual data E, we use a reference scale suggested by Jeffreys 30 .In this scale, log K 10 < 1 corresponds to inconclusive evidence for H 1 against H 0 or the evidence for H 1 is not worth more than a bare mention and corresponds to a posterior odds ratio of less than about 3 to 1 in favor of H 1 .On this same scale, log K 10 > 5 corresponds to very strong evidence for H 1 against H 0 associated, as such, with a posterior odds ratio greater than about 150 to 1 .In view of this scale for interpreting K 10 , we choose the threshold K * for the posterior odds ratio required to be specified in Step 2 of the algorithm summarized in Section 2 to be log K * 1 or, equivalently, K * ≈ 2.7 .

Example: An Application of the Methodology
In this section, we apply the proposed methodology for inverse dispersion modeling of an unknown number of sources to a real dispersion data set obtained from a specific experiment conducted under the FUsing Sensor Information from Observing Networks FUSION Field Trial 2007 FFT-07 .The scientific objective of this field campaign was to acquire a comprehensive meteorological and dispersion dataset that can be used to validate methodologies developed for source reconstruction.Details of the instrumentation deployed and the experiments conducted in FFT-07 are given in Storwold 31 , so only a brief summary of FFT-07 will be presented here.In particular, only the relevant details of the particular experiment that are required for the interpretation of the results in this paper are emphasized.
The experiments in FFT-07 were carried out in September 2007 at Tower Grid on US Army Dugway Proving Ground, Utah about 2 km west of Camel Back Ridge on the Great Salt Lake Desert.The easterly through southerly drainage flows that predominate during the early morning hours at this site originate on the higher terrain to the southeast and are channeled by Camel Back Ridge.Generally, the terrain was flat, uniform, and homogeneous consisting primarily of short grass interspersed with low shrubs with a height between about 0.25 and 0.75 m.The momentum roughness length, z 0 , was estimated to be z 0 1.3 ± 0.2 cm.
In all the experiments, propylene C 3 H 6 was used as the tracer gas.The concentration detectors used were fast-response digital photoionization dPID detectors.In the experiments, a plume was formed in the atmospheric surface layer by releasing propylene from one or more up to a maximum of four purpose-designed gas dissemination systems.The network or array of concentration detectors consisted of up to 100 dPIDs, arranged in a staggered configuration of 10 rows of 10 detectors, with the rows of detectors separated by 50 m and the detectors along each row spaced 50 m apart.The concentration detectors along the ten sampling lines in the array were placed at a height, z d , of 2.0 m.
In the experiment used to test the inverse dispersion modeling methodology proposed herein, tracer gas was released continuously over a period of 10 min from four source locations at a height, z s , of 2.0 m: 1 source 1 is at x s , y s 33.0, 171.0 m; 2 source 2 is at x s , y s 33.8, 240.7 m; 3 source 3 is at x s , y s 30.0, 312.9 m; 4 source 4 is at x s , y s 26.0, 384.4 m.The coordinates reported here for the source locations are referenced with respect to a local Cartesian coordinate system.Unfortunately, in this experiment, only the mass flow controller for source 3 functioned properly.The mass flow controllers for sources 1, 2, and 4 failed to properly regulate the flow, so the emission rates from these sources were unknown.
A three-dimensional sonic anemometer, placed at the 2 m level of a lattice tower located upwind of the array of concentration detectors, was used to characterize the background micrometeorological state of the atmospheric surface layer.For this experiment, the horizontal mean wind speed S 2 at the 2 m level, the friction velocity u * , the atmospheric stability Obukhov length L , and the standard deviations of the velocity fluctuations in the alongwind σ u , crosswind σ v , and vertical σ w directions were S 2 3.61 m s −1 , u * 0.282 m s −1 , L −27.3 m, σ u /u * 2.33, σ v /u * 1.86, and σ w /u * 1.10.The mean wind direction in the experiment was normally incident on the detector array namely, the mean wind direction corresponded to a wind from the x-direction and, hence, is perpendicular to the sampling lines of detectors along the y-axis .
The wind velocity and turbulence statistics were used in conjunction with Monin-Obukhov similarity theory relationships 32 as input to the atmospheric dispersion model namely, the backward-time LS particle trajectory model described briefly in Section 3 for provision of the predicted concentration C cf. 3.4 .More specifically, the backward-time LS model of 3.5 was applied with the Kolmogorov constant C 0 4.8, a value which was recommended by Wilson et  max 0.05, 0.02d i ppm parts per million by volume where the lower limit of 0.05 ppm represents the precision in the concentration measurements using the dPIDs.The model error standard deviation is estimated to be 20% of the predicted value of the concentration C i θ namely, σ m,i ≡ s m,i 0.20 C i θ .In consequence, the estimated noise error standard deviation is given by σ i ≡ s i s 2 The algorithm described in Section 2 was used for the source reconstruction, applied to the concentration data D ≡ {d 1 , d 2 , . . ., d N m } obtained from 62 detectors in the array N m 62 see Figure 1 which shows the locations of the 62 detectors as filled squares .The hyperparameters defining the prior π θ have been chosen as follows: Q j ∼ U 0, Q max with Q max 100 g s −1 and x s,j , y s,j ∼ U R with R 0, 100 m × 0, 500 m for j 1, 2, . . ., N. After Steps 2 and 3 of the algorithm were completed, we found log K 10 log Z m 1 /Z m 0 64.0 ± 0.2 implying that the concentration data strongly support the source model m 1 one localized source or "signal" against the alternative source model m 0 no localized source or "no signal" .In other words, at this point in the algorithm, a source for the concentration has been detected in the data D. Owing to the fact that log K 10 > log K * 1, the algorithm continues to the iterative loop involving Steps 4 to 6.This iterative loop executes four times before it is terminated.The results of the execution of each of these four iterations are exhibited in Figure 1

Table 1:
The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parameters x s,j m , y s,j m , and Q j ≡ q s,j g s −1 for j 1, 2, 3, and 4 calculated from samples of source model m N with N 4 drawn from the posterior distribution p θ|m N , D, I cf. Figure 1 d 1. Comparing these estimated values of the source parameters, it can be seen that the proposed algorithm has adequately recovered the true parameters for each source if these are known to within the stated uncertainties.

Conclusions
In this paper, we have proposed a Bayesian inference approach for addressing the inverse dispersion of an unknown number of sources using model comparison selection .The necessary integrations to compute the evidence Z m N for N > 1 can be very computationally demanding as well as technically difficult , and, as a consequence, we developed an efficient and robust algorithm for model comparison that recursively removes the influence of a source model m N from the measured concentration data D and tests the resulting residual data E to determine if the residual data are consistent with the estimated noise level.This test requires nothing more than the computation of the evidence Z m N for N 1 which is usually computationally simple.The procedure finds the minimum number of localized sources necessary to represent the concentration signal in the data D down to the estimated noise level.Furthermore, the uncertainty in the estimated noise level which includes contributions from both measurement and model errors is treated by using an integrated likelihood, which was implemented using a specific prior inverse gamma distribution to represent the uncertainty in the estimated noise level.Nested sampling is used for the evidence computation, as well as for sampling from the posterior distribution p θ | m N , D, I in the proposed algorithm for inverse dispersion modeling.
The new algorithm has been applied successfully to a real concentration dataset obtained from an atmospheric dispersion experiment conducted in the FFT-07 field campaign consisting of the simultaneous release of a tracer from four sources.It is shown that the proposed algorithm performed well for this example: the number of sources was determined correctly N * 4 , and for each of the identified sources the corresponding parameters e.g., location, emission rate were estimated reliably along with the determination of the uncertainty in the parameter in the form of either a standard deviation or a credible interval that encloses a prespecified posterior probability mass .The methodology proposed herein offers a simpler alternative for the inverse dispersion modeling of an unknown number of sources that was addressed previously 13, 14 as a generalized parameter estimation problem, the latter of which involves necessarily the complexities in the design of an appropriate reversible-jump MCMC sampling algorithm to allow transdimensional jumps in the generalized parameter space namely, jumps between the parameter space of source models m N of different dimensions involving different numbers of sources N .

1 ,
p θ | m N , D, I ≡ p θ is the posterior probability distribution of the parameters θ, p θ | m N , I ≡ π θ is the prior probability distribution of the parameters θ before the data D was made available, and p D | m N , θ, I ≡ L θ is the likelihood function and defines the probabilistic model of how the data were generated.Finally, p D | m N , I is the evidence and, in the context of parameter estimation, ensures proper normalization of the posterior probability distribution assuming θ ∈ R n : . , N max to give p m N | D, I p m N | I p D | m N , I p D | I .2.3 Here, p m N | I is the prior probability of source model m N given only the information I, and p D | m N I is the global likelihood of the concentration data D given the source model m N and I and represents the goodness-of-fit of the model to the data.Note that the global likelihood of 2.3 is identical to the evidence Z m N defined in 2.2 .Finally, the marginal probability of the data given the background information, p D | I , is a normalization constant over all source models: p D | I N max N 0 p m N | I p D | m N , I .2.4 2.2 or Z m 1 Step 5, cf.4.8 .The most popular method that can be used for sampling from p θ | m N , D, I see 2.1 is a stochastic algorithm referred to as Markov chain Monte Carlo MCMC .There has been considerable effort expended to improve the efficiency and the chain mixing efficacy of MCMC sampling procedures in order to allow rapid chain convergence to and efficient sampling from the posterior distribution of θ.Improvements to MCMC sampling methods include: introduction of an adaptive Metropolis algorithm byHaario et  al.21 , formulation of the differential evolution Monte Carlo algorithm by Ter Braak 22 , development of a differential evolution adaptive Metropolis algorithm DREAM by Vrugt et al. 23 , and design of multiple-try Metropolis algorithms by Liu et al. 24 .

a N 1 ,
Figure 1 b N 2 , Figure 1 c N 3 , and Figure 1 d N 4 .The caption in Figure 1 summarizes the values for log K 10 log Z m 1 /Z m 0 obtained from testing the alternative hypothesis H 1 additional source present in the residual data against the "null" hypothesis H 0 no additional source present in the residual data .Each panel in Figure 1 also displays a density plot of samples of the source model m N drawn from the posterior distribution p θ | m N , D, I and projected onto the x, y horizontal plane namely, each point in the plot corresponds to the x, y location of a sample of the source model m N drawn from the posterior distribution .Note that the algorithm terminates with N * N 4 with log K 10 0.92 ≤ log K * 1, providing the correct number of localized sources in the source distribution N * 4 in this case .The samples of source models m N N 4 drawn from the posterior distribution p θ | m N , D, I in the last iteration of the algorithm see Figure1d can be used to determine the characteristics location, emission rate of the four localized sources.Towards this objective, Figure2displays the univariate diagonal and bivariate off-diagonal marginal posterior distributions for the parameters of each source.For each univariate marginal posterior distribution of a parameter for a given source, the solid vertical line indicates the true value of the parameter for this source , and the dashed vertical line corresponds to the best estimate of the parameter for this source obtained as the posterior mean.For the bivariate −500 −400 −300 −200

Figure 1 :
Figure 1: Four iterations of the algorithm of Section 2 showing the samples of source model m N drawn from the posterior distribution p θ | m N , D, I and the computation of K 10 for a N 1, log K 10 62.4 ± 0.2; b N 2, log K 10 57.3 ± 0.3; c N 3, log K 10 39.1 ± 0.3; d N 4, log K 10 0.92 ± 0.33.A solid square denotes a detector.

Figure 2 :
Figure 2: Univariate diagonal and bivariate off-diagonal marginal posterior distributions of the source parameters, namely, alongwind location x s , crosswind location y s , and emission rate q s that characterize sources 1 a , 2 b , 3 c , and 4 d obtained from samples of source model m N N 4 drawn from the posterior distribution p θ | m N , D, I cf. Figure 1 d .
, 2 , ..., N s } and the associated uncertainty in this data {s 1 , s 2 , ..., s N m } in accordance to 2.7 and 2.8 , respectively.5Using the information from the previous step, compute K 10 in accordance to 2.9 .If K 10 ≤ K 1 Input concentration data D {d 1 , d 2 , . . ., d N m } and associated estimated uncertainties {s 1 , s 2 , . . ., s N m }. 2 Set N 1 and specify a threshold K * for the posterior odds ratio.Draw N s samples of source models encoded by the parameter vectors θ k m N , k 1, 2, . . ., N s from the posterior distribution p θ | m N , D, I given by 2.1 .3 Compute the evidences Z m N and Z m N−1 , and use these values to determine K 10 in accordance to 2.5 .If K 10 ≤ K * , stop no sources are present in the domain, so N * 0 .4 Use samples θ k m N k 1, 2, . . ., N s to estimate the residual data E { 1 * , set N * N, output θ k m N k 1, 2, . . ., N s , and stop.6 Increase N → N 1.Draw N s samples θ k m N k 1, 2, . . ., N s of source models from p θ | m N , D, I given by 2.1 .Continue from Step 4.
5 where s ≡ s 1 , s 2 , . . ., s N m are the estimated standard deviations for the noise e 1 , e 2 , . . ., e N m and dς ≡ dσ 1 dσ 2 • • • dσ N m .Now, substituting the form for L θ | D, ς from 4.2 and 4.3 into 4.5 and performing the integration with respect to ς, we obtain the integrated likelihood function given by al. 33 from a calibration of the model against concentration data obtained from the Project Prairie Grass atmospheric dispersion experiments.The example used here to illustrate the inverse dispersion methodology involve continuously emitting sources, so the relevant source parameters for source model m N are θ ≡ θ m N x s,1 , Q 1 , . . ., x s,N , Q N .Furthermore, it is assumed that the height of the sources ISRN Applied Mathematics above ground level z s 2.0 m is known a priori, so the only unknown location parameters are the x s , y s coordinates of the sources.Owing to the horizontal homogeneity of the mean flow and turbulence statistics in the current example, the adjunct concentration C * x s | x d in 3.4 can be precalculated for one detector location x d for the known source height z s , with the adjunct concentration C * considered as a function of x s at all other detector locations obtained simply by a linear translation of C * x s | x d in the horizontal x, y -plane.In order to calculate the likelihood function given by 4.6 , we need to provide values for the estimated standard deviations s i for the noise e i .The noise error variance σ 2 i includes the sensor sampling error variance, σ 2 d,i , in the measurement of the concentration datum d i and the model error variance, σ 2 m,i , in the prediction of C i .The two contributions are combined in quadrature to give the noise error variance σ 2 .