Bayesian inversion in resin transfer molding

We study the Bayesian inverse problem of inferring the permeability of a porous medium within the context of a moving boundary framework motivated by Resin Transfer Molding (RTM), one of the most commonly used processes for manufacturing fiber-reinforced composite materials. During the injection of resin in RTM, our aim is to update our probabilistic knowledge of the per- meability of the material by inverting pressure measurements as well as observations of the resin moving domain. We consider both one-dimensional and two-dimensional forward models for RTM. Based on the analytical solution for the one-dimensional case, we prove existence of the sequence of posteriors that arise from a sequential Bayesian formulation within the infinite-dimensional framework. For the numerical characterisation of the Bayesian posteriors in the one-dimensional case, we investigate the application of a fully-Bayesian Sequential Monte Carlo method (SMC) for high-dimensional inverse problems. By means of SMC we construct a benchmark against which we compare performance of a novel regularizing ensemble Kalman algorithm (REnKA) that we propose to approximate the posteriors in a computationally efficient manner under practical scenarios. We investigate the robustness of the proposed REnKA with respect to tuneable param- eters and computational cost, and display advantages of REnKA compared with SMC with a small number of particles. We further apply REnKA to investigate, in both the one-dimensional and two-dimensional settings, practical aspects relevant to RTM which include the effect of pressure sensors configuration and the observational noise level in the uncertainty in the log-permeability quantified via the sequence of Bayesian posteriors.

1. Introduction. In this paper we study the Bayesian inverse problem within the moving boundary setting motivated by applications in manufacturing of fiber-reinforced composite materials. Due to their light weight, high strength, as well as their flexibility to fit mechanical requirements and complex designs, such materials are playing a major role in automotive, marine and aerospace industries [4,3,26]. The moving boundary problem under consideration arises from Resin Transfer Molding (RTM) process, one of the most commonly used processes for manufacturing composite materials. RTM consists of the injection of resin into a cavity mold with the shape of the intended composite part according to design and enclosing a reinforced-fiber preform previously fabricated. The next stage of RTM is curing of the resin-impregnated preform, which may start during or after the resin injection. Once curing has taken place, the solidified part is demolded from the cavity mold. In the present work we are concerned with the resin injection stage of RTM under the reasonable assumption that curing starts after resin has filled the preform. Though the current study is motivated by RTM, the results can be also used for other applications where a moving boundary problem is a suitable model.
We now describe the (from the inverse problem prospective, forward) model (see further details in [3,46,37]). Let D * ⊂ R d , d ∈ {1, 2}, be an open domain representing a physical domain of a porous medium with the permeability κ(x) and porosity ϕ. The boundary of the domain D * is ∂D * = ∂D I ∪ ∂D N ∪ ∂D O , where ∂D I is the inlet, ∂D N is the perfectly sealed boundary, and ∂D O is the outlet. The domain D * is initially filled with air at a pressure p 0 . This medium is infused with a fluid (resin) with viscosity µ through an inlet boundary ∂D I at a pressure p I and moves through D * occupying a time-dependent domain D(t) ⊂ D * , which is bounded by the moving boundary Υ(t) and the appropriate The forward problem for the pressure of resin p(t, x) consists of the conservation of mass x ∈ D(t), t > 0, (1.1) where the flux v(x, t) is given by Darcy with the following initial and boundary conditions p(x, t) = p I , x ∈ ∂D I , t ≥ 0, (1.3) ∇p(x, t) · n(x) = 0, x ∈ ∂D N , t ≥ 0, (1.4) V (x, t) = − κ(x) µϕ ∇p(x, t) · n(x, t), x ∈ Υ(t), t ≥ 0, (1.5) p(x, t) = p 0 , x ∈ Υ(t), t > 0, (1.6) p(x, t) = p 0 , x ∈ ∂D O , t > 0, (1.7) p(x, 0) = p 0 , x ∈ D * , (1.8) Here V (t) is the velocity of the moving boundary Υ(t) in the normal direction, n(x) and n(x, t) are the unit outer normals to the corresponding boundaries. We note that in the considered one and two dimensional cases of this problem, we can view the velocity of the moving boundary as the following derivative: V (x, t) = dΥ(t) dt · n(x, t). (1. 10) We remark that for definiteness we have assumed that at the initial time the moving boundary Υ(0) coincides with the inlet boundary ∂D I and that the constant pressure condition is imposed at the inlet. It is not difficult to carry over the inverse problem methodology considered in this paper to other geometries and other conditions on the inlet (e.g. constant rate). Further, in two (three) dimensional RTM settings one usually models permeability via a second (third)-order permeability tensor to take into account anisotropic structure of the media [3,37] but here for simplicity of the exposition the permeability κ(x) is a scalar function. Again, the developed methodology is easy to generalize to the tensor case.
Let us note that in the one-dimensional case the nonlinear problem (1.1)-(1.9) is analytically simple and admits a closed form solution (see Section 2 and [37]) but the two and three dimensional cases are much more complicated and analytical solution is in general not available. We remark that in two and three dimensional cases the resin can race around low permeability regions and the front Υ can become discontinuous creating macroscopic voids behind the main front (see further details in [3,37]) but in this paper we ignore such effects which deserve further study.
It has been extensively recognized [13,14,30,29,36,41] that imperfections in a preform that arise during its fabrication and packing in the molding cavity can lead to variability in fiber placement which results in a heterogenous highly-uncertain preform permeability. In turn, these unknown heterogeneities in permeability of the preform give rise to inhomogeneous resin flow patterns which can have profound detrimental effect on the quality of the produced part, reducing its mechanical properties and ultimately leading to scrap. To limit these undesirable effects arising due to uncertainties, conservative designs are used which lead to heavier, thicker and, consequently, more expensive materials aimed at avoiding performance being compromised. Clearly, the uncertainty quantification of material properties is essential for making RTM more cost-effective. One of the key elements in tackling this problem is to be able to quantify the uncertain permeability.
In this work we assume that D * , ∂D O , ∂D I , ∂D N , p I , p 0 , µ and ϕ are known deterministic parameters while the permeability κ(x) is unknown. Our objective is within the Bayesian framework to infer κ(x) or, more precisely, its natural logarithm u(x) = log κ(x) from measurements of pressure p(x, t) at some sensor locations as well as measurements of the front Υ(t), or alternatively, of the time-dependent domain D(t) at a given time t > 0. We put special emphasis on computational efficiency of the inference, which is crucial from the applicable point of view.
1.1. Practical approaches for permeability estimation in fiber-reinforced composites. While the estimation of preform permeability during resin injection in RTM is clearly an inverse problem constrained by a moving boundary PDE such as (1.1)-(1.9), most existing practical approaches pose the estimation of permeability in neither a deterministic nor stochastic inverse problems framework. For example, the very extensive review published in 2010 [40] reveals that most conventional methods for measuring permeability assume that (i) the material permeability tensor is homogenous and (ii) the flow is one-dimensional (including 2D radial flow configurations). Under these assumptions the resin injection in RTM can be described analytically, via expressions derived from Darcy's law, which enable a direct computation of the permeability in terms of quantities that can be measured before or during resin injection. These conventional methods suffer from two substantial practical limitations. First, they do not account for the heterogenous structure of the preform permeability, and although they provide an estimate of an effective permeability, this does not enable the prediction of the potential formation of voids and dry spots. Second, those conventional methods compute the permeability in an off-line fashion (i.e before RTM) with specific mold designs that satisfy the aforementioned assumptions intrinsic to those methods (e.g. rectangular flat molds). This second limitation is not only detrimental to the operational efficiency of RTM but also neglects the potential changes in permeability that can results from encapsulating the preform in cavities with complex designs. Some practical methodologies for online (i.e. during resin injection) estimation of heterogenous permeability have been proposed in [35,48]. While these approaches seem to address the aforementioned limitations of conventional methods, they also use a direct approach for the estimation of permeability which faces unresolved challenges. As an example, let us consider the recent work of [48] which uses an experimental configuration similar to the one described in Figure 1.1 and which, by using pressure mea-surements from sensors located within the domain occupied by the preform, computes a finite-difference approximation of the normal flux to the front ∇p(Υ(t), t) · n. In addition, by means of images from CCT cameras, seepage velocity of the resin front is computed in [48]; this velocity is nothing but V (x, t) defined by (1.5) in the context of the moving boundary problem (1.1)- (1.9). Under the assumption that µ and ϕ are known, the approach proposed in [48] consists of finding κ(Υ(t)) = arg min θ V (x, t) + θ ϕµ ∇p(Υ(t), t) · n (1. 11) with V (x, t) and ∇p(Υ(t), t) · n computed from measurements as described above. This approach offers a practical technique to estimating κ on the moving front and can then potentially infer the whole permeability field during the resin injection in RTM. However, from the mathematical inverse problems perspective, this ad-hoc approach is not recommended as it involves differentiating observations of pressure data for the computation of ∇p(Υ(t), t) · n. Indeed, it is well-known [22] that differentiation of data is an ill-posed problem that requires regularization. In addition, rather than an inverse problem, the least-squares formulation in (1.11) is a data fitting exercise that excludes the underlying constraint given by the moving boundary problem and which entails a global effect induced by κ. As a result, the estimate of permeability obtained via (1.11) has no spatial correlation and thus fails to provide an accurate global estimate of the permeability field.
The recent work of [28] demonstrates considerable advantages of using systematic data assimilation approaches to infer permeability during the resin injection of RTM. By means of a standard ensemble Kalman methodology for data assimilation, the approach of [28] uses measurements from visual observations of the front location to produce updates of the preform permeability within the context of a discrete approximation of the moving boundary problem (1.1)-(1.9). While the methodology used in [28] is focused in producing deterministic estimates, the standard Kalman methodology can be potentially used to quantify uncertainty in preform permeability. However, it has been shown that standard Kalman methodologies, such as the one used in [28], could result in unstable estimates unless further regularisation to the algorithm is applied [18].
In addition to the lack of an inverse problem framework that can lead to unstable and ultimately inaccurate estimates of the permeability in resin injection of RTM, most existing approaches (i) do not incorporate the uncertainty in the observed variables and (ii) do not quantify uncertainty in the estimates of the permeability of preform. It is indeed clear from our literature review that the estimation of permeability of preform during resin injection deserves substantial attention from an inverse problems perspective capable of quantifying uncertainty inherent to the fabrication and packing of the preform.
1.2. The Bayesian approach to inverse problems. In this paper we propose the application of the Bayesian approach to inverse problems [45] in order to infer the logarithm of the permeability u(x) = log κ(x), from observations {y n } N n=1 collected at some prescribed measurement/observation times {t n } N n=1 during the resin injection in RTM. At each time t n we observe a vector, y n , that contains noisy measurements of resin pressure from sensors as well as some information of the moving domain (or alternatively front location) observed, for example, via CCT cameras or dielectric sensors [31]. In the Bayesian approach, the unknown u(x) is a random function that belongs to a space of inputs X. A prior probability measure µ 0 (u) = P(u) on u must be specified before the data are collected; this enables us to incorporate prior knowledge which may include design parameters as well as the uncertainty that arises from preform fabrication (i.e. prior to resin injection). In our work we consider Gaussian priors which have been identified as adequate for characterizing the aforementioned uncertainty in log-permeability from the preform fabrication [49,30,29] (see also references therein).
At each observation time t n during the infusion of resin in RTM, we then pose the inverse problem in terms of computing, µ n (u) = P(u|y 1 , . . . , y n ), the (posterior) probability measure of the log-permeability conditioned on measurements y 1 . . . , y n . Each posterior µ n then provides a rigorous quantification of the uncertainty in the log-permeability field given all available measurements up to the time t n . Knowledge of each of these posteriors during RTM can then be used to compute statistical moments of the logpermeability under µ n (e.g. mean, variance) as well as expectations of quantities of interest that may be needed for the optimization of controls (e.g. pressure injection) in RTM.
Although the proposed application of the Bayesian formulation assumes Gaussian priors, the nonlinear structure of the PDE problem, that describes resin injection in RTM, gives rise to a sequence of non-Gaussian Bayesian posteriors {µ n } N n=1 which cannot be characterized in a closed form. A sampling approach is then required to compute approximations of these posteriors. Among existing sampling methodologies, Sequential Monte Carlo (SMC) samplers [33,8,1,21] are particularly relevant for the formulation of the above described inverse problem as they provide a recursive mechanism to approximate the sequence of Bayesian posteriors {µ n } N n=1 . Starting with J samples from the prior u that approximates µ n as the new data y n collected at time t n become available. The weights {W n > 0) and the empirical measure converges to µ n as J → ∞ (δ w denotes the Dirac measure concentrated at w). Moreover, if f (u) denotes a quantity of interest of the unknown log-permeability u(x), the weighted particles {W can be easily used to compute the sample mean (1.13) which converges (see for example [33]) to the expectation (under µ n ) of the quantity of interest E µn (f (u)).
The recursive computation of the weighted particles in SMC is suitable for the proposed application in RTM as it allows us to update, potentially in real time, our knowledge of the uncertainty in the log-permeability. However, producing accurate approximations of the Bayesian posteriors {µ n } J n=1 in the context of the inference of preform log-permeability in RTM represents a substantial computational challenge that arises from the fact that these posterior measures are defined on a (infinite-dimensional) functional space. Upon discretization, these posteriors could be potentially defined on a very highdimensional space. Unfortunately, it has been shown [9,5] that standard Bayesian sampling methodologies such as standard SMC do not scale well with the dimension of the (discretized) unknown; this leads to unstable and ultimately inaccurate algorithms.
The recent works [9,21] developed scalable (dimension independent) sampling algorithms for the approximation of the Bayesian posterior that arises from high-dimensional inverse problems. While these algorithms have a solid theoretical background that ensures their stability and convergence properties, achieving a desirable level of accuracy often comes at extremely high computational cost. More specifically, Bayesian methodologies, that provide approximation of the form (1.12) and that converge asymptotically to the underlying posterior measure µ n , often involve solving the forward model thousands or even millions of times. In the context of the inverse problem for RTM, the numerical solution of the moving boundary (forward) problem in 2D or 3D settings is computationally very intensive. Therefore, the sequential approximation of the Bayesian posteriors of preform's log-permeability must be conducted with scalable computational efficiency so that it can be realistically used within a near real-time optimization loop for RTM. In the proposed work we develop a computational inverse framework that possess such computational efficiency with the ultimate aim of the real-time uncertainty quantification of the reinforced preform's log-permeability.
1.3. Contributions of this work. The contributions of this article are the following: (A) A Bayesian formulation of the inverse problem to infer log-permeability from sequential data collected during resin injection in RTM. Both the 2D forward model described by (1.1)-(1.9) as well as the corresponding 1D version are considered. For the 1D case, we show that application of the infinite-dimensional Bayesian framework of [45] leads to well-posedness of the sequence of Bayesian posteriors. (B) Application of a state-of-the-art SMC framework [21] for the approximation of the sequence of Bayesian posteriors {µ n } J n=1 that arises from the Bayesian formulation. From this SMC framework, we motivate a novel regularizing ensemble Kalman algorithm (REnKA) that aims at approximating this sequence of posteriors in a computationally efficient manner, thus suitable for its implementation in a practical setting of RTM. (D) Numerical investigation of the accuracy and robustness of the proposed REnKA scheme in the 1D case; this involves constructing, via the SMC sampler of [21], accurate approximations of the posteriors that we use as Benchmark against which we compare the proposed REnKA. The advantages of REnKA in terms of accuracy vs computational cost are showcased by comparing it with the implementation of a low-resolution SMC whose computational cost is comparable to REnKA's. (E) Application of REnKA for further investigations of the Bayesian inverse problem in both 1D and 2D. In particular for the 1D case we conduct a numerical investigation of the added value of assimilating the front location relative to the number of pressure sensors. Since the number of pressure sensors that can be physically deployed for preform permeability monitoring in RTM is usually limited, this investigation aims at providing practitioners with guidelines for the number of sensors that can accurately infer preform permeability alongside with its uncertainty. In addition, for the 1D case we study the effect of the frequency of the observations, as well as the observational noise level on the inference of the log-permeability. We further apply REnKA to the 2D forward model and, analogous to the 1D case, we study the effect that the number of pressure sensors have on the inferred log-permeability. The rest of the paper is organized as follows. In Section 2 we introduce the Bayesian inverse problem of inferring the permeability of a porous media in a 1D moving boundary problem for resin injection in RTM. In Section 3 we discuss and apply SMC to approximate the Bayesian posteriors that arise from the Bayesian approach. In Section 4 we introduce REnKA and conduct a numerical investigation of its approximation properties relative to its computational cost. In Section 5 we apply REnKA to further investigate relevant practical aspects of the inverse problem in both 1D and 2D; this includes the study of the effect of the number of pressure sensors as well as the noise level on accuracy of the inferred log-permeability and its uncertainty. Some conclusions are presented in Section 6.
2. Bayesian inversion of a one-dimensional RTM model. In this section we apply the Bayesian approach to infer log-permeability in the context of the one-dimensional version of the forward problem defined in (1.1)-(1.9). The corresponding 1D moving boundary problem induces a sequence of forward maps that we define in Section 2.1 and that we aim at inverting with the Bayesian formalism that we introduce in Section 2.2. This sequence of 1D forward maps admits a closed-form solution that can be numerically approximated at a very low computational cost. This will enable us in Section 3 to obtained accurate numerical approximations of the solution to the Bayesian inverse problem; we use this accurate approximations as a benchmark for assessing the approximation properties of the ensemble Kalman algorithm that we introduce in Section 4 .
2.1. The Forward 1D RTM model. Let us consider a one-dimensional porous media with physical domain D * ≡ [0, x * ] ⊂ R. As before, we denote by κ(x) (x ∈ D * ) and φ > 0 the permeability and porosity of the porous medium, respectively. Resin with viscosity µ is injected at x = 0 at a pressure p I . The pressure at the moving front (outlet) Υ(t) is prescribed and equal to p 0 . The initial pressure distribution before injection is also set to p 0 . For convenience of the subsequent analysis, we parameterize the permeability in terms of its natural logarithm u(x) ≡ log κ(x). The pressure p(x, t) and the moving front Υ(t) are given by the solution to the following model d dx The solution to (2.1)-(2.5) can be obtained analytically by the following proposition (see [3,46,37]).
The unique solution Υ(t), p(x, t) of (2.1)-(2.5) is given for t ≥ 0 by The quantity of interest arising from the RTM injection model is the so-called filling time: the time it takes the front Υ(t) to reach the right boundary of the domain of interest [0, x * ]. Filling time, denoted by τ * , is defined by Υ(τ * ) = x * . From (2.7) and the definition in (2.6) it follows [37] that τ * is given by Note that the parameters p 0 and p I are prescribed control variables and thus known. In addition, we assume that µ and φ are known constants. As stated earlier, we are interested in the inverse problem of estimating the permeability, or more precisely its natural logarithm u(x) = log κ(x) given time-discrete measurements of the front location as well as the pressure from M sensors located at {x m } M m=1 ⊂ [0, x * ]. We denote by {t n } N n=1 the set of N observation times. For fixed (assumed known) parameters p I , p 0 , φ and µ, the solution to the PDE model (2.1)-(2.5) induces the nth forward map G n : Given u(x) = log κ(x) ∈ X, the evaluation of the forward map G n (u) predicts the location of the front and the pressure at the sensor locations at the time t = t n . Since observation times are prescribed before the experiment, there is no assurance that for a given u, the corresponding filling time satisfies t n ≤ τ * for all n = 1, . . . , N . In other words, the front could reach the right end of the domain before we observe it at time t n . In a real experimental setting, the process stops at time τ * . However, in the inverse problem of interest here, observation times are selected beforehand, and the search of optimal u's within the Bayesian calibration of the nth forward map can lead to filling times greater than some observation times. In this case (t n > τ * ), the definition (2.10) yields G n (u) = [Υ(τ * ), {p(x m , τ * )} M m=1 ] T . The following theorem ensures the continuity of the forward map, which is necessary for justifying the application of the Bayesian framework in Section 2.2.2.
For the proof of this theorem, see Appendix A.
In the following subsection we apply the Bayesian framework for inverse problems in order to invert observations of G n (u).
Remark 2.1. We note that for the present work the porosity ϕ is an assumed known constant; our objective is to infer the log-permeability u(x) = log κ(x). However, the Bayesian methodology that we apply can be extended to the case where the unknown is not only log κ(x) but also ϕ, and can include the case where ϕ = ϕ(x) is a spatial function defined on the physical domain D * .
2.2. The Bayesian Inverse Problem. Suppose that, at each observation time t = t n , we collect noisy measurements of the front location as well as pressure measurements from sensors. We denoted these measurement by y Υ n ∈ R + and y p n ∈ R M , respectively. Our aim is to solve the inverse problem of estimating the log permeability u(x) = log κ(x) given all the data y p 1 , y Υ 1 , . . . , y p n , y Υ n up to time t = t n . We assume that the aforementioned observations are related to the unknown u(x), via the forward map (2.10), in terms of y p n = G p n (u) + η p n , (2.11) where η Υ n and η p n are realizations of Gaussian noise with zero mean and covariance Γ Υ n and Γ p n , respectively, i.e. η Υ n ∼ N (0, Γ Υ n ) and η p n ∼ N (0, Γ p n ) (i.i.d.). For simplicity we assume that both measurements of the front location and pressures are uncorrelated in time. We additionally assume that η Υ n and η p n are uncorrelated for all n = 1, . . . , N .
Note that (2.11)-(2.12) can be written as (2.14) Remark 2.2. Due to the nature of the RTM problem, we have that the pressure p(x m , t) at each sensor x m should increase with time as well as the fact that G Υ n+1 (u) ≥ G Υ n (u). However, the Gaussian noise in (2.11)-(2.12) can make the observations y p n and y Υ n "unphysical". In practice, observations need to be post-processed before using them for the Bayesian inverse problem and unphysical y p n , y Υ n should be excluded. We leave the question of how to incorporate such a post-processing framework for future study.
Here we follow the traditional point of view on data modeled via (2.11)-(2.12) and choose sufficiently small Γ Υ n and Γ p n so that the probability of y p n , y Υ n being unphysical is very low.
We adopt the Bayesian framework for inverse problems where the unknown u(x) = log κ(x) is a random field and our objective is to characterize the sequence of distributions of u conditioned on the observations which we express as u|y 1 , . . . , y n . In other words, at each observation time t = t n we aim at computing the Bayesian posterior µ n (u) = P(u|y 1 , . . . , y n ). From this distribution we can obtain point estimates of the unknown that can be used in real time to, for example, modify controls (e.g.. p I ). More importantly, as we stated in the Introduction, the aforementioned distribution enables us to quantify uncertainty not only of the unknown but also of quantities of interest that may be relevant to an optimization of resin injection in RTM.
Even though for the illustrative purposes the model presented in this section is discretized on a relatively low dimensional space (e.g. 60 cells), our aim is to introduce a general computational framework independent of the size of the discretized domain. We therefore consider an infinite-dimensional formulation of the Bayesian inverse problem for which the unknown u belongs to a functional space X. The discretization of the Bayesian inverse problem will be conducted at the last stage of the computational algorithm, when the posteriors are sampled/approximated. Thus, we are aiming at robust mesh-invariant computational algorithms.

The Prior.
For the Bayesian approach that we adopt in this work, we require to specify a prior distribution µ 0 (u) = P(u) of the unknown, before the data are collected. This distribution comprises all our prior knowledge of the unknown and may include, for example, the regularity of the space of admissible solutions to the inverse problem. For the present work we consider Gaussian priors which have been used to characterize the uncertainty in the (log) permeability that arises from the preform fabrication [49,30,29] (see also references therein). In particular, here we consider stationary Gaussian priors µ 0 = N (u, C) with covariance operator C that arises from the Wittle-Matern correlation function defined by [27,24,38,42]: where Γ is the gamma function, l is the characteristic length scale, σ 2 0 is an amplitude scale and K ν is the modified Bessel function of the second kind of order ν. The parameter ν controls the regularity of the samples. It can be shown [11,42] that, for any ν > 0, if u ∼ µ 0 , then u ∈ C[0, x * ] almost-surely, i.e. µ 0 ([0, x * ]) = 1. This requirement, together with the continuity of the forward map ensures the well-posedness of the Bayesian inverse problems as we discuss in the next subsection. In the context of composite preform's permeability, it is natural to choose the mean u according to the log-permeability intended by the design of the composite part [37].
For computational purposes we use the prior to parametrize the unknown u in terms of its Karhunen-Loeve (KL) expansion [2]: with coefficients u k and where λ k and v k are the eigenvectors and eigenfunctions of C, respectively. A random draw from the prior u ∼ N (u, C) can then be obtained from (2.16) with drawing u k ∼ N (0, 1) i.i.d.

The
Posterior. From (2.13) and our Gaussian assumptions on the observational noise, it follows that for a fixed u ∈ X, we have y n = G n (u) + η n ∼ N (G n (u), Γ n ). Therefore, the likelihood of y n |u is given by (2.17) l n (u, y n ) ∝ exp − At a given time t = t n , the Bayesian posterior µ n (u) = P(u|y 1 , y 2 , . . . , y n ) is defined by the following infinite-dimensional version of Bayes's rule. [45]). Let {G s } N s=1 be the sequence of forward maps defined by (2.10) and let {l s (u; y s )} N s=1 be the corresponding likelihood functions (2.17). Let µ 0 = N (u, C) be the prior distribution with correlation function (2.15). Then, for each n ∈ {1, . . . , N }, the conditional distribution of u|y 1 , · · · , y n , denoted by µ n , exists. Moreover, µ n µ 0 with the Radon-Nikodym derivative

Theorem 2.3 (Bayes Theorem
Proof: The proof follows from the application of Theorem 6.31 in [45] and the continuity of the forward maps (Theorem 2.2) on a full µ 0 -measure set X.
Note that from our assumption of independence of η 1 , . . . , η n , the right hand side of (2.18) is the likelihood of y 1 , . . . , y n |u. Remark 2.3. Due to the assumption of independence between front location and pressure measurements, the likelihood (2.17) can be expressed as This enables us to define two particular cases of the inverse problem. The first case corresponds to the assimilation of only pressure measurements y p n , while in the second case only front location measurements y Υ n are assimilated. Similar arguments to those that led to Theorem 2.3 can be applied (with l β s (u, y β n ) instead of l s (u, y n )) to define the Bayesian posteriors µ p n and µ Υ n associated to these two Bayesian inverse problems. In Section 5.1.1 we will study these two particular cases with an eye towards understanding the added value of assimilating observations of the front location with respect to assimilating only pressure measurements.
3. Approximating the posteriors via Sequential Monte Carlo method. In the previous section we have established the well-posedness of the Bayesian inverse problem associated to inferring the logpermeability in the one-dimensional moving boundary problem (2.1)-(2.5). The solution of this inverse problem is the sequence of posterior measures {µ n } N n=1 defined by Theorem 2.3. As we discussed in Section 1, these posteriors cannot be expressed analytically and so a sampling approach is then required to compute the corresponding approximations. Note that the sampling of each posterior µ n (n = 1, . . . , N ) can be performed independently by, for example, Markov chain Monte Carlo (MCMC) methods. However, we reiterate that, for the present application SMC samplers are rather convenient as they exploit the sequential nature of the considered inverse problem by enabling a recursive approximation of the posterior measures as new data (in time) become available. Such recursive approximations of the posterior could enable practitioners to update their probabilistic knowledge of preform's log-permeability which is, in turn, essential to develop real-time optimal control strategies for RTM under the presence of uncertainty.
Recognizing that the inverse problem under consideration involves inferring a function potentially discretized on a very fine grid, it is vital to consider the application of SMC samplers such as the one introduced in [21], carefully designed for approximating measures defined on a high-dimensional space. In this section we review and apply this scheme for the approximation of the Bayesian posteriors {µ n } N n=1 that we defined in the previous section. The aims of this section are to (i) provide a deeper quantitative understanding of the accuracy of the fully-Bayesian methodology of [21] with respect to its computational cost under practical computational conditions; (ii) provide a motivation for the proposed REnKA that we propose from this SMC sampler in Section 4; and (iii) define accurate approximations of {µ n } N n=1 which we use as a benchmark for testing our REnKA scheme.
In Section 3.1 we briefly discuss the essence of the standard SMC that we then use in Sections 3.2-3.3 to review methodological aspects of the adaptive-tempering SMC sampler for high-dimensional inverse problems of [21]. We then apply this SMC in Section 3.4 for the solution of the Bayesian inverse problem in the 1D case defined in the previous section. In Section 3.5 we assess practical limitations of the SMC.
3.1. Standard SMC for Bayesian inference. As we discussed in the Introduction, starting with the prior µ 0 , the objective of SMC is to recursively compute an approximation of the sequence of Bayesian posteriors {µ n } N n=1 in terms of weighted particles. More specifically, assume that at the observation time t n , we have a set of J particles {u , which provides the following particle approximation of µ n−1 (u) = P(u|y 1 , . . . , y n−1 ): The objective now is to construct a particle approximation of µ n (u) = P(u|y 1 , . . . , y n ), which includes the new data y n collected at time t n . In a standard SMC framework [32,8,10], this particle approximation is constructed by means of an importance sampling step with proposal distribution µ n−1 . To illustrate this methodology, let us first note formally that where we have used which can be obtained directly from Theorem 2.3. An approximation of (3.2) can be obtained by From (3.4) we see that the importance (normalized) weights W (j) n assigned to each particle u (j) n−1 define the following empirical (particle) approximation of µ n : However, the accuracy of such empirical approximation relies on µ n−1 being sufficiently close to µ n ; when this is not the case, after a few iterations (observation times) the algorithm may produce only a few particles with nonzero weights. This is a well-known issue of weight degeneracy that often arises from the application of empirical (importance sampling) approximations within the context SMC samplers [5]. Weight degeneracy is routinely measured in terms of the Effective Sample Size (ESS) statistic [23]: which takes a value between 1 and J; ESS = J when all weights are equal and ESS = 1 when the distribution is concentrated at one single particle. A common approach to alleviate weight degeneracy is, for example, to specify a threshold for the ESS below which resampling (often multinomially) according to the weights {W (j) n } J j=1 is performed. Resampling discards particles with low weights by replacing them with several copies of particles with higher weights. The approximation of a sequence of measures via the combination of the importance sampling step followed with resampling leads to the Sequential Importance Resampling (SIR) scheme [10].
It is important to note that the aforementioned resampling step in SIR can clearly lead to the lack of diversity in the population of resampled particles. This is, in turn, detrimental to the approximation of the sequence of posteriors. The general aim of the standard SMC approach is to diversify these particles by a mutation step with involves replacing them with samples from a Markov kernel K n with invariant distribution µ n . In the following subsection we provide a discussion of the aforementioned mutation in the context of the SMC sampler for high-dimensional inverse problem [21]. We refer the reader to [33,32,8,10] for a thorough treatment of more standard SMC samplers.
3.2. SMC for high-dimensional inverse problems. The weight degeneracy in the importance sampling step described above is more pronounced when the two consecutive measures µ n−1 and µ n differ substantially from each other. This has been particularly associated with complex (e.g. multimodal) measures defined in high-dimensional spaces. When the change from µ n−1 to µ n is abrupt, the importance sampling step can result in a sharp failure, whereby the approximation of µ n is concentrated on a single particle [5]. Recent work for high-dimensional inference problems has suggested [21,1] that further stabilization of the importance weights is needed by defining a smooth transition between µ n−1 and µ n . For the present work, we consider the annealing approach of [34,33], where q n intermediate artificial measures {µ n,r } qn r=0 are defined such that µ n,0 = µ n−1 and µ n,qn = µ n . These measures can be bridged by introducing a set of q n tempering parameters denoted by {φ n,r } qn r=1 that satisfy 0 = φ n,0 < φ n,1 < · · · < φ n,qn = 1 and defining each µ n,r as the probability measure with density proportional to l n (u, y n ) φn,r with respect to µ n−1 . More specifically, µ n,r satisfies Note that when q n = 1, φ n,1 − φ n,0 = 1 and so expression (3.9) reduces to (3.3). We now follow the SMC algorithm for high-dimensional inverse problems as described in [21].

Selection
Step. The first stage of the SMC approach of [21] is a selection step which consists of careful selection of the tempering parameters which define the intermediate measures {µ n,r } qn r=0 ; these are in turn approximated by the application of the SIR scheme described above. Let us then assume that at an observation time t n and iteration level r −1, the tempering parameter φ n,r−1 has been specified, and that a set of particles u (j) n,r−1 provides the following approximation (with equal weights) of the intermediate measure µ n,r−1 : From (3.9) we can see that the new tempering parameter φ n,r must be selected to ensure that φ n,r −φ n,r−1 is sufficiently small, so that the subsequent measure µ n,r is close to µ n,r−1 thus preventing a sharp failure of the empirical approximation of µ n,r (3.6). In particular, once the next tempering parameter φ n,r is specified, we note from expression (3.9) that the importance weights for the approximation of µ n,r are given by Recognizing that the ESS in (3.7) quantifies weight degeneracy in SIR, the approach of [21] (see also [20]) proposes to define on-the-fly the next tempering parameter φ n,r by imposing a fixed, user-defined value J thres on the ESS. More specifically, φ n,r is defined by the solution to the following equation: which may, in turn, be solved by a simple bisection algorithm on the interval (φ n,r−1 , 1]. An approximation of µ n,r is then given by the weighted particle set {u n,r } J j=1 . If at the r − 1 level, we find that ESS n,r (1) > J thresh , it implies that no further tempering is required and thus one can simply define φ n,r = 1. We note that the number of tempering steps q n is random. While the tempering approach described above is aimed at preventing ESS from falling below a specified threshold J thres and thus avoiding a sharp failure of the empirical approximation of µ n,r , resampling is still required to discard particles with very low weights. Let us then denote byû (j) n,r (j = 1, . . . , J) the particles, with equal weights, that result from resampling with replacement of the set of particles u (j) n,r−1 according to the weights W (j) n,r . 3.2.2. Mutation Phase. As stated in the preceding subsection, at the core of the SMC methodology is a mutation phase that adds diversity to the population of the resampled particlesû (j) n,r . In the context of the tempering approach described above, this mutation is conducted by means of sampling from a Markov kernel K n,r with invariant distribution µ n,r . Similar to the approach of [21], here we consider mutations given by running N µ steps of an MCMC algorithm with µ n,r as its target distribution. More specifically, we consider the preconditioned Crank-Nicolson (pcn)-MCMC method from [9] with target distribution µ n,r and reference measure µ 0 . Formally, these two measures are related by dµ n,r dµ 0 ∝ l n,r (u, y n ) ≡ l n (u, y n ) φn,r n−1 s=1 l s (u, y s ).
The pcn-MCMC method for sapling µ n,r is summarised in Algorithm B.1 (see Appendix B). Under reasonable assumptions this algorithm produces a µ n,r -invariant Markov kernel [21]. The resulting particles denoted by u n,r , ·)) then provide the following particle approximation of µ n,r : where the convergence is proven in a suitable metric for measures [1]. Note that at the end of the iteration r = q n , the corresponding particle approximation µ J n = µ J n,qn ≡ 1 n,qn provides the desired approximation of the posterior that arises from the Bayesian inverse problem of interest. This SMC sampler is summarized in Algorithm B.2 (see Appendix B).
Remark 3.1. For simplicity, here we use the resampling step at every iteration of the SMC sampler. However, whenever ESS n,r (1) > J thresh (and so φ n,r = 1) the resampling step can be skipped; this involves using the corresponding weighted particle approximation and modifying the formula for the incremental weights as discussed in [21, Section 4.3].

A note on tempering.
Let us define the following inverse of the increment in tempering parameters: and note that 0 ≤ φ n,r ≤ 1 implies α n,r ≥ 1. In addition, expression (3.9) can be written as where we have used the definition of the likelihood in (2.17). Informally, we can then interpret each iteration of the SMC sampler (at a given observation time t n ) as the solution of a Bayesian inverse problem that consists of finding µ n,r given the prior µ n,r−1 and the data: From (3.15) and the fact that 0 ≤ φ n,r ≤ 1, it follows that α n,r ≥ 1. Therefore, (3.17) is nothing but the original problem (2.13) albeit with a noiseη n,r that has an inflated covariance α n,r Γ n .
We also note that α n,r plays the role of a regularization parameter in the sense that it controls the transition between µ n,r−1 and µ n,r . The larger the α n,r the smoother this transition. Alternatively, we can see that α n,r can be interpreted as a "temperature" in the tempering scheme which, in turn, flattens out the likelihood function at the observation time t n . Clearly, more tempering will be required whenever ||(Γ n ) −1/2 (y n − G n (u))|| 2 is large; this can for example happen if the observational data are accurate (i.e small Γ n ) and/or many observations are available.
The amount of tempering is controlled by the number of parameters obtained via (3.12). The greater the number of tempering parameters, the larger the α n,r 's which in turn indicates that more regularization is needed to ensure a stable transition between those measures. This has also, in turn, an associated increase in iterations and thus in computational cost.
3.3.1. Computational aspects of SMC. The main computational cost of the SMC sampler previously discussed is attributed to the mutation step for which N µ steps of the pcn-MCMC algorithm are performed. At each observation time t n and iteration r, the SMC sampler then requires J N µ evaluations of the nth forward map G n . Therefore, the computational cost of computing µ n is q n g n J N µ , where g n denotes the computational cost of evaluating G n which, in turn, corresponds to solving the moving boundary problem from time t = 0 up to time t n . The total computational cost of computing the full sequence of posteriors {µ n } N n=1 is then which is expressed in terms of g N , the cost of evaluating G N (i.e. solving the forward model from time zero up to the final observation time).
The work of [21] has suggested that accurate approximations of the posterior via SMC samplers require, for example, values of N µ = 20 and J = 10 4 . If we assume for a moment that only one observation time N = 1 is considered and that only one tempering step q 1 = 1 is required to compute µ 1 , the computational cost in this case would be approximately 10 5 times the cost of solving the forward model from time t = 0 up to time t 1 . Such cost would be clearly computationally prohibitive for practical applications, where the aforementioned forward simulation may take several minutes of CPU time. In particular, for the 2D or 3D version of the RTM process, the high computational cost of the SMC sampler becomes impractical. While reducing the values of J and N µ may result in a more affordable computational cost, this is substantially detrimental to the level of accuracy of the SMC sampler as we show via numerical experiments in Section 3.5.

Numerical examples with SMC.
In this subsection we report the results from the numerical application of the SMC sampler discussed in the previous subsection. The objective is to approximate the sequence of Bayesian posteriors that arise from the 1D moving boundary problem defined in Section 2 for the experimental set-up described in Section 3.4.1. In Section 3.4.2 we discuss the numerical results obtained via the SMC sampler with a very high number of particles which results in accurate approximations of the Bayesian posteriors. These approximations are then used in Section 3.5 to assess the practical limitations of the scheme under certain choices of tunable parameters and number of particles. These limitations motivate the approximate methods that we propose in Section 4.

Experimental set-up.
We consider a dimensionless version of the one-dimensional model (2.1)-(2.5) which together with its numerical approximation is described in Appendix C. The dimensionless values for the control variables are p 0 = 1 and p I = 2. We use a Gaussian prior distribution µ 0 = N (u, C) with the covariance operator C that arises from the covariance function defined in (2.15). We numerically solve (off-line) the eigenvalue problem associated to the matrix that results from discretizing C; the corresponding eigenvector/eigenvalues are then stored for subsequent use in the parameterization of the log-permeability in the SMC sampler. The KL expansion (2.16) becomes a truncated sum with a number of elements equal to the the total number of eigenvalues of this matrix; these are, in turn, equal to the number of cells used for the discretization of the domain D * = [0, 1]. No further truncation to this KL expansion is carried out. A few samples from the prior are displayed in In order to generate synthetic data, we define the "true/reference" log permeability field u † whose graph (red curve) is displayed in Figure 3.2 (top-left); this function is a random draw from the prior described above. We use u = u † in the numerical implementation of (2.1)-(2.5) in order to compute the true pressure field p † (x, t) as well as the true front location Υ † (t). The plot of p † (x, t) is shown in where η p n and η Γ n are Gaussian noise (see subsection 2.2) with standard deviations equal to 1.5% of the size of the noise-free observations. Synthetic pressure data {y p n } 5 n=1 are superimposed on the graphs of {p(x, t n )} 5 n=1 in Figure 3.1 (middle). Synthetic front locations {y Υ n } 5 n=1 are 0.21, 0.39, 0.59, 0.74, 0.86. In order to avoid inverse crimes, synthetic data are generated by using a finer discretization (with 120 cells) than the one used to approximate the posteriors (with 60 cells).

Application of SMC.
In this subsection we report the application of the SMC sampler of [21] (see Algorithm B.2 in Appendix B) which, as described in the preceding section, provides a particle approximation of each posterior that converges to the exact posterior measure µ n as the number of particles J goes to infinity. In order to achieve a high-level of accuracy we use J = 10 5 number of particles which is substantially larger compared to the number of particles (e.g. 10 3 to 10 4 ) often used in existing applications of SMC for high-dimensional inverse problems [8,21]. In addition, we consider the selection of tunable parameters N µ = 20 and J thresh = J/3 similar to the ones suggested in [21]. For each observation time t n , we store the ensemble of particles {u (j) n } 10 5 j=1 that approximates the corresponding posterior µ n . From this ensemble, we compute the 0.02, 0.25, 0.5, 0.75, 0.98 posterior percentiles displayed in Figure 3.2 (top-middle to bottom-right), where we also include the graph of the true log-permeability (red curve). The vertical line in these figures indicate the true location of the front Υ † (t n ) at each observation time t n . We can clearly appreciate that the uncertainty band defined by these percentiles is substantially reduced as more observations (in time) are assimilated. In fact, the main reduction of the uncertainty is observed in the region of the moving domain D † (t n ) = [0, Υ † (t n )] at the corresponding observation time t n . It is then clear that at each observation time t n , measurements collected from pressure sensors with x m ∈ D * \ D † (t n ) are not very informative of the log-permeability field. This comes as no surprise when we recognize that the pressure field given by (2.8) depends on the permeability field only in the region of the moving domain D(t). In other words, the values of the permeability in the region defined by D * \ D(t n ) have no effect on p(x, t n ); hence the nth likelihood function is independent of u in this region. We can indeed observe from Figure 3.2 that the percentiles of the log-permeability in this region (see domain to the right of the vertical lines) is similar to those from the prior. However, due to the regularity of the log-permeability enforced in the prior µ 0 , there is a smooth transition in the uncertainty band at the interface defined by the front location Υ(t n ).
The number of intermediate tempering distributions that SMC adaptively computed to approximate the sequence of posteriors {µ n } 5 n=1 were the following: We use these numbers in (3.18) to compute the total computational cost of approximating the sequence {µ n } 5 n=1 . The values of g n (i.e. cost of evaluating each G n ) are estimated by the average CPU time from 1000 simulations computed with different log-permeabilities sampled from the prior. We obtain that total cost is approximately 1.5 × 10 7 times the cost of evaluation the 5th forward map G 5 (i.e. at the final observation time). Clearly, this computational cost is prohibitive for the two and three dimensional problems where, as stated earlier, evaluating the forward map can take several minutes of CPU time. For the present one-dimensional case we are able to afford this cost due to the relatively low cost associated with solving the 1D moving boundary problem. 3.5. Reducing the cost of SMC by adjusting tunable parameters. Given the high computational cost of computing accurate approximations of the posteriors with SMC, it is reasonable to ask whether its computational cost can be reduced by adjusting the tunable parameters in (3.18). By reducing either the number of particles J and/or the number of MCMC steps N µ , we can achieve a substantial decrease in the computational cost. The selection of J thresh also determines the computational cost as it, in turns, defines the number of tempering steps for each posterior. However, it is essential to understand the effect of decreasing these tunable parameters on accuracy of the SMC sampler. In this subsection we aim at understanding this effect by comparing the application of the SMC sampler with smaller number of particles J and different choices of the tunable parameters N µ and J thresh . This requires creating a Benchmark against which we can compare performance of SMC. The Benchmark is obtained by the highly-resolved characterization of the posteriors that we computed in the preceding section by using the SMC sampler with large number of particles (J = 10 5 ). In Appendix B we provide further discussions of the performance and diagnostics of the SMC sampler applied to approximate these posterior measures.
These diagnostics offer evidence that the SMC sampler has been successfully applied, thereby providing accurate characterization of the posterior that we may use as a Benchmark to compare against the posteriors computed via algorithms with lower resolution/accuracy. The numerical investigation below is aimed at assessing SMC with different selections of ensemble size J as well as the tunable parameters N µ and J thresh .
For the reasons stated above, through the rest of the this and the following sections, we refer to the aforementioned highly-resolved SMC particle approximations (with J = 10 5 ) as the "exact" sequence of posteriors {µ n } 5 n=1 that we use for subsequent comparisons purposes. Moreover, for these comparisons we assume that the sample mean and variance of these SMC samples are exact approximations of the mean E µn and variance V µn of the posterior µ n . In other words, we assume E µn = u n,10 5 , V µn = σ 2 n,10 5 , Let us now consider the application of the SMC sampler for the following choices of small number of particles: J = 50, 100, 200, 400, 800, 1600, 3200, 6400. We also consider three choices of the tunable parameter J thresh (J thresh = J/3, J/2, 2J/3) and two choices of N µ (N µ = 5, 20). In Figure 3.3 we show percentiles of the log-permeability posteriors µ n (for n = 1, 3 and 5) obtained using the aforementioned SMC sampler for some of those choices of the number particles J, and with the same selection of tunable parameters N µ = 20, J thresh = J/3 that we used for the highly-resolved SMC with large particles; percentiles from the latter are included in the right column of Figure 3.3 for comparison purposes. We can see that as the ensemble of particle increases, the approximation of SMC improves when compared to the one provided by the highly-resolved SMC. Note that very small number of particles results in very poor approximations of these percentiles. In order to quantify the level of approximation obtained with SMC with the aforementioned selections of parameters, we compute the L 2 (D * )-relative errors of the mean and variance with respect to the posterior measure approximated with the highly-resolved SMC computed as described in the preceding subsection. More precisely, let us define  Quantities E J n , V J n and n J are random variables that depend on the initial ensemble of particles that we generate from the prior µ 0 . We thus report these quantities (for each n = 1, 3, 5) averaged over 15 experiments corresponding to different selections of the initial ensemble of particles. In Figure 3.4 we display E J n (top), V J n (middle) and J n (bottom) for (from left to right) n = 1, 3, 5 as a function of the aforementioned selections of ensemble size J. For brevity we omit the results for n = 2, 4 as they display similar behaviour. The total computational cost of computing the full sequence of posteriors (i.e. C SM C from (3.18)) is shown in Figure 3.5 (left). We reiterate that this cost is expressed in terms of the number of evaluations of the 5th forward map G 5 .
While the numerical analysis of the convergence of the SMC sampler is beyond the scope of this work, the results presented in this section are aimed at understanding the level of accuracy of SMC with relatively small number of particles and for a selection of tunable parameters which may enable the use of this method in more practical scenarios. From these results it is clear that the selections of J thresh have no substantial effect on the accuracy of the scheme in terms of approximating the mean and variance of each posterior. Similarly, the computational cost with respect to our selections of J thresh does not seem to vary significantly. It is evident that the main effect in terms of accuracy is the choice of MCMC steps (i.e. parameter N µ ). Indeed, note that the error obtained with N µ = 5 is considerably larger than the one with N µ = 20 although the computational cost of the former is one quarter of the computational cost of the latter. We conclude that even though decreasing N µ can offer computational affordability, it is detrimental to the approximation properties of the scheme. This comes as no surprise as it is well known that the mutation step that involves running MCMC is crucial for the accuracy of any SMC methodology.
The behavior of the SMC sampler with respect to the number of particles J is as expected. On the one hand, an increase in J corresponds to a decrease in the error with respect to the mean and variance. On the other hand, the computational cost, C SM C , increases with J. Note that there is a clear linear relationship between these two variables which is, in turn, obvious from (3.18) provided that q n is invariant with respect to J. Indeed, for the cases considered here, the number of intermediate tempering distributions (not reported) computed at each observation time, is invariant with respect to our choices of J. This is somewhat an expected outcome since our choice of J tresh in (3.12) is always a fraction of J. It is also worth mentioning that the effect of J is less noticeable when we look at the error with respect to the truth. At each observation time, we notice that the n seems to converge to a nonzero value as J increases. Note that convergence to the truth is not ensured due to the limited number of measurements inverted and the potential lack of identifiability of the log-permeability.
The results reported in this subsection suggest that achieving a reduction in the computational cost by reducing N µ has a severe detrimental effect in the accuracy of the SMC sampler with small number of particles. In addition, J thresh does not seem to have a substantial effect in either the accuracy or computational cost. Clearly, we are only then limited to the number of particles J to control the computational cost of the sampler without severely compromising accuracy of the approximate posteriors. In this section we propose a regularizing ensemble Kalman algorithm (REnKA) that aims at providing an accurate approximation of the sequence of Bayesian posteriors at a much lower computational cost. In Section 4.1 we introduce REnKA as a Gaussian approximation from the SMC sampler of [21] discussed in the preceding section. The proposed REnKA in the context of existing ensemble Kalman methods is discussed in Section 4.2. A numerical investigation of the convergence properties of REnKA is reported in Section 4.3.
For the subsequent development of the proposed scheme, we extend the domain of definition of the sequence of forward maps G n introduced in (2.10). More specifically, we assume G n : X → R M +1 where X is a Hilbert space such that X = C[0, x * ] → X (compactly). We denote by < ·, · > X and < ·, · > the inner products in X and R M +1 , respectively. In addition, we define Z ≡ X × R M +1 with inner product denoted by < ·, · > Z . 4.1. Motivation for REnKA. Motivated by the SMC tempering approach described in the previous section, we now propose an ensemble Kalman algorithm whose aim is to approximate {µ n,r } qn r=1 by a sequence of Gaussian measures {ν n,r } qn r=1 which are, in turn, characterised by a set of particles with equal weights. Suppose that, at time t = t n we have an ensemble {u (j) n,r−1 } J j=1 of J samples from a Gaussian measure ν n,r−1 that approximates µ n,r−1 , and a prescribed tempering parameter φ n,r−1 . We may then solve (3.12) for the new φ n,r and define the regularization parameter α n,r in (3.15). We now wish to make a transition from ν n,r−1 to a Gaussian measure ν n,r that approximates µ n,r . To this end, let us define the augmented variable z = (u, G n (u)) T ∈ Z (4.1) and note that, in terms of this variable, we may rewrite (3.17) as y n = Hz +η n ,η n,r ∼ N (0, α n,r Γ n ), (4.2) where H = (0, I) and I is the identity operator. One can see that by reformulating the inverse problem in terms of the augmented variable, the resulting forward map (i.e. H) acting on this variable is linear.
From (4.1) we define the following augmented particles   The Gaussian measureν n,r−1 (z) is used to approximate the measure, denoted byμ n,r−1 (z), that arises from pushing forward µ n,r−1 (u) under (4.1). By using this Gaussian approximation ofμ n,r−1 (z), we then provide a Bayesian formulation of the inverse problem given by (4.2). More specifically, we wish to computeν n,r (z) ≡ P(z|y n ) givenν n,r−1 (z) and the data from (4.2). A formal application of Bayes theorem yields dν n,r dν n,r−1 (z) ∝ exp − 1 2 ||(α n,r Γ n ) −1/2 (y n − Hz)|| 2 . Moreover, from (4.4) and the linearity of the forward map H (on the augmented variable), it follows by standard arguments [25] that ν n,r (z) = N z n,r−1 + K n,r (y n − Hz n,r−1 ), (I − K n,r H)C n,r−1 , (4.8) where K n,r ≡ C n,r−1 H T (HC n,r−1 H T + α n,r Γ) −1 . (4.9) Let us then note that C n,r−1 in (4.6) can be written as  Informally, we use the block structure of (4.10) and define ν n,r , the approximation of µ n,r , as the marginal measure ofν n,r given by ν n,r (u) ≡ N u n,r−1 + K u n,r (y n − G n,r−1 ), C uu n,r−1 − K u n,r (C uw n,r−1 ) T , (4.14) where K u n,r = C uw n,r−1 (C ww n,r−1 + α n,r Γ) −1 . Although the measure (4.14) is fully characterised by its mean and covariance, for the subsequent tempering step we need a particle approximation of ν n,r (u). We can obtain those particles by updating the current set of particles u (j) n,r−1 via the formula n,r−1 + C uw n,r−1 (C ww n,r−1 + α n,r Γ) −1 (y (j) n,r − G n (u (j) n,r−1 )), (4.16) where y (j) n,r ≡ y n + η (j) n,r , η (j) n,r ∼ N (0, α n,r Γ n ). (4.17) Indeed, under the standard assumption that the noise η (j) n,r is independent from the variables u (j) n,r−1 and G n (u (j) n,r−1 ), it can be shown by the standard arguments in Kalman-based methods (see for example [7]) that n,r (z) → ν n,r (u), J → ∞. Expression (4.16) and selection of the regularisation parameter α n,r based on the adaptive tempering approach discussed in Section 3.2 constitutes the proposed scheme summarized in Algorithm 4.1.
Remark 4.1. Note that the key assumption for the proposed scheme is the Gaussian approximation ofμ n,r−1 (z) provided by (4.8). It is clear that the measureμ n,r−1 (z) is, as a rule, non-Gaussian and the aforementioned assumption will result in a methodology that will, in general, not converge to the posteriors µ n as the ensemble size J → ∞. Nevertheless, we will show via numerical examples that this approximation provides reasonably accurate estimates using only a small number of particles.
It is not difficult to see that the main computational cost of REnKA, in terms of the cost of evaluating the forward model at the final observation time, is given by where, as before, g n denotes the computational cost of evaluating the G n -forward map. As we will demonstrate via numerical experiments, for the moving boundary problem of Section 2.1, REnKA offers a computationally affordable and thus practical approach to approximate the solution to the Bayesian inverse problem that arises from RTM. n,r−1 ) needed below). Compute the tempering parameter φ n,r : if min φ∈(φ n,r−1 ,1] ESS n,r (φ) > J thresh then set φ n,r = 1. else compute φ n,r such that ESS n,r (φ) ≈ J thresh using a bisection algorithm on (φ n,r−1 , 1]. end if Construct C uw n,r−1 , C ww n,r−1 defined by expressions (4.11)-(4.12). Update each ensemble member: n,r−1 )), (4.19) where α n,r = (φ n,r − φ n,r−1 ) −1 , y (j) n,r = y n + η (j) n,r , with η (j) n,r ∼ N (0, α n,r Γ n ). end for end while Set u

end for
It is important to mention that, at a given observation time t n and iteration level r, the value of J j=1 l n (u (j) n,r−1 , y n ) 1−φ n,r−1 may be zero to machine precision. In this case, the tempering parameter φ n,r is not be computable, via a bisection scheme on (φ n,r−1 , 1], as stated in Algorithm 4.1. This computational issue is more likely to arise at the early iterations of the scheme for which the value of φ n,r−1 is not sufficiently close to one. This can be overcome, for example, by simply adapting the bisection algorithm in order to first compute a φ * such that J j=1 l n (u (j) n,r−1 , y n ) φ * −φ n,r−1 > 0. If min φ∈(φ n,r−1 ,φ * ] ESS n,r (φ) > J thresh we then set φ n,r = φ * ; otherwise, we find φ n,r by solving (3.12) via a bisection algorithm on (φ n,r−1 , φ * ]. For the numerical experiments reported in the present work, zero values to machine precision for J j=1 l n (u (j) n,r−1 , y n ) 1−φ n,r−1 were only encountered where a large number of measurements were inverted in the 2D setting of Section 5.2.

REnKA in the context of existing ensemble Kalman methods for inverse problems. Ensemble
Kalman methods for inverse/calibration problems have been widely used in the last decades [15]. More recently, using iterative Kalman methods with a regularization parameter (e.g. α n,r in (4.19)) have been proposed for a wide class of applications. In particular, the proposed REnKA scheme can be related to the recently developed regularizing ensemble Kalman method introduced in [19] for solving classical (deterministic) inverse problems. More precisely, Algorithm 4.1 is nothing but a sequential version of the iterative scheme presented in [19] except for the selection of the regularization parameter α n,r . While in the present work we have motivated Algorithm 4.1 from the SMC framework, the algorithm in [19] was obtained as an ensemble approximation of the regularizing Levenberg-Marquardt scheme initially developed in [16] for solving ill-posed inverse problems. In the context of the proposed work, REnKA aims at providing a derivative-free approximation to the solution of ||Γ −1/2 n (y n − G n (u))|| → min (4.21) in a stable (regularized) fashion. The regularization parameter in [19] was selected according to the discrepancy principle in order to regularize the inverse problem posed by (4.21). In contrast, the present work uses the adaptive tempering approach of [21] for the selection of this regularization parameter in the context of SMC. It is clear that tempering can be understood as a regularization to the Bayesian inverse problem; the effect of α n,r is to flatten out the posterior and allow for a more controlled/regularized transition between the sequence of measures. Other works highlighting the connection between ensemble Kalman methods and SMC approaches include [43,44,39]. In addition, by noticing from (3.15) that, for each n = 1, . . . , N , qn r=1 α −1 n,r = 1, the proposed REnKA can be also understood as a sequential version of the ensemble smoother with multiple data assimilation proposed by [12]. However, it is important to reiterate that the adaptive selection of α n,r proposed here is inherited from the SMC approach of [21]. This selection differs substantially from the strategy proposed in [12] for which the number of intermediate tempering distributions q n is fixed and selected a priori. We can clearly observe that the approximations provided by REnKA improves as the ensemble size J increases. More importantly, we can note that the uncertainty band defined by these percentiles provided better approximations than those from SMC with the same number of particles (see Figure 3.3).
We quantify the level of accuracy of REnKA with respect to the Benchmark defined by the highlyresolved SMC sampler reported in Section 3.4.2. In ≡ 24%), we need to apply the SMC sampler (say with J thresh = J/3, and N µ = 20) with J = 6400 particles for which the computational cost is approximately 5 × 10 5 G 5 -forward model evaluations.
The results above not only demonstrate that, when a small number of particles is used, the per-formance (accuracy vs computational cost) of REnKA outperforms SMC, but also these results show that REnKA is robust for reasonable selections of the tunable parameter J thresh . Similar to SMC, this parameter determines the number of tempering distributions at each observation time and thus has an impact on the computational cost of the scheme. It is also important to remark that even though the relative errors of REnKA with respect to the mean and variance decrease as the ensemble size increases, these errors seem to converge to a non-zero value thereby indicating that REnKA does not provide an asymptotic convergence to the posterior measures as J → ∞. Nevertheless, our results clearly showcase the advantages of REnKA for approximating these measures in an accurate and efficient fashion for a limited and realistic computational cost.

Numerical investigations of the Bayesian Inverse problem in RTM.
Having numerical evidence from Section 4.3 that REnKA is a robust and computationally feasible approach for addressing the Bayesian inverse problem defined in Section 2, in this section we use REnKA to provide further practical insights in the RTM Bayesian inverse problem. Section 5.1 is devoted to the 1D case studied earlier which includes the study of Section 5.1.1, concerning the effect of the number/type of measurements on the solution to the Bayesian inverse problem. The effect of the number of observation times as well as the observational noise level are investigated in Sections 5.1.2 and 5.1.3, respectively. The application of REnKA for the 2D forward model stated in Section 1 is reported in Section 5.2. For the 1D (resp. 2D) case we consider REnKA with a fixed number of J = 1000 (resp. J = 150). In all cases we select the tunable parameter J thresh = J/3. For clarity of the notation, in the expression for the ensemble mean and variance (3.21), we then omit the index J as appropriate. We will focus our numerical investigations in terms of the following quantities defined at each observation time t n (n = 1, . . . , N ): (A) The L 2 (D * )-relative error with respect to the truth, n defined in (3.23) in terms of the ensemble

Further investigations of the 1D case.
In this section we study how the number of observation times and the observational noise level influence the solution to the Bayesian inverse problem in the 1D case.

Effect of the number and type of measurements.
Since the number of pressure sensors available for real RTM processes is usually limited, in this section we study the effect of the number of pressure measurement locations in the solution to the Bayesian inverse problem. In addition, we wish to understand the added value of front location measurements in reducing uncertainty characterized by the posterior of the log-permeability. Let us then consider the same experimental set-up and relevant model parameters as described in Section 3.4.1. We recall that the pressure measurement configuration consisted of M = 9 sensor locations (see where (resp. X ) indicates that measurements of front are included (resp. excluded) in the application of REnKA.
In Figure 5.2 we show posterior percentiles of the log-permeability for some of the combinations described in (5.3). These posteriors were obtained using the same fixed initial ensemble. We note that the effect of inverting the front (Front: ) in reducing the uncertainty band defined by the posterior percentiles is mainly noticeable at the first observation time at which the front has not yet reached most of the pressure sensors. Therefore, pressure measurement from locations that have not been reached by the front provide little information of the log-permeability. As more observations (in time) are assimilated and more pressure sensors are reached by the front, the effect of inverting the front becomes less significant. We can also observe that, even at earlier observation times, the effect of inverting the front has little effect on reducing this uncertainty band when larger number (M = 20) of pressure sensors are used for the inversion. These results are confirmed in Figure 5.3 (middle) where we display the plot of Σ n . Indeed, Σ n decreases as more (fixed) observations in time are assimilated. Note that at the final observation time t 5 , the added value of the measurement of front location is only noticeable when a small number of pressure sensors M = 5 are considered. In fact, when no pressure data are inverted M = 0, we can note that Σ n is substantially reduced. Therefore, measurements of front location are indeed informative of the log-permeability when only a limited number pressure sensors are available; this also includes the case where pressure measurements are completely absent. However, increasing the values of M has no significant effect on the reduction of Σ n , thus suggesting that adding more pressure sensors will no further reduce the posterior uncertainty in the log-permeability.
The plots of n are displayed in the left panel of Figure 5.3. Although n decreases as more (fixed) observations in time are assimilated, these results neither reveal the added value of the front location nor provide insight in the effect of the number of pressure sensors in terms of improving the accuracy of the infer log-perm (i.e. reducing n ). Indeed, note that for t = t 4 , the results from configuration (M = 20, Front:X ) yield smaller n than (M = 20, Front: ). Furthermore, at t = t 5 we note that M = 9 results in smaller errors than M = 20. It is then clear that, n , as defined in (3.23) does not offer sufficient evidence that increasing the number of pressure sensors results in more accurate inference of the truth. As we recall that the definition of this error (see (3.23)) involves the norm defined on the whole physical domain D * = [0, 1], we should also then reemphasize that, at a given observation time t n , the reduction of the uncertainty in the log-permeability field mainly occurs in the region of the moving domain D † (t n ) = [0, Υ † (t n )], where Υ † is the true front location. Therefore, when applying REnKA (or any Bayesian inversion approach), the estimator of the truth (ensemble mean for REnKA) at time t n , may not necessarily results in decrease of the error in the region D * \ D † (t n ) as no informative measurements have there yet been collected. In order to further understand this effect, let us then consider, Υ n , the error with respect the truth defined by (5.2). This error accounts for error of the estimator w.r.t the truth only in the region D † (t n ), where measurements at time t n (as well as earlier measurements) have been collected. In contrast to n , the error Υ n is defined on the moving domain D † (t n ). Therefore, we shall not necessarily expect a decrease of Υ n as function of the assimilation times t n . In Figure 5.3 (right) we display the log of Υ n from which we can now fully appreciate that more accurate estimates of the permeability (in the domain D † (t n )) are obtained by inverting both pressure and front location as opposed to only inverting pressure data. Moreover, Figure 5.3 (right) now enable us to see that increasing the number of pressure locations, M , yields more accurate estimates of the log-permeability. However, it also reveals that a further increase in the number of pressure measurements locations (e.g. M = 20) has no significant effect in the accuracy of the estimates when compared with smaller and more realistic (from the applicable point of view) measurement configurations (e.g. M = 9). The (log) total computational cost together with the final (log) error Υ n is displayed in Figure 5.4 (left). Computational cost is expressed in terms of forward model evaluations as discussed in subsection 3.5. Note that higher computational cost is associated with larger measurements (i.e. M = 20) as this requires more tempering for the computation of each posterior µ n (see discussion in Section 3.3).    number of total observation times (N = 1, 2, 3,4,5,8,10,12,14,16). For each of these N we select the final time t N = 0.36. In addition, for any two given N 1 and N 2 with N 1 < N 2 we select the observation times so that {t n } N 1 n=1 ⊂ {t n } N 2 n=1 . The space-time measurement configuration for some of these cases is displayed in Figure 5.5 (top) as well as the synthetic pressure data (bottom) generated as before; for these examples we consider the general Bayesian inverse problem where both pressure and front location are assimilated.
The final error w.r.t the truth Υ N as well as Σ N are both displayed in Figure 5.4 (right). These numerical results indicate that increasing the number of observation times does not imply a decrease in either the error w.r.t truth or the uncertainty of the log-permeability at the final observation time.
Although for N > 8 we observe a smaller Υ N , the uncertainty in terms of Σ n is not substantially reduced. In fact, Σ N exhibits a slight increase for N > 8, presumably due to the observational noise of these additional measurements. In Figure 5.6 we display the percentiles of the approximate posterior µ N corresponding to the final observation time t N obtained for N = 1, N = 3, N = 8 and N = 16. While from Figure 5.4 we notice that Σ N increases slightly for N > 8, from Figure 5.6 we still observe a minor decrease in the uncertainty band defined by the posterior percentiles. It is essential to emphasize that even though our experiments suggest that assimilating measurements more frequently does not substantially improve our inference of the truth, the benefit of more frequent observations is the earlier reduction of the uncertainty which can be, in turn, used for the real-time optimization of the RTM process. This is illustrated in Figure 5.7, where we show the percentiles of some of the log-permeability posteriors from one example with a total N = 16 observation times. Indeed, the more frequent the observations, the better the characterization of the log-permeability field in the domain of the true moving front D † (t n ).     as those used in Section 3.4.1. We generate both pressure and front location synthetic data as described in Section 3.4.1 with observational noise levels of 15%, 5%, 1% and 0.5%. Percentiles of the posterior for different noise levels are displayed in Figure 5.8. It comes as no surprise that for smaller noise levels, the posterior uncertainty bands defined by these percentiles are narrower and more concentrated around the true log-permeability. However, we find again that the reduction in the uncertainty is only observed in the region defined by the moving front D * (t n ). Similar to the results from the previous sections, the plot of n (not shown) reveals that, at each observation time, inverting more accurate measurements of pressure and front does not necessarily implies that the log-permeability estimate is more accurate in the regions where the front has not yet arrived. In contrast, the plot of Υ n displayed in Figure 5.9 (left) allows us to clearly observe that smaller observational noise yield more accurate estimates of the truth. In Figure 5.9 (middle) we display Σ n from which we can notice that smaller noise also results in lower uncertainty. Yet, measurements with noise levels below 1.5% have little effect on both the reduction of the uncertainty and the accuracy of the ensemble mean at recovering the truth.
In Figure 5.8 (right) we show the final error with respect to the truth as well as the total computational cost involved in the approximation of the full sequence of posteriors. Smaller noise levels (i.e. more accurate measurements) require more tempering and thus larger number of intermediate measures. This, as discussed in Section 3.3, is reflected in larger computational cost when more accurate measurements are assimilated.

The 2D case.
In this section we apply REnKA for the Bayesian inversion of the 2D moving boundary problem described by (1.1)- (1.9). In contrast to the 1D case, the numerical solution of the   2D moving boundary problem is much more computationally intensive. Therefore, the application of a fully Bayesian methodology such as the SMC sampler considered in Section 3 is impractical for an online computation of the Bayesian posteriors in the 2D case. In this subsection, we exploit the capabilities of REnKA for providing an accurate yet computationally tractable approach for inferring preform (log) permeability alongside with its uncertainty.

Formulation of the 2D Bayesian inverse problem.
Let us consider now the 2D moving boundary problem introduced (1.1)-(1.9) from Section 1. We recall that we are interested in the inference of the log-permeability u(x) = log κ(x) given noisy measurements of the pressure field {p(x m , t n )} M m=1 from M sensors located at points {x m } M m=1 ⊂ D * collected at a discrete observation times {t n } N n=1 . In addition, we wish to invert observations of the front location, or alternatively from the moving domain D(t) that can be potentially obtained from CCT cameras such as in [48,31]. While in the 1D case we can trivially define observations of the (single point) front, in 2D the front Υ(t n ) is a curve which defines the moving domain D(t). Therefore, rather than dealing with measurements of the front Υ(t) itself we may assume pointwise measurements of D(t) via its characteristic function defined by (5.4) χ(x, t) ≡ 1 x ∈ D(t), 0 x / ∈ D(t).
More precisely, we define observations of the form ⊂ D is an array of points for which the characteristic function of D(t) is observed. In practice this array should be considered dense when a high resolution camera is used for capturing the moving domain. Note that this mathematical description of front measurements is suitable as it enables its direct comparison with observations from digitalized images and also with the discrete formulation of (1.2)-(1.9) via the control volume finite element method (CV/FEM) in which the front location is characterized in terms of the filling factor (see details in [37]) rather than using a parameterization of Υ.
For specified (known) p I , p 0 , ∂D I , ∂D N , ϕ, the solution of (1.2)-(1.9) induces a map u = log κ → [p(x, t), D(t)] which enables us to define the sequence of forward maps G n : C(D) → R M +M Υ by means of To our best knowledge, uniqueness, existence and regularity theory for problem (1.2)-(1.9) with non constant κ(x) = e u(x) is an open problem for d > 1 (see a related discussion in [37]). However, for the present work we assume that the following condition holds. We now follow the same formulation of the Bayesian inverse problem as the one we described in Section 2.2. At each observation time t = t n , we collect noisy measurements of the front location y Υ n ∈ R M Υ as well as pressure measurements from sensors y p n ∈ R M . We assume that observations y n = [y Υ n , y p ] T are related to the unknown via expressions (2.13) -(2.14) with G n defined in (5.5). As before, both measurements of D(t) (via its characteristic function) and pressures are assumed to be uncorrelated in time and independent from each other. Furthermore, we consider Gaussian priors µ 0 = N (u, C) with a covariance operator C that arises from the Whittle-Matern correlation function defined in (2.15). The assumption of continuity of the forward maps as well as the fact that µ 0 (C(D * )) = 1 ensures existence of the sequence of posterior measures µ n = P(u|y 1 , . . . , y n ) given by Theorem 2.3 with the same definition of the likelihood functions introduced in (2.17). In the following section we use REnKA to to compute an ensemble approximation of {µ n } N n=1 . 5.2.2. 2D Numerical Experiments. In this subsection we apply REnKA for the solution of the 2D Bayesian inverse problem defined in the previous subsection. The forward model described by (1.1)-(1.9) is solved numerically with the MATLAB code developed in [37] and available from https://github.com/parkmh/MATCVFEM. This code is based on the interface-tracking control volume finite element method (CV/FEM) [3,17,47]. For experiments of this subsection, we consider the following fixed values: Samples from µ 0 are generated via KL parametrization as described in Section 2.2.1 with parameters σ 2 0 = 0.25, ν = 1.5, l = 0.1 and u(x) = 0.0 for all x ∈ D * . Some draws from the prior are displayed in Figure 5.10 (middle row). The log-permeability field is plotted in Figure 5.10 (top-left) is a random draw from the prior that we use as the truth u † for the present experiments.
We where (resp. X ) indicates whether the moving domain D(t n ) has been observed via (5.4). We use the true log-permeability to numerically solve the forward model (1.1)-(1.9) via the CV/FEM code described above. Then, for each of these measurement configurations, synthetic data with a realistic choice of 2.5% Gaussian noise are generated in a similar manner to the one described in Section 3.4.1. In order to avoid inverse crimes, synthetic data are computed on a finer grid than the one we use for the Bayesian inversion via REnKA. These are shown in the middle and right panels of From the plot of Υ displayed in Figure 5.14 (left), we note that, at the latest observation times (n = 6, 7), there is a clear improvement in the accuracy with increasing the number of pressure locations. Similar to the 1D case, we also observe that the benefit of inverting measurements from the moving domain is only noticeable when the number of pressure sensors is relatively small (M = 9). This configuration of sensors is more realistic in practical settings. It is also worth noticing that, at the earliest observation times (n = 1, 2) when the front has not reached most pressure sensors, inverting measurements from the moving domain provides additional information of the log-permeability to the one provided only by pressure measurements.
As more observations (in time) are assimilated, the reduction of the uncertainty in terms of the ensemble variance can be observed from the plot of Σ n displayed in Figure 5.14 (middle). From this plot we also note that the variance decreases as we increase the number of pressure sensors. The added value of measurements from the moving domain is also quite substantial and more noticeable for a small number of pressure sensors. In fact, note that smaller uncertainty has been achieved by inverting only the front (M = 0) compared to the inversion of only pressure data from M = 9 sensors. Here we also find that, at the earliest observation times, the additional inversion of measurements of the moving domain results in further reductions of the uncertainty in comparison to the inversion of only pressure data. While this investigation was conducted with a realistic choice of measurement noise (2.5%), further studies should be conducted to understand the effect of the noise level on the uncertainty estimates of log-permeability in the 2D case.
Finally, the computational cost (see expression (4.20)) of approximating the sequence of posteriors {µ 7 n=1 } via REnKA is displayed in Figure 5.14 (right). This cost is expressed in terms of the number of G 7 forward model evaluations which, in turn, correspond to solving the moving boundary problem from t = 0 to the last observation time t 7 . Furthermore, this cost has been normalised by the number of particles J used in REnKA. This normalisation enables us to provide a rough estimate of the scalable (with respect to J) computational cost of the REnKA (4.20) if each evaluation of the forward map (see step 2(b) in Algorithm 4.1) is conducted in parallel. As discussed in Section 3.3, increasing the number of measurements results in more tempering distributions (i.e. iterations) in the REnKA scheme. Therefore the computational cost increases with the number of measurements. However, for a realistic choice of pressure sensors M = 9, we note that cost of inverting measurements of both front location and pressure sensor results (in average) in a scalable cost of 21 iterations. Since the number of particles that we use for REnKA is relatively low (J = 150), such scalability with respect to the number of particles is reasonable with a high-end computer cluster and can be achieved within a few minutes. 6. Summary and conclusions. In this work we studied the Bayesian inverse problem that arises from inferring physical properties in a setting for porous media flow with a moving boundary. Our investigation is focused on the inference of log-permeability from measurements of pressure and observation of the (moving) domain occupied by resin during the resin injection stage of RTM relevant to the fabrication of composite materials. We adopted the infinite-dimensional Bayesian approach to inverse problems where the aim is to characterise, at each observation time, the posteriors that arise from conditioning the logpermeability on pressure/front measurements. The simplicity of the 1D RTM model enabled us to show existence of the Bayesian posteriors in the aforementioned infinite-dimensional setting. These posteriors were then probed numerically with the dimension-independent SMC sampler for inverse problems from [21]. Our numerical experiments indicated that very large number of particles were needed to accurately approximate the Bayesian posteriors. This resulted in a high computational cost unfeasible for practical RTM settings. In order to reduce computational cost of Bayesian inversions for practical RTM settings, we proposed a regularising ensemble Kalman algorithm (REnKA) that we motivated from the adaptive tempering SMC sampler of [21]. The proposed REnKA is based on Gaussian approximations of the sequence of Bayesian posteriors and thus, in general, asymptotic convergence of posterior expectations cannot be ensured. Nevertheless, our numerical results demonstrated that REnKA is robust with respect to tuneable parameters and provides reasonably accurate estimates of the posterior mean and variance with a computational cost affordable for practical RTM processes.
While measurements have been widely used with ad-hoc approaches to estimate permeability of preform in RTM, to the best of our knowledge, this work constitutes the first investigation that uses a Bayesian inverse modeling framework for the inference of preform permeability under the presence of uncertainty. From the numerical investigations that we conducted some conclusions and recommendations can be made with relevant implications to practical RTM settings. In particular, our synthetic experiments indicated that, when a small number of sensors (5 sensors in 1D and 9 sensors in 2D) are used to measure pressure, observing the front/moving domain can substantially reduce the uncertainty (variance) of the estimates of the log-permeability. This is particularly relevant in real experiments when the number of pressure sensors are usually limited. However, when the inversion is conducted with only measurements of the moving front, the reconstruction of the main spatial features of the true permeability are not recovered accurately.
Our results also display the benefit of the proposed sequential approach in updating our knowledge of the log-permeability as soon as measurements become available. Indeed, inverting measurements of pressure and front frequently in time, enabled us to reduce the uncertainty in the log-permeability. While the reduction of the uncertainty is mainly achieved within the region occupied by the resin at a given time, a decrease in the uncertainty (with respect to the prior) can also be observed in an unfilled region close to the front. Such a reduction of this uncertainty via the Bayesian posteriors can be valuable for decision-making purposes with the aim of an optimal control of RTM in real time. Finally, our numerical investigations show that the observation noise in pressure measurements and front location have a substantial effect on the estimates of log-permeability and its uncertainties. Indeed, it comes as no surprise that more accurate measurements (e.g. 1%) result in estimates of log-permeability concentrated around the truth. Again, giving the limitation of deploying many pressure sensors within a real RTM scenario, it is then essential to use high precision pressure sensors to achieve enough confidence in the posterior uncertainties of the inferred permeability.
Although the context of this work is the resin injection in RTM processes, the Bayesian framework at the core of the proposed methodology, and the corresponding Gaussian approximations emerged from the proposed REnKA are generic, flexible, and thus transferable to a wide class of inverse problems constrained with PDEs with moving boundaries. where t f is a reference time and u 0 is a reference (constant) log-permeability. For simplicity we choose which enable us to transform (2.7)-(2.8) into where, with a slight abuse in the notation, the variables p, u, Υ, x and t are now the dimensionless variables. Similarly, the dimensionless filling time (2.9) is given by Let us note from (A.2) and (2.6) that .
Recall that {t n } N n=1 is the collection of observation times with 0 < t 1 < t 2 < · · · < t N and {x m } M m=1 are the pressure measurement locations (0 < x 1 < · · · < x M ≤ x * ). We now prove two technical lemmas which are needed for the proof of Theorem 2.2.
for all t,t ∈ [t 1 , τ * ,u ]. The constant A u may depend only on ||u||, t 1 and x * .
Proof: From the Mean Value Theorem it is not difficult to see that where M u,v ≡ e max{||u||, ||v||} .

From (A.3) it follows
We use (A.35), (A.11) and the fact that F v (x m ) ≤ F v (x * ) to rewrite (A.34) as follows From similar arguments it is easy to see that We use (A.36)-(A.37) and Lemma A.2 to conclude that 38) which proves the continuity of u → p u (x m , t ∧ τ * ,u ) for all t ≥ t 1 . The continuity of G p n (u) then follows.
Appendix B. SMC sampler and pcn-MCMC algorithm.
In Algorithm B.1 we display the pcn-MCMC method that we use for the mutation step in the SMC sampler of [21] discussed in Section 3.2 and summarized in Algorithm B.2 below.
Note that F u defined in (2.6) can we written as where u s = u(x s ). Therefore, (C.1) can be approximated by where we have replaced H(x) with its smooth approximationĤ(x) = 1 2 + 1 2 tanh rx (with r = 300) [6]. We consider a temporal domain [0, 0.4] discretized with K points t k = k∆t, where ∆t ≡ t f /K. An implicit backward Euler scheme applied to the dimensionless version of (2.4) yields where Υ k ≡ Υ(t k ). For the approximation of (C.3), we use (C.2). The solution of the resulting nonlinear equation is implemented in MATLAB by means of the routine fzero. Once Υ k is computed, we evaluate the pressure field from p(x, t k ) = 2 − F u (x) F u (Υ k ) (C.4) at the mesh points x s+1/2 defined earlier.
Appendix D. Performance of the SMC sampler. In this section we discuss the performance of the SMC sampler (Algorithm B.2) applied in Section 3.4.2. The first measure of performance is ESS defined in (3.7) as a function of iterations given by n p=1 q p , with q p given in (3.19)). This is displayed in Figure D.1 (left). By design of φ n,r , ESS oscillates between the tunable parameter J thresh (when φ n,r < 1) and larger values when φ n,r > 1. Note that more iterations (i.e. tempering) is needed for computing µ 1 . This comes as no surprise due the substantial difference between µ 0 and µ 1 , which arises from the assimilation of data at the first observation time t 1 . By ensuring that the ESS does not fall below the threshold J thresh , we avoid a sharp failure in the SMC sampler as discussed in Section 3.2.
We now monitor performance of the mutation step (Step (e) of Algorithm B.2) of the SMC sampler that controls the diversity added to the population of particles. As suggested in [21], this can be monitored by the average acceptance ratio with respect to the number of particles in Algorithm B.1. We plot this ratio in Figure D.1 (middle). For the present work, we have tuned β in the proposal of Algorithm B.1 so  that the average acceptance ratio for the measures is no less than 30%. Additionally, it is suggested in [21] to monitor the quality of the efficacy of the mutation step by means of where u (j) k,n,r (0) (resp. u (j) k,n,r (M )) denotes the k KL coefficient of the jth particle u (j) n,r before the mutation step (resp. after M MCMC iterations). In Figure D.1 (right) we display the plots of max k J k,n,r , min k J k,n,r and the average with respect to k of J k,n,r . As suggested in [21], this quantity should not be less than 0.09, which is indicated by a horizontal dotted line in Figure D.1 (right). For further discussions of the role of (D.1) in assessing the mutation step in SMC we refer the reader to [21].
Finally, in order to assess stability of the scheme with respect to the prior ensemble, we consider 4 different runs of the SMC sampler with a different set of initial J = 10 5 particles. The errors (3.23) with respect to the truth for all these cases are displayed in Table D.1. These results indicate that there is indeed quite a good agreement in the estimator of SMC provided by the mean of the particles. Right: Plots of max k J k,n,r , min k J k,n,r and the average with respect to k of J k,n,r .