MODEL REDUCTION AND POLLUTION SOURCE IDENTIFICATION FROM REMOTE SENSING DATA

. We consider a source identiﬁcation problem related to determina- tion of contaminant source parameters in lake environments using remote sensing measurements. We carry out a numerical example case study in which we employ the statistical inversion approach for the determination of the source parameters. In the simulation study a pipeline breaks on the bottom of a lake and only low-resolution remote sensing measurements are available. We also describe how model uncertainties and especially errors that are related to model reduction are taken into account in the overall statistical model of the system. The results indicate that the estimates may be heavily misleading if the statistics of the model errors are not taken into account.


1.
Introduction. Inverse problems are typically encountered in situation in which the objective is to carry out inference on factors affecting the behavior of a system by making observations on the system state. Such factors of interest are often considered as the parameters of the underlying mathematical model such as material properties and object shapes and structures. Parameter identification problems arise in several medical and industrial imaging applications such as optical tomography [2], electrical impedance tomography [11], acoustic imaging [13], seismic exploration [41], magnetic induction tomography [18] and X-ray tomography [34].
An important subclass of model parameter estimation problems is known as source identification problems, or inverse source problems, in which the aim is to find estimates to model parameters that are understood as source terms related to different phenomena. The determination of source parameters plays an important role in many fields of science and engineering, and such problems are encountered e.g. in neurophysiology [7], seismology [30] and acoustics [33].
The identification of sources related to advective-diffusive mass transport is an increasingly important issue in occupational and public safety management and in environmental protection, see e.g. [39]. In a case of accidental or intentional release of a hazardous substance it is essential to obtain information on the source in order to assess the damages and conduct protection operations. For contaminant source identification it would of course be ideal to have an extensive monitoring network that would give information on the contaminant distribution itself. Combining these (direct) measurements with an evolution model such as a advection-diffusion model for the contaminant, would then facilitate the estimation of the source properties. Such measurement networks, however, rarely exist and indirect measurements can typically only be accessed. If the measurement data is sparse, has relatively low resolution or is partially missing, and/or the model between the measurements, the contaminant itself is uninformative and the concentration distribution is timedependent, we are dealing with a nonstationary inversion problem [23].
Several approaches have been proposed to solve contaminant source identification problems in different situations such as in groundwater, atmospheric and lake environments. Backward propagation methods can be used for reconstructing the history of a contaminant plume originating from instantaneous sources. The objective is to determine an estimate for the initial concentration distribution and thus determine the source location and the released amount of substance. The methods that have been considered include e.g. the quasi-reversibility method [45], the marching jury backward beam equation method and its variants [4,5], adjoint methods [28] and particle methods [6]. A problem of backward propagation approaches is that the time instant of the release is usually not known.
Another class of methods use forward propagation of governing models, and the basic principle is to determine source parameters so that an optimal fit with measured data is obtained. Such approaches can be understood as optimization approaches in which the solution is sought subject to the governing transport model. With simplified models and geometries the solution of the transport model can be expressed analytically. Analytic solution for an instantaneous point source has been used e.g. in [25]. In [24] an analytic solution was used in the estimation of the release history of a point source with a Tikhonov regularization approach. Furthermore, a correlation coefficient optimization approach relying on analytic solutions has been presented in [35]. Analytic solutions have also been exploited in the estimation of time-varying boundary concentrations that were considered as contaminant sources. Due to the ill-posed nature of the problem, Tikhonov regularization and minimum mean entropy approaches have been used e.g. in [43,37,42]. The quasi-reversibility approach has also been used to solve a similar problem [36].
Since analytic solutions exist only in simple situations that are not feasible in real world cases, the solution of the contaminant transport model has to be implemented numerically e.g. with finite element, volume or difference methods, or boundary element methods. Also, in such cases it is often necessary to use numerical methods to solve also the flow field. Optimization approaches with numerical schemes have been used with various source models such as point sources [31,32] and distributed sources [10,1] in steady state situations, as well as point sources at known locations [17] and distributed sources [1] in non-steady state cases. In addition to traditional optimization approaches, a Bayesian approach with Markov chain Monte Carlo (MCMC) estimation was proposed in [12]. MCMC methods have been used also in [40] for the estimation of initial concentration distribution. Naturally, the applicability of any approach depends also heavily on the observation modality, and typically the observations consist of point measurements with varying spatial coverage.
If we are dealing with an (nonstationary) inverse problem, it is evident that the accuracy of mathematical models bears a significant impact on the accuracy and feasibility of the estimates. With feasibility we refer to the fact that the estimates as well as the error estimates should be consistent with reality. When modeling errors and uncertainties exist, neglecting these in the overall model will not only result in poor estimates but also in error estimates that may be significantly overoptimistic. Typical modeling errors and uncertainties include approximate geometry, (partially) unknown coefficients (secondary variables), (partially) unknown boundary and initial data, inaccuracies due to numerical implementation and the physical model itself may be only a more or less crude approximation for the situation at hand. Until recently, model errors and their effects on the solutions of inverse problems have not been considered systematically. Recently, a theory for analyzing model errors and, in addition, for compensating their effect in inversion has been developed [23,22]. This theory is based on Bayesian statistics and on the modeling of all uncertainties as random variables (or operators). The theory has been successful with several applications in which it is either advantageous or practically unavoidable to use reduced or approximate models in inversion. By analyzing and quantifying the joint statistics of the uncertainties that are related to a numerical model and the other unknowns, it may be possible to use significantly reduced models with only small effects on the estimates, see e.g. [3,29,27]. As a result, the error estimates using the reduced model are slightly increased but, most importantly, these are consistent with the reality. This theory has been developed also for nonstationary inverse problems [20,19,21].
Model errors have not been considered extensively in the estimation of contaminant source properties. In [17] a reduced model was used, but the effects of model errors were ignored. Furthermore, the effects of different discretizations were tested in [10], and it was observed that discretization can be decreased to some extent without significant effect on the estimates but at some point the effects of discretization errors start to increase. However, the limit at which the discretization errors start to become more and more significant is usually hard to predict.
In this paper we consider the problem of detecting and estimating contaminant sources in lake environment based on remote sensing observations from a timevarying concentration distribution. We focus on the model reduction and its effects and take the boundary conditions and other secondary parameters to be (approximately) known. We present a simulated example case in which a reduced model obtained by using coarse discretization is employed in data processing. The focus is on model error analysis, and we use the approximation error approach in order to account for the model errors in data processing. The example shows that the use of a reduced model without a statistical model for the induced errors can lead to misleading (infeasible) estimates.
2. Models and model errors in source identification. In general, not restricting to any specific environment and application, the objective in contaminant source identification is to make inferences on the properties of pollution sources by making observation on the contaminant plume and utilizing a transport model. Contaminant transport is most commonly modeled by the advection-diffusion equation is the volume source term, and x is the spatial coordinate. The aim in source identification is to determine an estimate for the source term s(x, t).
The model (1) equipped with initial and boundary conditions defines the (deterministic) evolution model. If the initial and boundary conditions together with the model parameters v(x, t) and D(x) are known and the source term is the only unknown, the task is obviously a radically simpler one when compared to the case in which these are (at least partially) unknown. Note also that while it is possible to simulate the flow field (for example in a lake), the field is time-dependent and the actual flow field at the time of the measurements is unknown.
Irrespective of the environment in question, fluid flow can be described by the continuity equation and the Navier-Stokes equations (see e.g. [9]) where ρ = ρ(x, t) is the fluid density, p = p(x, t) is pressure, µ = µ(x, t) is the fluid viscosity, the coefficient λ is related to bulk viscosity and usually set as λ = − 2 3 µ, g is the gravitational acceleration and ε(v) is the strain tensor of the form Equations (2) and (3) are valid for compressible Newtonian fluids, i.e. for gases. In a case of incompressible fluids (such as water) the Navier-Stokes equations take a slightly simpler form [9]. With certain assumptions it is possible to derive approximate flow models from the general Navier-Stokes equations. For instance, macroscopic flow in saturated porous media can be described with the Darcy's law, which is used in groundwater flow models. Flow in unsaturated porous media can be modeled with the Richard's equation which is based on the Darcy's law [8].
It is obvious that model uncertainty issues are present also in the implementation of the above-mentioned models as discussed in Section 1. In addition, further uncertainties are also due to the temperature dependency of flow and diffusion phenomena. Therefore the heat balance model with appropriate boundary conditions and flow and advection-diffusion models should be coupled and solved simultaneously. Furthermore, in the case of a reactive contaminant, the entire transport model includes additional transport and reaction models for all substances reacting with the primary substance of interest. All in all, the implementation very accurate models can be a complicated task and lead to heavy numerical models that cannot be used for practical applications. In such cases the approximation error theory may enable the use of computationally feasible models without significant effects on the results.

Problem setup and numerical models.
3.1. Description of the problem. We consider a simulated example case in which the objective is to determine the properties of a contaminant source. We consider a fictitious lake (domain Ω), see Figure 1, in which the flow is driven by the input flow (pressure difference) and wind traction. At the bottom of the lake there is a pipeline. At certain time instant a leakage appears in the pipeline, and the substance starts spreading in the lake. The substance is optically active in the sense that it changes the reflectance of water. The lake is monitored with a low resolution spaceborne or airborne remote sensing system providing spectral reflectance data on the lake surface, and we have a model that connects the reflectances and surface concentrations. The location, release rate and the time instant at which the release is started are the quantities to be estimated on the basis of the remote sensing data.

Generation of simulated data.
3.2.1. Flow model. In order to simulate measurement data we fix the initial and boundary data and other parameters. The flow is considered time-invariant and incompressible, and the viscosity of water is independent from temperature. It is assumed that the Reynold's number is not very high and therefore there are no large-scale eddies in the lake. Thus the flow is approximated with the steady-state Navier-Stokes model The boundary conditions are set as follows: a) v = 0 at the lake bottom, b) input flow is parabolic-like and perpendicular to ∂Ω input , c) output flow is perpendicular to ∂Ω output , d) 2µ ε 1 (v) · n − pn 1 = 0 at the output boundary ∂Ω output (∂Ω output is perpendicular to x 1 axis and n is outward unit normal vector), e) v 3 = 0 on the lake surface and f) 2µ ε m (v) · n − pn m = T m (x 1 , x 2 ), m = 1, 2, on the lake surface.
Be defining an appropriate traction function T(x 1 , x 2 ), the velocity field can be solved with a finite element method by employing streamline upwinding Petrov-Galerkin (SUPG) stabilization and Picard's linearization scheme [14]. The mesh in which the data is computed consists of 11,521 nodes and 52,824 tetrahedra and is later denoted by M 0 . The simulated traction and the resulting velocity field are illustrated in Figure 2 3.2.2. Advention-diffusion model. Contaminant transport is described with model (1) assuming that the diffusion coefficient is constant and that the substance of interest is water soluble with density equal to that of water. We assume that the contaminant substance is not naturally present in lake environment and that there are no potential sources in upper water ecosystems. Thus, we set c(x, t = 0) = 0 and c(x, t) = 0 at ∂Ω input . In addition, we assume that diffusion through the lake surface and bottom is negligible, and that diffusive mass transport through the output flow boundary is very small in comparison with advective mass transport. Therefore we can write ∇c · n = 0 at lake bottom, surface and output flow boundaries. The leakage from the pipeline is modeled as a volume source so that the source function s(x, t) is spatially homogeneous in the source element Ω source while s(x, t) = 0 in Ω \ Ω source . The source term s(x, t) is mathematically defined as where γ(t) is the release rate, V Ωsource is the volume of the source element and χ Ωsource (x) is the characteristic function of the source element. The location of the source element Ω source and the release rate in arbitrary units are shown in Figure  3.
Numerical solution is sought in mesh M 0 with a SUPG stabilized finite element (FE) approach [15] using piecewise linear basis functions. This leads to a set of ordinary differential equations where the vectorc contains the concentrations at nodal points in Ω\∂Ω input . Matrices A and B depend on discretization and on selected basis functions, and they are independent of time. Vector S is related to the volume source. The set of ordinary differential equations is solved numerically using the backward Euler method which leads to the update scheme where the subscript ℓ refers to time instant t = t ℓ and ∆t = t ℓ+1 − t ℓ . Matrix B is constructed using diffusion coefficient D = 1 × 10 −6 , and the time step is ∆t = 400 (both in arbitrary units). Since the flow is time-invariant, the diffusion coefficient is selected so that it accounts for mass transport due to small-scale eddies.
3.2.3. Observation model. In optical remote sensing of substances in water environment the relation between the concentration and observations is typically modeled using semiempirical optical models or empirical regression models, see e.g. [26,38]. For data generation we employ a typical slightly nonlinear regression model R λq , and R λp and R λq are observed reflectances at two different wavelength bands characterized by mean wavelengths λ p and λ q , respectively,c is the mean surface concentration of the substance and ν is observation noise due to instrument noise and varying atmospheric conditions.
The lake surface is divided into n pixels = 39 non-overlapping pixels and observations representing different time instants are generated with equation (10) using parameters α 1 = −0.00110, α 2 = 0.128 and α 3 = 0.450. Observation noise is Gaussian with zero mean and standard deviation σ = 0.18. In addition, noise  3.3. Models in nonstationary inversion. The computational models in the inversion are different from those with which the data are generated. The coefficients and boundary data are set only slightly off their actual values. The main variable here is, however, the accuracy of the numerical approximation, that is, the density of the mesh in the finite element model. This choice is reasonable since the most commonly used technique to reduce the computational cost is to decrease the density of discretization.
3.3.1. Flow model. Since exact modeling of the flow field in a real lake is very difficult or impossible in practice, it is unrealistic to use the same model in inversion as in the data generation phase. Thus we use a slightly different model in which the governing equation is the same as earlier but flat flow profile (constant flow speed in all nodal points but zero speed on the lake bottom) is used at the input flow boundary. Also the flow speed is selected so that the resulting mean flow speed at the input flow boundary is 1.77% smaller than the actual mean flow speed. In addition, we use slightly modified surface traction in order to account for the uncertainties in wind conditions. When evaluated at the same points as in Fig. 2  (upper figure) the relative error in the magnitude of the wind traction vary from 0% to 13.4% (mean relative error 1.79%), and the absolute error of angle is in the range of 0.0 • -4.6 • (mean deviation 0.8 • ). Mesh M 2 is used in inverse models and therefore the discretization is significantly coarser than in data generation, see Table 1. Similar SUPG stabilization method is used as in Section 3.2.1.

3.3.2.
Mass transport and source models. The evolution of the contaminant plume is described with the advection-diffusion model (1). Boundary and initial conditions are the same as in data generation but the value of the diffusion coefficient is set 1% smaller than the actual value. The FE implementation and backward Euler scheme lead to the following update equation (8) and (9). Matrices A and B are now related to the mesh M 2 . In addition, in this model we use coarser discretization also with respect to time, and the time step is ∆t = 500. See [20,19,21] for an exact formulation of the sparse mesh -long time step related approximation errors.
In principle, the evolution model (11) could be used in the estimation of source properties. However, the time step in equation (11) is too small for practical purposes leading to a large number of matrix-vector multiplications and significant computational burden. In order to avoid this, we construct a multistep update scheme in which one-step update equations are combined to obtain an update equation over a longer period of time. The equations of the lth multistep leap (c N (l−1) →c N l ) consisting of N one-step updates can be written as +2 + · · · + F r N l−1 + r N l Assuming that the source terms associated to the lth multistep update do not vary with time, i.e. r N (l−1)+1 = r N (l−1)+1 = . . . = r N l ≡ ρ l , a general multistep update equation can be written with the time indexing of inverse computations as (12)c τ +1 = Kc τ + q τ +1 , τ = 0, 1, . . . where K = F N and q τ +1 = (F N −1 + F N −2 + . . . + F 2 + F + I)ρ τ +1 . The necessary number of steps corresponding to the data acquisition rate is N = 80. By including the concentrationsc τ and the concentrations at the Dirichlet boundary to the same vector c τ , the deterministic part of the evolution model can be written in a general form as (13) c τ +1 = f τ (c τ ).
Since the initial and boundary conditions are fixed and the concentration distributions depend only on source parameters, we can write where θ is the source parameter vector. The source model is constructed by assuming that there are no sources of the substance other than the one from the pipeline. Since the location of the leakage is not known, we select N k = 31 elements from the bottom of the lake to describe potential leakage locations, see Figure 5. In addition, it is assumed that the leakage comes from a single element (index k) and it starts at certain time instant (τ * ) and continues with constant release rate (γ). Thus the source parameter vector is θ = (τ * , γ, k) T .
Using the above assumptions on the contaminant source, the source related term ρ τ , and further q τ , are computed using the model where (15) δ τ * (τ ) = 1, when τ ≥ τ * 0, otherwise , the constant γ is the release rate (mass per unit time) andŜ (k) corresponds to unit release rate from the k'th candidate element for the source location. More specifically,Ŝ (k) is constructed by defining the source term s in equation (1) as where V (k) is the volume of the kth source candidate element and χ (k) (x) is characteristic function of the same element.

3.3.3.
Observation model. We assume that calibration measurements relating the remote sensing observations and mean concentrations are available. The model (10) is used to generate the calibration data which is shown in Figure 6. We assume that we have remote sensing data available for time instants τ = 0, 1, . . . , τ max . In the simulation study the final time instant is either τ max = 28 or τ max = 24.
A linear regression model is used to describe the relation between mean concentration and observation in the ith pixel. The model estimated on the basis of data shown in Figure 6 is of the form  Figure 6. Calibration data, the linear regression estimate (bold line) and the nonlinear regression model (see equation (10)) used in the data generation (thin line).
with zero mean and mutually independent with concentrations and noise terms in other pixels. The actual errors are spatially correlated.
The relation between concentration distributions and observations at all time instants can be generally written as

Approximation error analysis and estimation of source parameters.
In optimization approaches targeting at minimizing the (regularized) output error, the possibly nontrivial statistics of the measurement errors are sometimes (but not always) taken into account. However, the discrepancy between the model and reality induces additional errors and if their effects are ignored, the estimates for the source parameters may be inaccurate with respect to the associated error estimates, see [23,22]. In the following we review briefly the approximation error approach and explain how the approximation error statistics were computed in this example.

4.1.
Computational forward models. Assume that the actual (nonrealizable) evolution of concentration distribution is where c true 0:τmax contains the concentration distributions at time instants τ = 0, 1, . . . , τ max and ϑ represents actual source parameters. The observations are related to concentrations and the relation can be written as where H • F denotes function composition, H(c true 0:τmax ) is the actual observation model (physics) and W 0:τmax is the observation error that is related directly with the measurement system. However, since accurate models cannot be used or implemented in practice, only approximate models can be used in inversion.
Let us choose a model (H • F )(θ) as an approximate model for the measurement and evolution models. Then, we rewrite the equation (21) in the form where ω(ϑ) = (H • F)(ϑ)) − (H • F • P )(ϑ) is the error due to the approximate models and P : ϑ → θ is a projector from actual model parameters to the parameter space of the approximate model. The term ω(ϑ) is called the approximation error. The calibration data shown in Figure 6 is not appropriate for evaluating the noise correlations between adjacent pixels. Thus, in the computations we replace the actual error W 0:τmax by W 0:τmax that is described in section 3.3.3.
Naturally, since we do not have access to (H • F)(ϑ), we cannot compute ω(ϑ). The approximation error approach can be thought of as the following two steps. The approach adheres to the Bayesian framework and thus an explicit statistical model for all unknowns (here only ϑ, prior model) is to be constructed first. Then, since ω(ϑ) is a function of a random variable, it is a random variable, too. Furthermore, in principle, if all mappings were known, the statistics of ω(ϑ) could be determined and it would be treated as an additional noise term along with the term W 0:τmax .
The second step is the realization (or experience) that we do not necessarily need the physical model, but rather two realizable models H 1 • F 1 and H 2 • F 2 , of which H 1 • F 1 could be described loosely as "an accurate approximation" for the actual physical reality, and H 2 • F 2 is the model intended for the use in inverse computations. Then, we would compute an approximate statistics for ( Typically the actual source model is approximated with a simpler model with parameters ϑ a that are not necessarily compatible with the model H 2 • F 2 . Thus, the model error is approximated as where the model H 1 • F 1 is now written for model parameters ϑ a . In our case both ϑ a and θ = P (ϑ a ) contain the time instant when the leakage is started (τ * ), the release rate (γ) and the index of the source element (k). Due to different discretizations M 1 and M 2 the source element indices are not comparable, and therefore the projection P is defined simply to map the source element from the M 1 to the nearest one in M 2 while for the leakage start time and release rate the projection is identity.
It has turned out that a normal distribution approximation for this is usually sufficient to allow for quite significant model reduction in both stationary [23,22,3] and nonstationary inverse problems [20,19,21]. This approach and its robustness have also been confirmed with real electrical impedance tomography measurements [29]. As noted above, the model reduction is the only considered aspect here, and errors are thus due only to the (limited) accuracy of the forward model that is to be used in the inversion. The mesh information used in the generation of the synthetic data, inversion and in the generation of the model for the approximation error statistics, and the related time steps are given in Table 1.

4.2.
Prior model π(ϑ a ) and the approximate statistics for ω(ϑ). The prior model π(ϑ a ) is chosen so that the elements (of M 1 ) with bottom face crossing the centerline of the pipe are taken to be source elements with equal probabilities. The prior density for the release rate is the uniform density γ ∼ U (0, 4 × 10 −3 ), and the prior probability for the time instants τ * = 0, 1, . . . , τ max are also set as equal. All three variables are set to be mutually independent.
For the computation of the second order statistics (Gaussian approximation) of ω(ϑ a ), we draw a set of samples {ϑ (i) a , i = 1, 2, . . . , 2400} from π(ϑ a ). The forward solution using the models M 1 and M 2 are computed for ϑ (i) a and θ (i) = P (ϑ (i) a ), respectively, and further ω (i) . The sample mean η ω and covariance Γ ω are then computed and we use the Gaussian approximation ω ∼ N (η ω , Γ ω ) in the sequel.
The mean and the element variances are shown in Figure 7. The elements of ω are ordered so that the measurements obtained at the same time form consecutive blocks. Both the mean and variances tend to increase with vector element index (i.e. with time). This is due to the cumulation of discretization errors with time. The significance of these errors depends on their relative strength with respect to the measurement errors.
The likelihood density π(R 0:τmax |θ) is specified by equation (23) and the approximation and observation error models therein. Since ω = ω(ϑ), it is clear that the approximation error ω and the parameters ϑ (and ϑ a and θ) are correlated. The independence of the measurement errors and the approximation errors, however, is usually warranted even though it may have some effects on the results, see e,g, [23,22]. In this paper we neglect the cross-correlations between the parameters and the approximation error. We take the measurement errors as Gaussian where Γ = Γ W0:τ max + Γ ω .
In our case, with the posterior density (26), the maximization problem (27) is identical to the (constrained) minimization problem where L T L = Γ −1 and the constraints are posed by the prior model. The minimization problem (28) can further be written as (29) θ MAP = arg min Inverse Problems and Imaging Volume 3, No. 4 (2009), 711-730 Since there are discrete variables in θ, optimization methods requiring the evaluation of derivatives are not feasible for the computation of MAP estimates. However, they can be computed e.g. with a blind random search algorithm (see [44]), in which the basic principle is to analyze the posterior density by random jumps in the parameter space. The jumps giving increasing values of the posterior density are accepted while the jumps with decreasing values of the posterior density are rejected.
The conditional mean (CM) in our case is defined as The conditional covariance of the random variable θ ∼ π(θ|R 0:τmax ) is In general, the conditional mean and covariance are typically computed via the Markov chain Monte Carlo (MCMC) sampling methods, see [16] for an overview of MCMC methods. In this paper the source parameter space is three-dimensional and τ * and k are discrete variables. Therefore, we are able to evaluate the integrals in (30) and (31) using a numerical quadrature.

5.
Results. Given the observations (see Figure 4) and the mean and covariance of the approximation error (see Figure 7), the objective is to construct the posterior density for the source parameters θ and find estimates and credibility intervals for them. Due to relatively simple parametrization of the source, the posterior density can be easily visualized. The posterior density corresponding to final time instant τ max = 28 is shown in Figure 8a. For comparison, the posterior density was determined also by neglecting the approximation error, i.e. we set Γ ω = 0 and η ω = 0, and the resulting posterior density is illustrated in Figure 8b. Figure 8 shows that the posterior densities differ significantly form each other; the MAP estimates are clearly different and the variances are smaller in the case where the approximation error is neglected.
The actual source parameters and the estimates from the posterior densities which are obtained with and without accounting for the approximation errors are shown in Table 2. The MAP estimates τ * MAP and γ MAP for the leakage start time and release rate with and without approximation error analysis are close to each other. The estimate for the source location, however, has significant difference.
The CM estimate for the leakage start time and the release rate are more accurate than the MAP estimates. Again, for these variables, the adoption of the approximation error model does not have a seemingly major effect. Since the increase of the likelihood variances naturally also increases the posterior variances, the ±1 STD error estimates are smaller in the case in which the approximation error model is not employed. For the estimates that are based on the approximation error analysis, the estimates are slightly better but the significant thing to notice is that the true values lie within ±1 STD limits while the estimates without the approximation errors lie just within ±2 STD limits. This is not, however, a major issue with the τ * and γ variables. The significant problem lies with the estimate for the source location. With the model neglecting the approximation errors, the ±1 STD credibility interval is 11.00 ± 0.01 which essentially claims that the probability of the leakage   11.00 ± 0.01 coming from any other element than the element #11 is vanishing, while the actual source was near the element #16. We refer to such estimates as infeasible. This has been observed in most earlier approximation error related studies. The ability to yield feasible credibility interval estimates is perhaps the most important benefit of employing the approximation error statistics in the overall model.
It is evident that the number of available measurements has an effect on the posterior density and on the estimates. Figure 9 shows the posterior densities corresponding the final time instant τ max = 24, that is, the last 4 observations are not at our disposal. The same data until time τ max = 24 were used as in the previous case. Since the amount of data is decreased, the variances become larger than in Figure 8. The MAP and CM estimates are shown in Table 3. Now the MAP estimates obtained with the approximation error modeling are very close to the true values. When the effect of approximation error is not taken into account, the MAP estimate for γ is a bit poorer and the estimated location is relatively far from the true location. The CM estimates are again rather close to the MAP estimates but, as seen also from Figure 9, the ± 1 STD intervals were wider than  in the previous case. When the approximation error is modeled, the true values of the source parameters are again within the θ CM ± 1 STD intervals. On the other hand, when the approximation error is neglected, the true values of τ * and γ are within the intervals τ * CM ± 1 STD and γ CM ± 1 STD, while the true value of k is well outside the estimated interval k CM ± 1 STD. 6. Conclusions and discussion. In this paper we have considered the identification of a contaminant source in a lake environment by using remote sensing measurements. We proposed a Bayesian approach for the estimation of the properties of contaminant sources. As an extension of previously proposed approaches, we described how the inaccuracies related to mathematical models, in this paper, the model reduction related errors, can be analyzed and taken into account in the construction of the posterior model for the source parameters. Approximation error modeling may play a key role in ensuring the quality and reliability of the estimates especially if heavily reduced models are used.
We implemented a numerical 3D test case in which the objective was to determine the location, release rate and the time instant at which the release was started. We analyzed the discretization error of the numerical models used in inversion, and the resulting error model and an observation error model were used in inversion. The effects of other uncertainties related to the transport and observation models were not considered. The results showed that the differences between the posterior densities obtained with and without approximation error analysis may be significant. Especially, the credibility interval estimates (error estimates) may be infeasibly optimistic when the approximation errors are neglected.
It should be noted that the models used in inversion typically contain also other uncertainties than discretization errors. The effects of these errors can be approximated in a similar way as the discretization errors. For a joint error analysis it is necessary to specify feasible statistical models for all uncertain quantities and generate sets of samples from these models. The samples are then used to solve the forward problem and further the statistics of their effects on the observations. It is obvious that in inverse computations the variances of the resulting posterior densities tend to become larger as the number of uncertain factors increases in the error analysis.