Approximation errors and truncation of computational domains with application to geophysical tomography

Numerical realization of mathematical models always induces errors 
 to the computational models, thus affecting both predictive 
 simulations and related inversion results. Especially, inverse 
 problems are typically sensitive to modeling and measurement errors, 
 and therefore the accuracy of the numerical model is a crucial issue 
 in inverse computations. For instance, in problems related to 
 partial differential equation models, the implementation of a 
 numerical model with high accuracy necessitates the use of fine 
 discretization and realistic boundary conditions. However, in some 
 cases realistic boundary conditions can be posed only for very large 
 or even unbounded computational domains. Fine discretization and 
 large domains lead to very high-dimensional models that may be of 
 prohibitive computational cost. 
Therefore, it 
is often necessary in practice to use coarser discretization and 
smaller computational domains with more or less incorrect 
boundary conditions in order to decrease the dimensionality of the model. 
In this paper we apply the recently proposed approximation error approach 
to the problem of incorrectly posed boundary conditions. 
As a specific computational example we consider the imaging of conductivity 
distribution of soil using electrical resistance tomography. 
We show that the approximation error approach can also be applied to domain truncation 
problems and that it allows one to use significantly smaller scale 
forward models in the inversion.


Introduction
As is common to all inverse problems, the accuracy of the reconstructions is drastically affected by errors associated with the computational forward model that is used to predict the observed variables. Even in cases in which it is possible to construct an accurate mathematical model to describe the behavior of the system of interest, numerical realization of the model introduces modeling errors, which are called the approximation errors [15]. For instance, in problems related to partial differential equations, the approximation errors originate from the discretization and/or an inappropriately selected computational domain for which it is difficult or impossible to impose accurate boundary conditions. The implementation of a numerical model of adequate accuracy may require dense spatial discretization of a very large computational domain. The resulting numerical model is usually very high-dimensional and the computational cost may be excessive for inverse computations. Therefore, it is necessary to implement reduced forward models by using coarser discretization and smaller computational domain. While it is obvious that forward models are approximations of the real system behavior and may thus suffer from various systematic and random modeling errors, the impact of these errors on the estimates is often ignored or accounted for only in an ad hoc manner (e.g., by increasing the variances for the observation errors [28]). Neglecting approximation errors results in a potentially significant estimation bias or other inversion artifacts.
Recently, a novel approach for treating approximation errors has been proposed [15,16]. The basic principle of the approximation error theory is to construct a statistical model describing the approximation errors, and include this information in the inverse calculation. So far, the focus of approximation error modeling has been on the impact of discretization errors. For example, the approach has been successfully applied to inverse problems in optical diffusion tomography [2], demonstrating that modeling approximation errors and taking them into account in data processing allows for the use of a significantly coarser discretization without loss of estimation accuracy.
In this paper the objective is to extend the idea of approximation error theory to modeling the errors due to domain truncation and incorrect boundary conditions. A statistical model for these errors can be constructed in a way that is similar to the approach used for handling discretization errors. As mentioned above, truncating the domain requires specifying appropriate boundary conditions for the subdomain. These boundary conditions typically depend on the solution obtained with the unbounded model, are difficult to determine, and are thus likely to introduce unwanted effects that detrimentally impact the inverse solution. Techniques such as infinite element methods and perfectly matched layer methods have been proposed to overcome the difficulties in handling the domain truncation effects [5,4]. Unfortunately, these methods are not applicable in a general case.
Apart from problems that are induced by partial differential equations, similar types of problems arise in other fields. For example, in image processing data may be missing from either the external boundary or from the inside of the image area. The related deconvolution and in-painting problems have been recently solved using a similar approach as the one proposed in this paper [6,7].
As an example application for modeling the errors due to domain truncation we consider electrical resistance tomography (ERT) employed to geophysical imaging. The objective of geophysical ERT is to obtain information on the subsurface structures and/or hydraulic properties of the soil with electrical measurements [11,20,10]. The problem in constructing the ERT forward model is to determine sufficiently large computational domain so that the standard assumption of negligible currents across the boundary holds true. The effects of domain truncation on estimates from electrical impedance tomography have been examined in [27] in which it was found that the computational domain has to be extended far beyond the region of interest. We construct a statistical model for the errors due to the incorrect boundary condition and discretization in a specific measurement set-up. The benefits of the utilization of the error model are examined in inverse computations.
In geophysical applications, potential sources of approximation errors are numerous. Kowalsky et al. [17] analyzed the systematic errors related to the use of an oversimplified straight-ray model for the simulation of ground-penetrating radar data. They accounted for this error by estimating a constant bias parameter along with hydrogeological and petrophysical parameters during the joint inversion of hydrological and geophysical data. Singha and Gorelick [22] used ERT for the imaging of a saline tracer plume. They noted that the use of standard tomographic inversion techniques yielded a biased estimate of the concentration mass and plume size due to inversion artifacts from regularization and variable measurement sensitivity as well as limitations in Archie's law, which relates electrical conductivities to solute concentrations. Errors in petrophysical models are usually addressed by estimating site-specific relationships [1]. The problem of correlation loss incurred by upscaling core-scale petrophysical relationships to the pixel scale used in ERT and radar-based methods is discussed by Day-Lewis et al. [12]. Linde et al. (2006) [18] analyzed the impact of unaccounted borehole deviations on inverted groundpenetrating radar data. Cooley and Christensen [9] incorporated a second-moment matrix of the model errors into the weight matrix used in weighted nonlinear regression methodologies, thus reducing the estimation bias and improving the assessment of model uncertainties. Finally, Drecourt et al. [14] implemented a bias error covariance matrix into a colored noise Kalman filter and demonstrated its performance in a one-dimensional, saturated groundwater flow problem where a constant and drifting bias was introduced by specifying false boundary conditions. The paper is organized as follows: A brief summary of the ERT observation model is given in Section 2. In Section 3, we review the fundamentals of statistical inversion theory and summarize the approximation error theory. The results from a synthetic ERT inversion are presented in Section 4.

Electrical resistance tomography
Electrical resistance tomography is an imaging method in which the internal structure of the target of interest is imaged based on measurements of electrical current and the related potential. Alternating currents are injected into the target zone through electrodes placed on the boundary or within the object of interest. The currents generate a potential distribution, and potential differences between electrode pairs are measured. The objective of ERT is to reconstruct the electrical conductivity distribution within the target based on the measured potential differences.
The inverse problem needs to be based on a forward model that appropriately relates the conductivity distribution and voltages between the electrode pairs. We use the complete electrode model (CEM) [8,23] as an accurate and realistic representation of the system. The CEM is governed by the following elliptic partial differential equation and boundary conditions: ∂u( r) ∂ n dS = I ℓ , r ∈ e ℓ , ℓ = 1, 2, . . . , N el (3) where r is the position vector, u( r) is the potential distribution in Ω, σ( r) is the conductivity distribution, N el is the number of electrodes, e ℓ is the patch on the (internal or external) boundary ∂Ω corresponding to the location of the ℓ th electrode, z ℓ is the contact impedance between the ℓ th electrode and the object to be imaged, U ℓ is the potential on the ℓ th electrode, and I ℓ the current injected through electrode e ℓ . The currents must satisfy the following charge conservation law: In addition, a reference level for the potential is defined by setting Given the injected currents, the ERT forward model can be formally written using the CEM as where the vector V ∈ R NV contains the measured voltages between N V electrode pairs, and v is the measurement error resulting from non-ideal operation of the data collection system. The inverse problem consists of estimating the conductivity distribution σ( r) given the data V . Since no analytical solutions exist for the geometric conditions encountered in practical applications, the solution of the inverse problem is sought numerically. For this numerical solution, the mapping R(σ( r)) is approximated with a finite-dimensional model Here, σ ∈ R N is a high-resolution discrete approximation for σ( r) using a very large model domain, where N is the number of conductivity values. The approximate model R(σ) : R N → R NV can be constructed using, e.g., the finite element method (FEM), boundary element method (BEM) [26,25], or any other suitable numerical scheme. In cases where the domain Ω is large, accurate numerical implementation may lead to very high-dimensional problems which are difficult or very time-consuming to solve. The problem size can be reduced by either using a coarser discretization and/or by truncating the model domain to a smaller subdomain Ω c ⊂ Ω. In ERT, the parts of Ω where the current density j( r) = −σ( r)∇u( r) is believed to be sufficiently small are typically neglected, and the problem is solved on a truncated subdomain. However, if the selected subdomain is too small, the assumed boundary condition (4) is not accurately fulfilled. Replacing (4) with a realistic boundary condition is difficult or impossible as it depends on the unknown conductivity distribution outside the modeled subdomain. On the other hand, increasing the size of the subdomain may make the problem computationally intractable. In the following section, we propose an approach to account for the errors induced by (a) domain truncation (i.e., imposing incorrect boundary conditions) and (b) by the finite discretization used in the numerical forward model.

Methods
In this section we provide the fundamentals of modeling approximation errors and its implementation. We begin with statistical inversion theory [15,24], which is the basis of the approximation error theory. As mentioned before, we assume that the infinite-dimensional model R(σ( r)) in (5) can be approximated with sufficient accuracy using the finite-dimensional model R(σ) : R N → R NV .
3.1. Statistical interpretation of inverse problems. In statistical inversion theory all variables are considered as random variables. Given the measurements, the primary objective in statistical inversion is to determine the posterior probability density of the quantity of interest. The posterior density can be understood as the solution of the inverse problem, providing the basis for computing various point estimates as well as spread and interval estimates.
Let us assume that the continuous random variables σ ∈ R N and V ∈ R NV have a joint probability density π(σ, V ). According to the rule of conditional densities, the joint density can be expressed as where the marginal densities are From (7) we obtain the conditional probability density which is the well-known Bayes' rule. Once observations are obtained, Bayes' formula allows us to calculate the posterior density. Above we assumed that the joint probability density π(σ, V ) = π(V |σ)π(σ) is known. The density π(V |σ), which is called the likelihood density, is specified by the observation model R(σ) and the statistical properties of the observation noise v. In addition, the marginal density π(σ) must be specified. In statistical inversion, the density π(σ) is called the prior density. The model π(σ) is constructed on the basis of knowledge on the quantity of interest before any measurements are carried out. The construction of the prior density is a critical step; the choice may strongly affect the posterior density, and may introduce some unwanted bias if there are conceptual inconsistencies between the prior information and the corresponding parameters in the observation model.
In the case of additive Gaussian observation noise where N denotes a Gaussian density function, the likelihood density is of the form Further, assuming that π(σ) = N (σ * , Γ σ ) and assuming that v and σ are mutually independent, the posterior density is In high-dimensional problems, posterior densities are difficult to illustrate; the solution is thus reported by computing point estimates. For example, the maximum a posteriori (MAP) estimate is given bŷ Alternatively, the conditional mean (CM) can be calculated aŝ 3.2. Modeling approximation errors. As mentioned in Section 2, minimizing discretization and boundary errors by using an accurate high-resolution model R(σ) with a very large model domain may not be practical because of computational limitations. It is often necessary to use a reduced model, i.e., a model with a smaller computational domain and/or coarser discretization. Let the reduced model be where the subscript c refers to the computational subdomain Ω c ⊂ Ω. Assume that P : R N → R Nc is a linear mapping from R N to the reduced space, i.e., σ c = P σ. Then, the exact model can be written as The term ǫ(σ) describing the discrepancy between the accurate model and the reduced model is usually neglected when constructing the posterior density. Depending on the discretization and computational domain associated with the reduced model, the model error ǫ(σ) may have a significant effect on the solution of the inverse problem [15]. In a statistical framework, the term ǫ(σ) is considered a random variable because it depends on σ ∼ π(σ), which is a random variable. It is evident that by analyzing the error ǫ(σ) and taking it into account in the reconstruction, the accuracy of the estimates can be improved. Given the probability densities π(σ) and π(v), the probability density of the random variable η(σ, v) = ǫ(σ) + v is Since the forward models R : R N → R NV and R c : R Nc → R NV are nonlinear, the general expression (8) is difficult to implement, and is thus approximated by a Gaussian distribution, π(η) ∼ N (η * , Γ η ). Assume that a prior probability density π(σ) and draw a set of sample conductivities S = σ (k) ∈ R N , k = 1, . . . , M . The set S is used to approximate the mean and covariance of η as Given the Gaussian prior density π(σ c ) ∼ N (σ c * , Γ c ) and the approximation π(η) ∼ N (η * , Γ η ), which is called the enhanced error model [15], and assuming that the error η is statistically independent of the conductivities σ c , the posterior density is

Model development
The approach is applied to geophysical ERT, where we construct a statistical model for the errors due to domain truncation and discretization. In addition, we demonstrate the advantages of the enhanced error model for numerical simulations and inversions. For simplicity, we consider two-dimensional model domains; the implementation in three dimensions is straightforward.

4.1.
Simulation study set-up and computational domains. The synthetic ERT measurement system consists of two vertical boreholes 8 m apart. An array of electrodes is placed into each borehole, and the distance between the electrodes in each borehole is 1.2 m. The electrodes are represented as internal boundaries applied to small "holes" in the computational domain; the boreholes themselves are not explicitly discretized, as they would block electric current in this two-dimensional model (see Figure 3   domain Ω 0 is considered to be large enough for simulating the ERT measurements using the boundary condition (4). Finally, a finely discretized domain Ω 4 of size 150 m × 70 m is created and used for the generation of synthetic observation data for the inversion. The computational domains are shown in Figure 1. The current density on the boundary ∂Ω 1 is visualized in Figure 2.
The domains are discretized using triangular elements. Two different meshes are created for the domains Ω 0 , Ω 1 , Ω 2 and Ω 3 : a fine, unstructured mesh M F i , (i = 0, 1, 2, 3) for the simulation of the potential distribution, and a coarser, more uniform mesh M C i , in which electric conductivities at each nodal point are to be estimated. A measure of relative discretization density compared to the "accurate" mesh Ω 0 is given as where N F i is the number of nodal points of mesh M F i , and N F 0 Ωi is the number of nodal points of mesh M F 0 in the region covered by subdomain Ω i . The details on the domains and meshes are given in Table 1. The mesh M F 3 is shown in Figure 3 and mesh M C 3 is shown in Figure 4.

4.2.
The Gaussian model of the approximation error. We followed the approach described in Section 3.2 to assess the effect of domain truncation and coarse discretization on the ERT observations and the resulting reconstruction of the electric conductivity field. First, a set of samples from the prior density needs to be generated. A commonly used prior model in geophysics is the smoothness  prior model [13,19]. The smoothness prior density can be written as    where D is the first-order difference matrix, α > 0 is the regularization parameter that is used to control the level of smoothness, and αD T D = Γ −1 σ is the inverse prior covariance. Unfortunately, the smoothness prior density (12) cannot be used for sample generation since Γ −1 σ is rank-deficient and thus improper. However, by conditioning the density (12) on the conductivity at certain nodal points it is possible to obtain a proper density for generating samples for the approximation error analysis [15].
By reordering the indices such that the conditioning points are assembled at the end, the discrete conductivity distribution can be written as σ = (σ 1 ; σ 2 ), where σ 1 contains the sampled conductivities, and σ 2 represents the conductivity values at the conditioning nodes. Accordingly, matrix D T D can be partitioned as and the conditional smoothness prior can be written as By setting the proper joint density can be written as To evaluate approximation errors, multiple samples of the conductivity field are drawn from the density (14). A first-order difference operator on the uniform mesh M C 0 is created, as described in [21]. We operate in a uniform mesh in order to keep the smoothness properties of the samples uniform throughout the domain Ω 0 . We then randomly choose 20 nodes (shown in Figure 5) for which we specified a Gaussian model (13) with σ 2 * = 0.03 and Γ σ2 = (0.15 × 10 −5 )I. The smoothness parameter α is set to 10 4 . A total of 500 samples are drawn from the proper prior density; four samples are shown in Figure 5. Each draw σ (k) is mapped onto the accurate mesh M F 0 and to the truncated, coarser meshes M F i , i = 1, 2, 3. For each of the domains and each of the 500 realizations of the conductivity field, current injection and voltage observations are simulated using a FEM approximation of the complete electrode model, with the homogeneous Neumann condition (4) imposed at the boundaries of the truncated models. Injection of current is simulated for 16 pairs of electrodes, and the corresponding voltages are recorded. As shown in Figure 6, the first eight injections and voltage measurements are made in a crosshole pattern between the electrodes in the two boreholes, followed by measurements among four electrode pairs within each of the two boreholes.
The approximate, truncated observation models were evaluated and compared to the accurate model M F 0 at the sample points, thus evaluating the approximation error ǫ(σ) for each realization of the underlying conductivity field: (15) ǫ(σ (k) ) = R(P 0 σ (k) ) − R c (P c σ (k) ), k = 1, 2, . . . , 500, where P 0 is an interpolation matrix from the mesh M C 0 to M F 0 , P c is an interpolation matrix from M C 0 to the associated fine mesh in the truncated domain Ω c = Ω i , i = 1, 2, 3, and R c is the approximate forward model. The mean and covariance of the approximation errors related to each truncated mesh were obtained by evaluating (9) and (10).  The approximation errors are illustrated in Figure 7, which shows the mean and covariance of the error related to current injection No. 8 (see Figure 6) for the successively larger subdomains Ω 1 , Ω 2 , and Ω 3 . Truncating the infinite-acting system (represented by domain Ω 0 ) by imposing no-flow conditions on boundary ∂Ω 1 leads to a significantly larger voltage difference between the injection electrode pair No. 8. The approximation errors as calculated by (15) and the corresponding mean η * (9) show a relatively strong negative bias of −20 V for the eighth electrode pair used for current injection. The bias is somewhat diminished for the other electrode pairs, which are used for voltage observations only, and becomes slightly positive for the four electrode pairs within the left borehole (pairs No. 9-12). As expected, the mean approximation error decreases as the size of the truncated subdomain increases from Ω 1 to Ω 3 . Moreover, the covariance matrix, which reflects the uncertainty induced by the variability of the underlying conductivity field and its impact on the current density near the boundaries, also decreases with increasing size of the subdomain.

Results and discussion
In this section we demonstrate the benefits of using the approximation error model for subsurface imaging. In the first example, the conductivity distribution to be reconstructed consists of a circular region with higher conductivity in an otherwise homogeneous background. The circular region is located in the center between the two boreholes, which is the region most difficult to image using ERT. To avoid committing an "inverse crime" (see [15]), two different meshes are used for data generation and the modeling of the approximation error. The rectangular domain Ω 4 of size 150 m × 70 m is discretized using triangular elements, and the mesh is denoted as M F 4 . The relative number of elements and nodes in this mesh is larger than in the mesh M F 0 , which represents the "accurate" model in the construction of the statistical models for the approximation errors. Using a different mesh M F 4 allows us to account for the minor differences between the "actual" behavior of the system and the computational model used in the evaluation of approximation errors.
The measurements are simulated in the mesh M F 4 using the same current injection and measurement patterns as described in Section 4.2. The computed voltages are corrupted with additive noise, which consists of two components, both being Gaussian with zero mean. The standard deviation of the first component is 1 % of the absolute noiseless voltage, and the standard deviation of the second component is 0.1 % of the difference between the maximum and minimum voltages.
Using the simulated observation data, the estimates for the actual conductivity distribution are computed with approximate ERT forward models. In the following notations, subscript c refers to the models and quantities associated to the truncated domains Ω c = Ω i , i = 1, 2, 3. In the case of a smoothness prior density, the posterior density is of the form where L T η L η = Γ −1 η and D c is a first-order difference matrix. In (16), all models and quantities are related to the fine mesh M F i . The MAP estimate can be obtained using the Gauss-Newton algorithm.
For the inversion, conductivities are estimated for the uniform meshes M C i , i = 1, 2, 3 (see Table 1) to reduce the dimensionality of the problem. However, in the Gauss-Newton recursions, the forward models were evaluated using the fine mesh. Thus, by denoting the estimate in the coarse mesh asσ c , the recursion is of the form  Figure 7. The sample mean and the associated covariance of the approximation error corresponding to current injection between the eighth electrode pair (see Figure 6).
The mapping R c : is the Jacobian matrix of the forward model, P is an interpolation matrix from the coarse to the fine mesh in domain Ω c , the matrix L η contains the Cholesky factors of the inverse covariance matrix Γ −1 η , andD c is the first-order difference matrix of mesh M C i . The step length parameter κ (j) is computed using the method of Barzilai and Borwein [3]. The covariance Γ v is estimated on the basis of the observations assuming that the noise level are approximately known. The standard deviation of the first component is assumed to be 1 % of the absolute noisy voltage, and the standard deviation of the second component was assumed to be 0.2 % of the difference between the maximum and minimum noisy voltages.
For comparison, we also compute conductivity estimates by neglecting the effect of the approximation errors, i.e., we set η * = v * and Γ η = Γ v in the posterior density (16). The actual conductivity distribution and the estimates are shown in Figure  8. The same voltage data are used in all inversions. In each case, the smoothness parameter α is selected by visual inspection so that the smoothness properties of the estimates are acceptable and, consequently, the conductivity on the circular region is reconstructed appropriately. Figure 8 shows the true conductivity field as well as the estimated distributions with and without modeling of the approximation errors for the three truncated subdomains Ω i , i = 1, 2, 3. Consistent with the previous results shown in Figure 7, the accuracy of the estimates improves as the size of the subdomain increases. This was expected since the homogeneous Neumann boundary condition (4) is more justified for a larger domain size.
Accounting for approximation errors due to domain truncation yields significantly improved estimates of the actual distribution. Neglecting the approximation errors leads to a substantial over-estimation of conductivity near the truncated boundaries, and a poor reconstruction of the circular conductivity anomaly in the center of the domain. The estimates obtained with the enhanced noise model are closer to the actual distribution near the domain boundaries. The circular region with larger conductivity can be relatively well located, but the conductivities in this region are somewhat underestimated as a result of the smoothness prior model used in these inversions.
In a second case, we study potential effects of a conductivity anomaly that is located immediately outside the truncated domain. The actual conductivity distribution constructed for data generation is similar to that used in the previous case, but an additional circular anomaly of higher conductivity is placed outside the right vertical boundary of the domain Ω 1 . The observations were generated as before, and the reconstruction are computed using the convectional and enhanced noise models. The actual conductivity distribution and estimates are shown in Figure 9.
While the figure shows the large domain with both anomalies, the inversion is restricted to the subdomain Ω 1 , which is indicated by the rectangle. The estimates are very similar to those shown in Figure 8, i.e., the more conductive region outside the computational domain does not seem to have a significant effect on the estimates.
In a third case, we study effects of number of draws from the proper smoothness prior model that is used to generate Gaussian distribution of approximation error. The actual conductivity distribution constructed for data generation is similar to that used in the first case. The observations are also generated as in the previous cases, and the reconstructions are computed using the enhanced noise models. The actual conductivity distribution and estimates are shown in Figure 10.

Conclusions
In solving inverse problems the accuracy of the computational forward model may have a significant effect on the estimates. The errors in the forward model may be due to the discretization required for numerical implementation, an inappropriately chosen computational domain, incorrect boundary conditions, or an inaccurate mathematical model that is used to describe the phenomena of interest. Implementation of accurate, high-resolution forward models usually leads to very  high-dimensional inverse problems and thus the computational burden may become prohibitive.
In this paper we focused on the errors due to the truncation of the computational domain and discretization with an application in geophysical ERT. The principles of the approximation error theory presented in this paper are not restricted to ERT, but are applicable to a wide variety of inverse problems. The fundamental principle is to construct a statistical model for the difference between the observations computed with the accurate and approximate forward models on the basis of prior knowledge of the target to be imaged.
We carried out an ERT simulation study which showed that by employing an appropriate Gaussian statistical model for the approximation errors the accuracy of the estimates can be significantly improved by reducing estimation bias and uncertainty. Modeling the approximation errors makes possible the use of less accurate and computationally less intensive forward models in an inversion without significant loss of estimation accuracy.
It should be pointed out that in this study the enhanced error models were constructed using a smoothness prior model. Thus, the samples drawn from the prior density were quite different from the conductivity distributions to be imaged, which consisted of a discrete anomaly with sharp boundaries. Nevertheless, the quality of the estimates was significantly improved, demonstrating the robustness of the approach even if prior information is limited.