Adaptive regularisation for ensemble Kalman inversion with applications to non-destructive testing and imaging

We propose a new regularisation strategy within the classical ensemble Kalman inversion (EKI) framework for estimating parameters in PDEs. The regularisation strategy consists of: (i) an adaptive choice for the regularisation/inflation parameter in the update formula in EKI, and (ii) criteria for the early stopping of the scheme. Our main contribution is the selection of the regularisation parameter which, in contrast to existing approaches, it does not rely on further parameters which often have severe effects on the efficiency of EKI. The proposed approach is aimed at providing EKI with computational efficiency and robustness for addressing problems where the unknown is a heterogeneous physical property with possibly sharp discontinuities arising from the presence of anomalies/defects. In these settings, for EKI to produce accurate estimates of the truth, the unknown needs to be suitably characterised via a parameterisation that often increases the complexity of the identification problem. We show numerically that the proposed approach can produce efficient, robust and accurate estimates under those challenging conditions which tend to require larger ensembles and more iterations to converge. We test our framework using various parameterisations including one that combines a truncated Whittle-Matern (WM) level-set function with other WM fields to characterise spatial variability of the physical property on each region. We use these parameterisations to solve PDE-constrained identification problems arising in (i) medical imaging where the aim is to detect the existence and properties of diseased tissue and (ii) non-destructive testing (NDT) of materials from data collected during manufacturing processes. We provide comparisons against a standard method and demonstrate that the proposed method is viable choice to address computational efficiency of EKI in practical/operational settings.


Introduction
Ensemble Kalman Inversion (EKI) [1,2,3,4,5] is a collection of algorithms that import ideas from the ensemble Kalman Filter (EnKF) developed in [6]. While EnKF was devised for assimilating data into models for numerical weather prediction and oceanography, EKI has been applied to solve parameter identification problems arising from multiple disciplines. The initial versions of EKI were proposed for the calibration of oil-reservoir models [7,8], and then transferred in [1,2] to more generic PDE-constrained inverse problems settings. The development of EKI as an iterative solver for parameter identification problems has lead to numerous application including the calibration of climate [9], turbulent flow [10], finance [11] and biomedical [12] models. EKI has been used for imaging and non-destructive testing including electrical impedance tomography [13], seismic inversion [14], characterisation of thermophysical properties of walls [15] and composite materials [16]. More recently, EKI has also been applied for the solution of machine learning tasks [17,18].
There is clear promise for the potential use of EKI as a practical and operational tool for PDEconstrained calibration and imaging arising from multiple applications in science and engineering. However, most EKI algorithms still rely on the appropriate selection of user-defined parameters that control the stability, accuracy and and computational efficiency of EKI. Unfortunately, the lack of a general theory for the convergence of EKI means that there is no principled approach to select these parameter in an optimal fashion. The seamless version of EKI proposed in [19] and further developed in [20,21,22] has provided enormous theoretical advances for understanding EKI within the more general frameworks of Stochastic Differential Equations (SDEs). However, to the best of our knowledge, the selection of those crucial EKI parameters are still chosen empirically. Among those parameters, the choice of the regularisation (inflation) parameter in EKI, or alternatively, the step size from the discretisation of the SDE formulation of EKI often depends on additional tuning parameters which can only be informed after problem-specific numerical testing which maybe computationally intensive.
The aim of this work is to introduce a simple, yet computationally efficient, regularisation strategy for EKI that does not rely on further tuning parameters. Our main focus is large-scale/high resolution imaging/identification settings in which there are two fundamental challenges: (i) the forward model is computationally costly and (ii) the unknown must be parameterised in a highly-complex manner so that key features from the truth can be extracted via EKI. Both challenges are intertwined since the latter means that EKI requires large ensembles and more iterations to achieve accurate identifications; thus the total computational cost of EKI can become unfeasible. An efficient and robust regularisation strategy within EKI is a key requirement in these practical settings.

The inverse problem framework with ensemble Kalman inversion.
We work on a generic setting where the properties that we wish to identify are functions which, in turn, play the role of inputs (e.g. coefficients) of Partial Different Equations (PDEs) describing the underlying experiment/process. The observable variables in these experiment/processes are functionals of the solution of these PDEs. Under this inverse problem setting, the forward problem can be written in terms of a nonlinear operator F : X → Y that maps physical properties from an admissible space X to the space of observable quantities Y. Our work is focused on using EKI for solving the classical (deterministic) inverse problem that consists of estimating the underling physical property of the medium, denoted by κ † , given noisy measurements of F(κ † ), which we assume are given by where η is the unknown measurement error. We further assume η is drawn from a centred Gaussian measure with known covariance Γ. In order to address this inverse problem, we define a class of suitable parameterisations P : U → X that enable us to characterise our estimate of the unknown as κ = P(u), i.e. in terms of an input (function) u ∈ U which we calibrate within the EKI framework. We pose the parameterised (in terms of u) inversion via the following weighted least-squares problem u * = arg min where G = F • P is the forward (parameter-to-output) map, and S 0 ⊂ U is a user-defined subspace of admissible solutions. We use EKI to approximate (2) and to provide an estimate of κ † via κ * = P(u * ).
For the PDE-constrained identification problems that we wish to address, G is often a compact operator which leads to the ill-posedness of (2) in the sense of stability [23,24]. Although numerous regularisation techniques [25,26] can be used to address ill-posedness, most of them require the computation of the Frechét derivative of the forward map, as well as the corresponding adjoint operator. This constitutes a substantial limitation in many practical applications where the map F is only accessible via commercial software simulations with no adjoint built-in functionalities. Such a practical limitation has given rise to a large body of work on EKI techniques that steam from EnKF [6] and which, in turn, does not require derivatives of the forward map.
While there are several versions of EKI algorithms to approximate (2), here we focus on the classical, perturbed-observation EKI displayed in Algorithm 1. This generic version of EKI involves selecting an initial ensemble of J particles {u (j) 0 } J j=1 ⊂ U. Then, each particle is iteratively updated according to the following formula u (j) n+1 = u (j) n + C uG n (C GG n + α n Γ) −1 (y + √ α n ξ n − G(u (j) n )) where α n is a tuning (regularisation) parameter, ξ n ∼ N (0, Γ) is perturbation of the data, and C uG n and C GG n are empirical covariances defined in (7)- (8). The running estimate of the unknown is obtained via the ensemble mean Using informal arguments it can be shown (see Appendix A and the work in [2]) that the ensemble u n+1 can be seen as an ensemble approximation of m n+1 = arg min where DG denotes the Frechét derivative of G, C n is a covariance operator that we define in Appendix A, and S 0 ≡ span{u (j) 0 } J j=1 is the subspace generated by the initial ensemble. If C n is the identity operator, (5) is the standard Levenberg-Marquardt (LM) scheme [27] applied for the computation of (2). Note that (5) can also be interpreted as iterative Tikhonov regularisation applied to the linearisation of G. The link between EKI and (5) is very useful because (i) it motivates EKI as a derivative-free solver for (2) and (ii) it reveals the role of α n in (3) as a Tikhonov regularisation parameter. According to the theory in [28], α n must be carefully selected, together with the stopping criteria, in order to ensure the stability of the LM scheme.
Existing approaches for selecting α n in EKI [2,16], some of them which are adapted from the LM work of [28]), relies on tuning parameters that, unless chosen very carefully, can lead to unnecessary large number of iterations n * . Since the main computational cost of Algorithm 1 is n * × J, it is clear that a large n * is detrimental to the computational efficiency of EKI.

Our contribution.
In this work we propose an adaptive choice of α n that does not require any tuning parameters. This new selection of α n , to which we refer as the data misfit controller (DMC), only depends on the particles data misfit We motivate our data misfit controller from the Bayesian perspective of EKI in which α −1 n is the difference between consecutive tempering/annealing parameters. The stopping criteria or, alternatively, the number of iterations in EKI, n * , is determined from the constraint n * n=1 α −1 n = 1 which arises naturally from the tempering/annealing setting. We show that this new choice of α n in EKI leads to a stable, robust and computationally efficient algorithm (EKI-DMC). We investigate the performance of EKI-DMC for (i) Electrical Impedance tomography (EIT) where we use EKI to identify electrical conductivity in the Complete Electrode Model (CEM) and (ii) Non-Destructive Testing (NDT) of composites where we estimate defective regions of high/low permeability from pressure data collected during the resin injection phase of a Resin Transfer Moulding (RTM) process. We demonstrate the computational and practical advantages of EKI-DMC over existing approaches in which α n is borrowed from the theory for the LM scheme. In contrast to most existing EKI work, here we test EKI with parameterisations that enable us to identify realistic, possibly discontinuous physical properties. In particular, we investigate the performance EKI-DMC to infer regions characterised via a level-set function which is, in turn, parameterised with anisotropic Whittle-Matern (WM) fields. We further consider the case of possibly heterogenous lengthscales (also parameterised via WM fields) which allows us to infer geometric features of complex geometries which are relevant for practical imaging and NDT settings.
In Section 2 we review the literature on existing choices of α n for EKI. We introduce our new EKI algorithm in Section 3. In particular, in subsection 3.1 we discuss the data misfit controller for selecting α n together with the stopping criteria. The test inverse problems, EIT and NDT, are discussed in subsection 3.2. The parameterisations of the unknown that we infer via EKI are introduced in subsection 3.3. Numerical examples for EIT and NDT are presented in Section 4 and Section 5, respectively. In Section 6 we discuss final remarks and conclusions. For completeness in Appendix A we motivate the classical EKI from the Bayesian tempering scheme, and from which we also draw links between the LM algorithm and EKI.

Literature review
We discuss some existing regularisation strategies for EKI within the inversion setting posed in terms of the unregularised least-squares formulation of (2), and which leads to the classical EKI formulation in (3). We highlight that there is a new alternative EKI methodology and algorithms proposed in [4,5] that arise from regularising (2) with a Tikhonov-like term. In addition, the review is focused on PDEconstrained inversion; for a review of modern Kalman methods in the context of the data assimilation framework we refer the reader to the recent work on [29] and references therein.
2) Measurements y and covariance of measurement errors Γ.
(2) Compute regularisation parameter α n and check for convergence criteria if converged then set θ = 1 and n * = n (3) Update each ensemble member n ← n + 1 end output: {u (j) n * } J j=1 converged ensemble 2.1. EKI as an iterative solver for identification problems.
The initial versions of EKI [1] for generic identification problems proposed to use the classical EnKF [6] update formula 1 as an iterative solver for (2) by introducing an artificial dynamical system. For various PDE-constrained identification problems, the work of [1] numerically showed that this early version of EKI approximated well the solutions of (2) (with S 0 = span{u (j) 0 } J j=1 ) within the first few iterations. However, they noted the algorithm became unstable if it was allowed to iterate after the data misfit (2) had reached the noise level δ defined by 1 the classical EnKF consists of eq (3) with α n = 1 fixed throughout the iterations where u † is the truth. This lack of stability led to the work of [2] where links between EKI and the regularising Levenberg-Marquardt (LM) scheme of [28] were first established and used to develop a regularising version of EKI. For the LM scheme, the work of [28] ensures that, under certain assumptions of the forward map G, the scheme in (5) converges to the solution of (2) (as δ → 0) provided that (i) the regularisation parameter α n satisfies where ρ < 1 is a tuning parameter, and (ii) that the algorithm is terminated at an iteration level n * determined by the following discrepancy principle where τ is another tuning parameter that must satisfy τ > 1/ρ. In [2] (see also Appendix A), these regularisation strategies from LM were adapted for the selection of α n in EKI via using derivative-free Gaussian approximations in (10)- (11). We refer to the approach from [2] as EKI-LM (see Algorithm 2).
The numerical results of [2] showed that EKI-LM enabled stability and accuracy for sufficiently large ensembles. Further work that has explored EKI-LM can be found in [30,3,31] as well as some practical applications including seismic tomography [14], modeling of intracraneal pressure [12] and time fractional diffusion inverse problems [32].
Despite of addressing stability in EKI, the approach EKI-LM suffers from a potential practical limitation that arises from the fact that it relies on the tuning parameters ρ and τ in (12). Larger ρ's yield larger α n and in turn more iterations to converge. Smaller ρ's, while desirable for computational efficiency, result in larger τ 's which can lead to larger data misfit and possible loss of accuracy from stopping too early (via (13)). Selecting these parameters in a computationally optimal fashion becomes even more urgent when EKI is combined with complex parameterisations of the unknown. As we discuss in subsection 3.3, there are cases for which, instead of using EKI directly on the physical property that we wish to infer, we need to parameterise the unknown to be able to capture properties that are not necessarily encoded in the initial ensemble. For example, in [2] the LM approach for EKI was applied with a level-set parameterisation of the unknown in order to characterise discontinuous properties (i.e. discontinuous conductivity in the context of EIT). In comparison to the simpler case in which EKI directly estimates a physical property of interest, [2] found that not only a larger ensemble size was needed but also more number of EKI iterations to achieve converge. The application of EKI-LM with level-set parameterisations for seismic imaging in [14] reported up to 40 iterations to achieve convergence. The numerical results of [3] also show that when EKI is combined with various others parameterisations of the unknown, the number of iterations of EKI can become large even for simple 1D and 2D forward modelling settings. We aim at addressing these very same issues with the selection of α n that we propose in Section 3.1.

EKI as Gaussian approximation in linearised Bayesian tempering.
Although the goal of most of the existing applications of EKI is to solve the deterministic problem in (2), the role of EKI within the Bayesian setting for parameter identification can be useful for identifying suitable choices of the regularisation parameter α n . In the Bayesian setting we put a prior measure, µ 0 (u) = P(u), on the unknown u that we wish to infer. Given measurements, y, the Bayesian inverse Algorithm 2: EKI with LM selection of (EKI-LM) Input: Same from Algorithm 1, ρ < 1, and τ > 1 ρ as well as the noise level δ.
(2) Compute α n such that if then set n * = n.

break
(3) Update each ensemble member using (3) (see also Step 3 from Algorithm 1) problem consists of approximating the posterior µ(u) ≡ P(u|y) which, from Bayes' rule [33] is given by where we have made the standard assumption that y = G(u) + η with η ∼ N (0, Γ). Modern computational approaches [34,35,36] to tackle high-dimensional Bayesian inverse problems, use the tempering approach that consists of introducing N intermediate measures {µ n } N n=1 between the prior and the posterior. These measures are defined by where {φ n } N n=1 are tempering parameters that satisfy: Note that n = 0 and n = N + 1 in (15) yields the prior, µ 0 , and posterior (µ N +1 = µ), respectively. From expression (15) we obtain the following recursive formula for the intermediate measures: where Here we use the same notation that we use for the regularisation parameter in EKI (see eq. (3)), because the ensemble computed at the nth iteration of EKI, is nothing but an approximation of a Gaussian measure ν n which, in turn, approximates the intermediate tempered distribution µ n in (15) (see [16] and Appendix B). Even when J is large, these approximations are only exact in the linear-Gaussian case (i.e. G linear and µ 0 Gaussian). Therefore, the final (converged) ensemble from EKI does not, in general, sample the posterior unless further modifications to the algorithm are made [37]. Nevertheless, under the tempering interpretation, controlling the regularisation parameter α n in EKI means to gradually transition between prior and posterior in order to facilitate more accurate sampling of the intermediate measures. For further details on tempering-based (fully) Bayesian algorithms we refer the reader to [34,35,36,38,39].
Although the Bayesian perspective of EKI as a Gaussian approximation within the tempering setting was initially mentioned in [19] and further developed in [16], the early work of [7,40] established a strong link between EKI and the Bayesian setting which led exactly to the same EKI scheme from Algorithm 1. However, instead of using tempering, they used a heuristic approach in which the data was inverted multiple times with noise inflated by √ α n . From algebraic manipulations they figured out that, in order to accurately sample from the posterior in the linear-Gaussian case, their inflation parameter α n must satisfy where n * is the total number of EKI iterations. Note that, in the tempering setting, (19), with N = n * intermediate measures, is simply a consequence of the definition of α n in (18) as well as the definitions of φ 0 and φ N +1 in (16). Moreover, in the general Bayesian tempering settings, expression (19) holds for the general nonlinear case and without any assumptions on the prior. Nonetheless, those considerations for the linear-Gaussian case from [40] enabled them to (i) justify the use of (19) in the nonlinear case and (ii) to propose a simple selection of α n given by α n = n * with n * selected a priori (which trivially satisfies (19)). This approach, referred to as Ensemble Smoother with Multiple Data Assimilation (ES-MDA), has been popularised in petroleum engineering applications (see [41] and references therein). However, from the insight gained from the link between EKI and the LM scheme, this selection of α n was discouraged in [42] since the stability of EKI, as shown in [42,2], requires α n to be large at the beginning of the iterations, and gradually decrease as the data misfit approaches the noise level. Further versions of ES-MDA [43] adopted selections of α n similar to those initially proposed in [2] based on the LM scheme. The link between EKI and the Bayesian tempering setting has been further explored in [16], where the selection of α n is borrowed from the adaptive-tempering Sequential Monte Carlo (SMC) method of [35]. With this approach α n is selected so that where J * is a tuning user-defined parameter and The left hand side of (20) is the effective sample size that, in SMC methods, is used to determine the quality in the population of particles that approximate each intermediate measure µ n . For further details we refer the reader to [16], where the selection of α n according to (20) was implemented in a batch-sequential EKI framework to sequentially solve a similar inverse problem to the one we describe in subsection 3.2.2. The same EKI methodology was applied in [15] for parameter identification of the heat equation, including identification of thermal conductivity and heat capacitance given boundary measurement of heat flux. Both time-dependent applications tackled in [15,16] involved the inversion of small number of measurements (e.g. < 30) at each observation time, and with relatively large noise informed by measurement protocols. However, the selection of α n via solving (20) becomes problematic when the number of observations is large and/or a the measurement noise is small. Indeed, note that (21) involves the computation of the likelihood between consecutive measures (17) which, in turn, can take very small values unless large α n 's are used for solving (20). This means that more iterations may be needed to satisfy condition (19). For the problems that we formulate in subsection 3.2, the number of measurements can be higher (> 100)) to those that can be addressed with these approaches. Furthermore, we note that (20) requires a tuning parameter which may substantially affect the efficiency or the stability depending of whether J * is large or small, respectively.

EKI as a discretisation of Stochastic Differential Equations
The pioneering work of [19] has shown that EKI can be derived as a discretisation of an SDE system for the ensemble of particles in EKI. This so-called continuos-time limit or seamless formulation of EKI has led to further theoretical understanding including its gradient-flow structure. In addition, using alternative discretisation schemes of the seamless formulation of EKI, leads to new EKI algorithms in both the context of optimisation [20,21,17] and sampling within the Bayesian approach [22].
In the seamless formulation of EKI, the regularisation parameter α n from the classical EKI becomes the inverse of the mesh-size/discretisation step. The work of [17] proposes to choose this parameter according to where α 0 and are user defined parameters, · F denotes the Frobenius norm and U is a matrix with entries U j,k = (G(u . To the best of our knowledge, this selection of α n has only been used for new different EKI algorithms including the ones arising from Forward-Euler [17] and the implicit split-step method [22]. While (22) is also a perfectly valid choice of α n for classical EKI, we emphasise that it relies on the choice of the tuning parameters α 0 and .

The proposed EKI framework
The aim of this section is to introduce a new approach to select α n in the classical EKI setting given in (3). In subsection 3.2 we introduce the identification problems that we use to test the proposed EKI algorithm, using the parameterisations of the unknown introduced in subsection 3.3.

Data misfit controller: an adaptive regularisation strategy
The proposed approach for selecting α n is motivated by the tempering setting discussed earlier and for which the expression in (19) must be satisfied. We assume that the dimension of the observation space is finite: Let us first make the following observation concerning the tempering scheme introduced in the previous section.
Remark 3.1. We can view (17) as an iterative application of Bayes rule, where at the nth iteration level we have a prior µ n (u) and a likelihood defined by the observational model In other words, (17) defines a sequence of Bayesian inverse problems similar to the original one 2 but with a Gaussian error √ α n η that has a covariance matrix Γ inflated by α n . We can then think of µ n+1 given by (17) as the distribution of u|y under the observational model (24) and prior µ n .
Let {u (j) n } J j=1 be EKI samples that approximate the intermediate measure µ n . The proposed approach is based on the assumption that it is possible to find an α n such that the ensemble {u is consistent with the observational model (24). In other words, there is α n such that, for all j ∈ {1, . . . , J}, In (25) y is the actual instance of the observed data, and which, as before, we assume is given by where u † is the truth and η is a particular realisation of N (0, Γ). If there is, indeed, an α n such that (25) holds, then where we have used the fact that each Γ −1/2 η (j) 2 is a chi-square random variable with M degrees of freedom. Furthermore, the expectation of the random variable in the right hand side of (28) is M . Hence, we propose to estimate α n from (28) using M in the right hand side of (28). In other words, which entirely depends on the squared data misfit averaged over the ensemble. The regularisation of EKI is thus controlled by the ensemble data misfit. In order to make sure that (19) is satisfied, we propose Let use denote by n * ∈ N ∪ {0}, the iteration such that α −1 n * = 1 − θ n * −1 for the first time, and for which the algorithm is terminated. Indeed, from (31) we note that and so (19) is, indeed, satisfied. Thus, our selection of α n yields a well defined collection of tempering parameters (from (18)). We use EKI-DMC to refer to the EKI algorithm that uses this data-misfit controller for selecting α n ; the algorithm is summarised in Algorithm 3. Note that, in contrast to EKI-LM, we still update the ensemble to obtain {u (j) n * +1 } J j=1 at the final iteration of EKI-DMC. Further theoretical aspects of the DMC within a seamless formulation of tempering are studied in [44,45].
with η (j) ∼ N (0, Γ). Equations (26) and (32) imply that, for this case, particles from the initial ensemble are within the same noise level from the truth. For most cases, however, the data misfit for the initial particles is several orders of magnitude greater than M , and thus, n * > 0. From the definition of n * we see that, for all j < n * , α j is given by (29) and so α j > 0. Furthermore, for n < n * , it follows from (31) that and so Therefore, before termination, the algorithm honours a discrepancy principle for the root mean squared error of the data misfit, provided that we estimate the noise level in (9) via δ ≈ √ M [24] (i.e. using (26) and the fact that Γ −1/2 η 2 is a chi-squared with M degrees of freedom).

Test problems
In this subsection we introduce the two identification problems that we use to showcase the capabilities of EKI-DMC. The aim in both problems is to identify a heterogenous (possibly discontinuous) physical property κ † of a medium with physical domain D ⊂ R 2 . For the first test problem the identification is based on boundary measurements from a stationary process. The second test problem is time depend and measurements are observed in the interior of D at several observation times.

Test Problem 1: Electrical Impedance Tomography.
As a prototypical imaging problem we use Electrical Impedance Tomography (EIT) with the Complete Electrode Model (CEM). Given electric currents {I k } me k=1 injected through a set of surface electrodes {e k } me k=1 placed on the boundary, ∂D, of D, the CEM consist of finding where v is the voltage in D and {V k } me k=1 are the voltages on the electrodes. The dependent variables [v, {V k } me k=1 ] are given by the solution to [46]: where κ is the electric conductivity of D and {z k } me k=1 are the electrodes' contact impedances. We consider an experimental setting consisting of n p current patterns For each of these patters {I j,k } me k=1 , we denote by {V j,k } me k=1 the prediction of voltages at the electrodes defined by the CEM (34)- (37). For simplicity we assume that the contact impedances of the electrodes are known. We define the map that for every conductivity, produces voltage measurements. EIT consist of finding the conductivity κ † of D given measurements of V † = F(κ † ). For a review of the EIT problem we refer the reader to [47]. Recent uses of EIT for medical imaging can be found in [48,49] 3.2.2. Test Problem 2: Estimating permeability of a composite preform during resin injection.
The second problem is in the context of Resin Transfer Moulding (RTM) process, which is commonly applied to manufacture composite materials [50]. Composites are high-performance materials widely used in automotive and aerospace industries [51,52]. The inverse problem consists of estimating the permeability of a composite preform from measurement collected during the injection of resin within the RTM process; see [16] and references therein for a review of existing approaches to solve this inverse problem. Preforms are designed to have homogeneous permeability so that resin can fill the preform uniformly. However, variations arising from preform fabrication can lead to regions of anomalous permeability that are detrimental to the material properties. A particular concern is the presence of possible channels of very high permeability that can substantially change flow patterns. The appearance of these disturbances, also known as race-tracking, can cause air entrapment and, in turn, material defects [53].
We model the preform as a porous medium with physical domain D, 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, initially filled with air at a pressure p 0 , is infused with resin with viscosity µ through an inlet boundary ∂D I at a pressure p I . The resin moves through D occupying a time-dependent domain D * (t) ⊂ D, which is enclosed by the moving boundary Υ(t) and the appropriate parts of ∂D. A diagram of the setting for this problem in 2D is illustrated in [16, Figure 1].
From conservation of mass and Darcy's law, it follows that resin's pressure p(t, x) satisfies [50,54,55]: with the following initial and boundary conditions In (62)-(65) V (x, t) is the velocity of the point x on the moving boundary Υ(t) in the normal direction at x, n(x) and n(x, t) are the unit outer normals to the corresponding boundaries. In the considered one and two dimensional cases of this problem, we can view the velocity of the moving boundary as the following derivative [54]: Note that p 0 and p I are control variables and hence known. For simplicity we assume ϕ, µ are known and so we focus the inverse problem on κ. The dependent variables from (38)- (45) are p(x, t) and Υ(t). We could potentially observe these two variables (see [16] for further details), but for simplicity, here we only consider pressure measurements from sensors placed at locations {x k } K k=1 observed at times {t n } Nt n=1 . From (38)- (45) we define the map and we pose the inverse problem of finding κ † given measurements of F(κ † ) = P † .

Parameterisations of EKI
In this subsection we introduce three maps P 1 , P 2 and P 3 that we use to parameterise the unknown physical property κ † that we wish to estimate using the EKI framework introduced earlier.
The objectives of these maps are to (i) enforce possible known features of the truth κ † (e.g. smoothness/discontinuities) and (ii) control/adjust the spatial variability in the values of κ(x) within its domain of definition. Selecting such a parameterisation is essential for the effectiveness of EKI for the practical applications that we discuss in the previous subsection. In most academical examples, the common choice is simply P(u) = κ, or more precisely P(u) = log κ (to ensure that κ is positive as required for the test models above). With this naive parameterisation, EKI is applied directly on the space of physical properties. Unfortunately, this poses severe restrictions in practical settings because of the invariant subspace property discussed in Remark Appendix B.3. Indeed, if EKI is applied directly on κ (or log κ), the regularity and geometric features of the truth, κ † , must be encoded in the initial ensemble so that EKI can produce an accurate identification of the truth 3 . Of course, when both truth and particles from the prior ensemble are all samples from the same distribution, reconstructions tend to be very optimistic albeit this is clearly a statistical inverse crime. When the truth is not selected in this fashion, the estimates computed via EKI, as shown numerically in [2] and more recently in [3], can be substantially inaccurate. This limitation has been recently addressed in [3] by means of adequate parameterisations of the underlying field which includes various parameters that can be used to characterise the unknown spatial variability of the physical property that we wish to identify via EKI. In this work we follow the parameterisation from [3] with additional components from the work of [56] that allows us to incorporate even more features (e.g. anisotropy, and variable lengthscales) of the spacial variability of the unknown physical property.

Parameterisation P 1 . Smooth properties
Let us first discuss the case in which prior knowledge suggest that the unknown physical property κ † is a smooth function. For simplicity we only consider the 2D case that we use in the numerical experiments of Sections 4-5 but we emphasise that the approach can be used for 1D and 3D settings.
Let us introduce what we call Whittle-Matern (WM) parameterisation of the unknown κ(x) that we wish to infer via EKI. The WM parameterisation involves a positive smoothness parameter denoted by ν, an amplitude scale σ, and two intrinsic lengthscales, L 1 and L 2 , along the horizontal and vertical direction, respectively. Given Θ = (ν, σ, L 1 , L 2 ) ∈ R + × R × R 2 + , we define an operator W Θ that maps with robin boundary conditions on ∂D: where ζ R is a tuning parameter. In (48) Γ denotes the gamma function, I is the identity operator, and The proposed WM parameterisation of κ is defined by where λ a positive scaling factor for κ. Note that (50) enforces that the physical property is positive, as required for test models presented in subsections 3.2.1-3.2.2. We can succinctly write (50) as, The motivation behind this parameterisation comes from the theory of Gaussian random fields (GRFs) [57,58]. In particular, the work of [58] shows that if ω ∈ H −1+ is Gaussian white noise (i.e. ω ∼ N (0, I)) then where C Θ is the covariance operator induced by the Matern autocorrelation function defined by where K ν is the modified Bessel function of the second kind of order ν, and It can also be shown that if log κ ∼ N (log λ, C Θ ), then almost surely log κ ∈ H ν− (D) [59] which further shows the role of the smoothness parameter ν.
Our choice for the boundary conditions (BCs) in (49) follows from the work of [56] that shows that Robin BCs are better suited, compared to Neumman and Dirichlet, to alleviate (via the appropriate choice of ζ R ) undesirable boundary effects which arise from the discretisation of GRFs. Let us reiterate that our goal here is to introduce parameterisations for EKI and then suitable initial ensembles on the parameters. Whether the initial ensemble yields Gaussian (or log-Gaussian) properties is not our main concern. However, it would not be advisable to select parameterisations that, for example, restrict the values of the physical property near the boundary which is something that we would expect if ζ R = 0 in (49).
The WM parameterisation in (50) allows us to incorporate the smoothness and lengthscales of the underlying field as part of the unknown that we can estimate with EKI. In contrast to the work of [3] where the lengthscales are estimated under isotropic assumptions, here we consider the anisotropic case which, as we show in the numerical experiments of Section 5, is essential in practical settings where the unknown have rapid changes along one particular direction. We focus on vertical/horizontal anisotropy but a rotation matrix can be further introduced in (48) to characterise properties with an arbitrary (and unknown) direction of anisotropy [56] A straightforward approach to generate the initial ensemble for EKI with the parameterisation from (50) is to specify (hyper prior) densities, π λ , π L 1 , π L 2 , π σ and π ν , and produce samples: Remark 3.3 (KL expansion). For some geometries of D, the eigenfunctions and eigenvalues of a Matern covariance C Θ are available in closed form [59]. In this case, an equivalent formulation of the WM parameterisations can be defined in terms of the spectral decomposition of C Θ (i.e. KL expansion under the prior). Obviously, for computational purposes WM parameterisations can be defined on a larger (simple) domain that encloses D and simply restrict the physical property to the domain of interest. While the spectral/KL approach is more standard, the approach we use here based on the work of [57,58] allows to naturally extend GRF that have non-constant lengthscales in (48). As we show in subsection sections, this is key to ensure our parameterisation capture key geometric features of the unknown. For the experiments that we discussed in the subsequent sections, the spatial variability in the unknown can be captured quite effectively only via ω(x). However, we recognise that including a variable σ can be beneficial in some cases as reported in [56]. Similarly, in the context of the experiments reported later, where the aim is mainly to recover medium anomalies, sensible (fixed) choices of ν are sufficient provided that the lengthscales is properly estimated via EKI. Nevertheless, as shown in [3], EKI can be used to estimate the smoothness parameter ν. Note that, in order to keep a parameter (i.e. a component of u) fixed within EKI, it suffices to produce the appropriate component equal to the desired value for all members of the initial ensemble; this comes from the fact that C u,G (see eq. (8)) vanishes for all those components.

Parameterisation P 2 . Piecewise-smooth functions (constant lengthscales).
In order to characterise piecewise-smooth functions we use the level-set approach initially proposed in [60] for deterministic inverse problems and, more recently, for fully Bayesian [61,59] and EKI [2] settings. Our main modelling assumptions is that the unknown property has a background value potentially heterogeneous, and that possible anomalies/defects consist of regions with (also possibly heterogeneous) higher/lower values than those in the background field. More specifically, let us first define the parameterisation where κ b , κ l and κ h denote the background, low-value and high-value fields that characterise the unknown physical property. The level-set function denoted by f (x) determines the background Ω b ≡ {x : ζ 1 < f (x) ≤ ζ 2 } as well as the region of low Ω l ≡ {x : f (x) ≤ ζ 1 } and high Ω h ≡ {x ∈ f (x) > ζ 2 } values. We assume the unknown functions {κ ι } ι∈{l,b,h} and f are smooth fields with some variability that we enforce via a second level of parameterisation in terms of P 1 introduced earlier. More specifically, we consider where P 1 is defined according to (50) and (53) with (54)- (55) we can write where The selection of the initial ensemble for each u ι can be done similarly to the one in (52).

Parameterisation P 3 . Piecewise-smooth functions (spatially variable lengthscales).
In some cases the accuracy of the estimation of anomalous/defective regions can be improved by using spatially variable lengthscales L 1,f (x) and L 2,f (x) in the level-set function parameterisation of (55). For example, if the physical property displays anomalies/defects with geometric features of different sizes and shapes, the constant lengthscale framework will not produce accurate identifications. To address this limitation, we assume that L 1,f and L 2,f in (55) are heterogeneous, and parameterised in terms of P 1 as before: If we substitute (58) in (57) we obtain where all inputs are comprised in The level-set parameterisations introduced here have some limitations including the number of unknown regions, and the fact that, by construction, the high-value and low-value regions cannot intersect. These limitations arise from our modelling assumptions which we made only for simplicity and clarity in the subsequent numerical investigation. Other settings can also be applied given the black-box flexibility of EKI. This include for example, dealing with various intersecting regions via methods that use multiple level-sets [62], as well as other characterisations of the unknown such as the pluri-Gaussian method of [63] 3.4. Implementation For the numerical experiments discussed in the following sections, Algorithms 2 -3 are implemented in MATLAB. Let us recall that for these algorithms we need to construct the forward map G = F • P where F is the operator induced by the forward model (e.g. those defined in subsections 3.2.1-3.2.2) and P is any of the parameterisations defined earlier. The numerical implementations of F for each of these test problems are discussed in subsequent sections. However, we reiterate that EKI uses F in a black box fashion. Thus, we only need pass input κ(x) = P(u) that we compute from the appropriate parameterisation, and retrieve the corresponding vector of model predictions G(u) = F(κ).
For the evaluation of the WM parameterisation κ = P 1 (u), we solve (48)-(49) using the techniques from [58] and which restrict us to the cases in which ν ∈ N. The discretisation of the operator I − ∇ · diag(L 1 , L 2 )∇ with BCs from (49) is performed via cell-centred finite differences. The PDE in (48)-(49) is solve in a square domain, equal to, or enclosing D (the domain of definition for the PDE encoded in F). The actual field κ(x) that we pass into F is an interpolation of P(u) on D. The parameterisations P 2 and P 3 are based on the truncation of the level-set function, f (x) so the implementation is straightforward once all the fields κ ι (ι ∈ l, b, h) and f (x) are computed. Given these construction of G, the rest of the steps in Algorithms 2 -3 are computed in a straightforward manner.

Measures of performance for EKI
Given the ensemble {u (j) n } computed via EKI at the n iteration of the scheme (n = 0 corresponds to the prior ensemble), our estimate of the unknown property κ † is given by We measure the accuracy in terms of the relative error with respect to (w.r.t) the truth defined by We often visualise some transformed ensemble members (mainly for the initial ensemble n = 0), i.e.
but note that our estimate κ n in (61) does not involve taking the average of the particles in (63); this would be particularly detrimental for P 2 and P 3 since averaging (63) will not preserve discontinuities. For these two parameterisations we also visualise an estimate of the level-set function given by We additionally monitor the following data-misfit quantities

Numerical Experiments for Electrical Impedance Tomography
The numerical implementation of the CEM from subsection 3.2.1 is conducted using MATLAB software EIDORS [64]. The experimental setting consists of (i) a circular domain of unit radius centred at the origin, (ii) 16 surface electrodes with contact impedances of values 0.01 Ohms, (iii) an adjacent injection pattern with an electric current of 0.1 Amps, and (iv) measurements at each electrode. All these parameters are assumed known and fixed for the inversions of this section. Synthetic data are generated using the mesh from Figure 1 (right) with 9216 elements, while a coarser mesh of 7744 elements is used for inversions. The total number of measurements is M = 16 2 = 256.

Experiment Exp EIT 1 . Continuous Conductivity.
For the first series of experiments the true conductivity, κ † , is a C ∞ (D) function that we specify analytically; the plot of κ † is displayed in Figure 1 (left). Synthetic data are constructed via y = V † + η where V † = F(κ † ) is computed with the CEM and η is a realisation from N (0, Γ). We chose Γ = diag(γ 1 , . . . , γ M ) where The first term in the right hand side of (68) corresponds to adding 1% Gaussian noise. The second term is added simply to avoid small variances from very small voltage (noise-free) measurements. For this experiment we use the parameterisation of smooth functions, κ = P 1 (u), from (50), with u = (λ, ν, σ, L 1 , L 2 , ω). The unknown consist of 5 scalars and 1 function, ω(x), that we discretise on a 100 × 100 grid. Upon discretisation, the dimension of the unknown is dim(U) = 10, 005. We follow the discussion of subsection 3.3.1 (see eq. (52)) for the selection of the initial ensemble. More specifically, we select J particles u , ω (j) ) ∼ µ 0 . We define µ 0 as follows: where U [a, b] denotes the uniform distribution on the interval [a, b], and δ(g − 3) denotes the Driact measure centred on g = 3.
In Figure 2 (top) we show plots for the logarithm of κ (j) 0 = P 1 (u (j) 0 ) for five of members of the initial ensemble. Note that our choice of smoothness parameter ν = 3 produces an initial ensemble that is quite smooth (recall from 3.3.1 that the draws belong to H 3− (D)). We also observe substantial differences in the degree of anisotropy which arises from our selection of a reasonably wide distribution of intrinsic lengthscales that we use to produce the initial ensemble.

Results from the inversion using EKI-DMC with various choices of J.
Synthetic data and initial ensembles produced as described earlier are used as inputs for EKI-DMC (Algorithm 3) with different choices of ensemble size J: 100, 200, 400, 800. For each choice of J, we conduct 30 experiments with different random selections of the initial ensemble. The plots of (log) κ n * = P 1 (u n * ) (i.e. upon convergence) from one of these experiments, computed for each J, are displayed in Figure 3 together with the truth (right panel). In Figure 4

Further results from one run of EKI-DMC with J = 200.
In Figure 2 (bottom) we displayed 5 members of the final (converged) ensemble of (log) κ (j) n * corresponding to the initial ensemble from Figure 2 (top). Plots of log κ n , at some of the intermediate iterations 1 ≤ n < n * = 10 can be found in Figure 5. In Figure 6 (right) we plot, as a function of n, the values of the means λ, L 1 , and L 2 , of the ensembles 2 } J j=1 , respectively 4 . Note that EKI produces a larger lengthscale in the vertical direction. This comes as no surprise since the truth κ † (see Figure 1 (left)) has two inclusions of lower conductivity with larger correlation along the vertical direction. Finally, in Figure 6 (right and middle) we display, for each of the 30 runs, the relative error w.r.t the truth as well as the (log) data misfit (65), as a function of the iteration number n. We note that EKI-DMC is very robust and accurate across ensembles.

Comparison between EKI-DMC and EKI-LM
We compare the performance of Algorithms 2-3 using the same set of 30 initial ensembles for each algorithm (with J = 200). We consider different choices of the input ρ in EKI-LM and, for simplicity, we set τ = 1/ρ + 10 −6 . In Figure 7 we show boxplots of the error w.r.t the truth (right) and the data misfit DM 1,n * (left) obtained with several choices of ρ for EKI-LM; the results from EKI-DMC for J = 200 are also included in these plots. Similar behaviour is observed for DM 2,n * and DM 3,n * and so these plots are omitted. The number of iterations for EKI-LM to achieve convergence, n * , is displayed in Table 1. For one of the 30 experiments described earlier, in Figure 8 we show the plots of log κ n * computed with both algorithms using the selections of inputs described above. In Figure 9 we display the behaviour of log α n as a function of n, computed from one run (same initial ensemble) of EKI-DMC and EKI-LM (with ρ = 0.6). For this particular run, both EKI-LM and EKI-DMC converged in 10 iterations so the cost of running both algorithms is the same.
The above results suggest that EKI-DMC is more accurate than EKI-LM in terms of the error with respect to the truth and data misfit. While EKI-DMC does not depend on any tuning parameter, the performance of EKI-LM is highly dependent on the selection of ρ (which is specified a priori). From Table 1 we see that, for ρ = 0.6, the computational cost of EKI-LM is similar to the cost of EKI-DMC; i.e. convergence in approximately 10 iterations (in average). For this ρ, the average error computed with EKI-LM (≈ 38%) is, however, larger than the error obtained via EKI-DMC (≈ 0.305%). The improvement in accuracy of EKI-DMC over EKI-LM may not be overly impressive 5 and, indeed, we could always find problem-specific tuning parameters ρ and τ that will approximately yield the same level of accuracy of EKI-DMC. Nevertheless, in practical (real) settings, we may not be able to afford such a thorough and computationally intensive investigation of tuning parameters. For these experiments, if instead of ρ = 0.6 we choose, say ρ = 0.8 (see Table 1), the computational cost of EKI-LM doubles without improving its accuracy with respect to EKI-DMC.

Experiment Exp EIT 2 . Discontinuous (piecewise constant) conductivity.
The true conductivity for the experiments of this section is the piecewise constant function with plot displayed in Figure 1 (middle). The background, low and high conductivity regions have constant values: respectively. Synthetic data are generated in the same way as described in the previous subsection, and with noise covariance from (68). In Figure 10 we show the plots of (log) κ n * computed from one run with EKI-DMC using the same parameterisation (for continuous functions) and initial ensemble from Exp EIT 1 . While these results show that we can identify the three different regions of different conductivity, the reconstruction of the interface between these regions is quite inaccurate because of smoothness enforced by the parameterisation used within EKI. We overcome this limitation by means of a level-set parameterisations within the EKI framework.

Level-set parameterisation and prior ensemble
For the series of experiment in this subsection we test EKI with the level-set parameterisation κ = P 2 (u) from (56) where, for simplicity, we use κ ι = λ ι , (ι ∈ {b, l, h}) constant. Hence the parameterisation (56) becomes The unknown consists of 8 scalars and one function ω f which we discretise also on a 100 × 100 grid. Hence, the dimension of the unknown u is dim(U) = 10, 008. We select the initial ensemble according to where we see the choices λ f = 1, ν f = 2 and σ f = 0.5 fixed across the ensemble. From (54) we know that this selection produces a centred ensemble of initial level-set functions. The selection for (fixed) λ f is made for simplicity; we expect to capture all the variability of the level-set function via the term W Θ f ω f in (54). It is worth noticing from (72) that there is no overlapping among the support of the distributions for the conductivity values (i.e. λ l , λ b and λ h ) on each region. Furthermore, each of these intervals contain the true values form (70). Therefore, we work under the assumption that (i) there is clear contrast between the (unknown) values on each region and (ii) we have good knowledge of the range of possible values for the (unknown) conductivity in each of those regions. In Figure 11 we show the plots of some log κ (j) 0 = log P 2 (u (j) 0 ) (top panels) and the corresponding level-set functions f (j) 0 = log(P 1 (u (j) f,0 )) (bottom panels) from five ensemble members. Our selection ν f = 2 imposes smoothness in the level-set function and, in turn, in the interface between the three regions of low, high and background conductivity. From Figure 11 (bottom) we notice that the ensemble of level-sets displays anisotropy induced by randomising the intrinsic lengthscales in the level-set function. This variability can be seen in the corresponding interface between regions of different conductivities (top row Figure 11). Note that the values of conductivity within each region are variable across particles.

Results from several choices of J in EKI-DMC.
In Figure 12 we show boxplots of the relative error w.r.t the truth as well as data misfits (65)- (67). As before, these are results from 30 EKI-DMC runs using different initial ensembles for each choice of J. We can clearly see a decrease in the error w.r.t the truth as we increase J, while the data misfits DM 1,n and DM 2,n achieve values close to the noise level (indicated via dotted red line) for J ≥ 200. We see that, again, J = 200 is a good compromise between accuracy and cost. Using J = 400 will double the cost with a marginal improvement in accuracy. The choice of J = 200 also yields reasonable visual results as we can verify from Figure 13 where, for one run of EKI-DMC, we display the log of κ n * = P 2 (u n * ) (top panels) and the level-set function f n * = log(P 1 (u f,n * )) (bottom panels).
Although we note that the average data misfit DM 3,n * in Figure 12 seems to increase with J, experiments (not shown) with larger J suggest that DM 3,n * eventually stabilises. The increase in DM 3,n * can be attributed to the fact that larger J's produces a better spread/coverage of the distribution of particles; some of these particles yield a larger data misfit within EKI 6 . The ensemble mean u n , however, is quite accurate (see Figure. 13, for J = 200) and so the corresponding κ n = P 2 (u n ) yields reasonable values of the DM 2,n * (i.e. close to the noise level).   Figure 14 shows some members from the final (converged) ensemble corresponding to the initial ensemble from Figure 11; the top panels are log κ f,n * )). This figure shows that there is, indeed, significant variability across particles of the converged ensemble κ (j) n * and, in turn, possible large spread in the values of the data misfit obtained using each of these conductivities (hence potentially large values of DM 3,n ). The average error w.r.t. the truth and (log) data misfit DM 1,n are shown in Figure 15, as a function of the iteration number n (these are results from our 30 runs). Note that, in contrast to the previous experiment, the error displays more variability across ensembles.
Plots of log κ n and the level-set function f n , at some of the intermediate iterations 1 ≤ n < n * , are displayed in the top and bottom panels of Figure 16, respectively. We can see that EKI not only estimates the shape (via the level-set function) of the regions with different conductivity but also the conductivity values on each region. In Figure 17 (top) we plot, as a function of n, the values of L 1,f , L 2,f , λ l , λ b and λ h , i.e. the means of the ensembles {L h } J j=1 (we reiterate that these variables are updated at each iteration of EKI). The ensemble mean for the intrinsic lengthscale of the leve-set function in the vertical direction is larger than the horizontal one. Again, this is due to the presence of the regions of low conductivity which are longer in the vertical direction. From the middle-right panel we can see that the true conductivity in the background region λ † b is recovered quite accurately. In contrast, the mean values λ l and λ h do not seem to vary much with respect to the mean.
Although we are mainly focus on the deterministic case here, to further appreciate the accuracy of the inversion for these variables, in Figure 17 (bottom) we show their probability densities approximated from the initial and final (converged) ensembles. For most of these variables we can see that the converged ensemble has a much smaller variance compared to the initial one. For the conductivity values we note that the lower and background values are identified accurately with the ensemble mean; the higher value is captured in the tail of the final ensemble.

Comparisons of EKI-DMC and EKI-LM
In Figure 18 we compare the performance of EKI-DMC and EKI-LM using 30 different initial ensembles. For EKI-LM we explore different choices of ρ. We note that EKI-DMC outperforms EKI-LM for all our choices of ρ. From Table 1 (for J = 200) note that for ρ > 0.7 the computational cost of EKI-LM is approximately twice the cost of EKI-DMC and the cost even triple if we choose ρ = 0.8. The plots of log κ n from one run with the three algorithms (same initial ensemble) are shown in Figure 19, where we see that all these runs perform well. Similarly conclusions to those in Exp EIT 1 are also drawn for this case. Namely, EKI-DMC is more accurate than EKI-LM in the chosen metrics although visually we achieve good performance from both algorithms and inputs. The main advantage of EKI-DMC over EKI-LM is that the former is robust as does not reply on tuning parameters, while for the latter, the computational cost can double or even triple with relatively small changes in the tuning parameters.  Figure 18: Exp EIT 2 . Error with respect to the truth (left), E n * (see (62)), and data misfit DM 1,n * (65) (right) computed at the final iteration n * using EKI-DMC and EKI-LM with various choices of ρ. The noise level estimated by δ = √ M is indicated with the dotted red-line in the right panel. Figure 19: Exp EIT 2 . Logarithm of κ n * computed using the same initial ensemble (J = 200) with EKI-LM using different selections of the parameter ρ. In the right panel with display the log of κ n * that we obtain using EKI-DMC.

Numerical Experiments for non-destructive testing of composites.
In order to solve the forward model described by (39) The plot of the (log) true κ † that we define for this experiment is displayed in Figure 20 (topleft). The permeability consist of the background (low perm.) region and defects of higher constant permeability with value λ † h = 10 −8 m. For the background region we produce heterogeneity via a realisation of a Gaussian random field with constant mean λ † l = 10 −9 m and Matern covariance with lengthscales L † 1,l = 0.65, L † 2,l = 0.05, and ν = 1. Note that we create anisotropy along the horizontal direction. We solve (39)- (45) for the κ † defined above and the plots of the pressure field that we obtain, at some of the observation times, are shown in the bottom panels of Figure 20. Synthetic data are computed by y = P † + η where P † = F(κ † ) (see eq. (47)) and η ∼ N (0, Γ) is 2% Gaussian noise. In other words we choose, Γ = diag(γ 1 , . . . , γ M ) with γ m = (0.02 P † m ) 2 .

Exp RTM 1 . Level-set parameterisation with constant lengthscale
We use the parameterisation κ = P 2 (u) from subsection 3.3.2 with only two regions: low/background and high permeability (i.e. ζ 1 = ζ 2 in (53)). We assume that the variability of permeability values in the region of high permeability (i.e. defects) is small so that we use κ h (x) = λ h (constant). The parameterisation in (53) reduces to where u = (u l , λ h , u f ), with u l = (λ l , ν l , σ l , L 1,l , L 2,l , ω l ), We note that u here consists of eleven scalars and two functions, ω l and ω f that we discretise on a 80 × 80 grid that we use for the WM parameterisations involved in P 2 . The dimension of the unknown u is dim(U) = 12811. We consider the following selection of initial ensemble: Plots of some of the (log) κ 0 = P 2 (u (j) 0 ) and the corresponding level-set function f 0 = log(P 1 (u f,0 )) are shown in the top and bottom panels of Figure 21. We use EKI-DMC with various choices of J and conduct 30 runs with different initial ensembles. For one of these runs, the log of our converged estimate κ n * = P 2 (u n * ) and its level-set function f n * = log(P 1 (u f,n * )) are displayed in the top and bottom panels of Figure 22, respectively. Boxplots of the error w.r.t the truth and data misfits from our 30 runs are displayed in the top panel of Figure 23. Although the quality of the identification seems to improve with increasing J, the average relative error w.r.t the truth is quite large (e.g. > 50% for J = 800). The number of iterations to achieve converge with EKI-DMC is reported in Table 1. Similar suboptimal results (not shown) are obtained using EKI-LM.
The results above show that P 2 , which uses constant lengthscales for the level-set function, cannot accurately identify the presence of the type of defects displayed in the true κ † . In order to further understand this limitation, let us inspect from Figure 24 the densities of the initial and final (converged) ensembles for the lengthscales of the level-set function (L 1,f and L 2,f ). These plots corresponds to one run for the case J = 400 shown in Figure 22. While the initial ensemble covers a very wide range of values, EKI produces ensembles which are fairly concentrated around the mean values L 1,f = 0.3 and L 2,f = 0.05. These values are consistent with what we observe in Figure 22, namely, that we can identify the thin horizontal defect centred around x 2 = 0.675 via a relatively large inferred L 1,f and a small inferred L 2,f , However those lengthscales are not consistent with the vertical defects (centred at x 1 = 0.275 and x 1 = 0.725) which could only be characterised via a level-set function with small L 1,f and large L 2,f . Note that the background permeability field κ l (x) displays higher values to compensate for the absence of a the high permeability region from the defects. For this experimental setting, the data are more informative of the horizontal channel which has substantial influence on the flow pattern.
In summary, κ = P 2 (u) is limited to cases in which we most defective areas will have relatively similar geometry such as those considered in subsection 4.2. Large ensemble sizes can partially alleviate this issue but the computational cost may become prohibitive. We aim at addressing this limitation by incorporating spatial variability in these intrinsic lengthscales of the level-set function.

Exp RTM 2 . Level-set parameterisation with spatially variable lengthscale
We repeat the same experiment from the previous subsection but we use the parameterisation κ = P 3 (u) discussed in subsection 3.3.3. We reiterate that using variable lengthscales WM fields for inverse problems was used in [56] albeit not in the context of level-set parameterisations. As before, we consider only two regions and assume the permeability in the anomalous region is constant but  unknown. The parameterisation in (53) simplifies to The unknown u consists of 19 scalars and 5 functions (i.e. ω f , ω l , ω L 1,f , ωL 2,f ) discretised, as before, on a 80 × 80 grid. The parameter space here is of dimension dim(U) = 25, 619 The initial ensemble is chosen as follows: where, for simplicity we keep fixed the lengthscales of each level-set function lengthscale (i.e. L 1,L i,f and L 2,L i,f ). The plot of (log) κ 0 = P 3 (u We conduct 30 runs with different initial ensembles using EKI-DMC with different choices of J. For one of these runs, in Figure 26 we show the log of κ n * = P 3 (u n * ) (top row), the level-set function f n * = log(P 1 (u f,n * )) (top-middle row) as well as the horizontal (middle-bottom row) and vertical lengthscale (bottom row) i.e. L i,f = P 1 (u L i,f ) (i = 1, 2), respectively. In contrast to the results we obtain from Exp RTM 1 , (See Figure 22), here the case J = 400 can accurately identify the defects from the truth κ † . Increasing to J = 800 produces relatively small improvements while doubling the cost of the inversion. Plots of κ n * at intermediate iterations are shown in Figure 27. From Figure 26 we notice (see middle-bottom row) that EKI infers large values for the horizontal lengthscale L 1,f (x) for points x in the upper half part of the domain. Approximately in the same region, the vertical lengthscale L 2,f (x) (bottom panels) takes quite small values. These values are consistent with the reconstruction of the horizontal defect in the upper part. Similarly, we note that the two vertical defects are identified because EKI infers a small (resp. large) horizontal (resp. vertical) lengthscale for the lower part of the domain.
In Figure 28 we show densities computed from the initial and final ensemble of (i) lengthscales of the background field (L 1,l ,L 2,l ), (i) mean of background permeability, λ l , and (iii) values of permeability on defects, λ h . In these plots we also include the corresponding values from the truth. We notice that the mean of the background field has been successfully identified by the final ensemble mean. In contrast, EKI over estimates the permeability value on the defects region. We notice that only the lengthscales of the background field in the horizontal direction, L 1,l , has been accurately recovered by the corresponding mean of the final ensemble.
In Figures 29 we compare the performance of EKI-DMC and EKI-LM over multiple runs. In Figure 30 we display plots of κ n * from one single run obtained with different choices of ρ in EKI-LM. In Table 1 we display the average number of iterations obtained with both methods. The results show that EKI-DMC and EKI-LM have, in average, the same level of accuracy when we chose ρ = 0.8 in EKI-LM. However, the cost of the latter is 2.6 times the cost of the former. Again, we see that EKI-DMC is robust and computationally efficient when used for the estimation of physical properties that parameterised in highly-complex manner.

Conclusions
We introduced the data misfit controller (DMC): a new adaptive regularisation strategy within the classical EKI setting. This led to an algorithm, EKI-DMC, that in contrast to existing EKI ap-  Logarithm of κ n * computed using the same initial ensemble (J = 400) with EKI-LM using different selections of the parameter ρ. In the right panel with display the log of κ n * that we obtain using EKI-DMC.
proaches, it does not require any tuning parameters. Although we focus on the solution of deterministic identification problem, the proposed DMC is motivated from the Bayesian perspective of EKI within the tempering setting, where the inverse of the regularisation parameter, α −1 n , controls the transition between two consecutive intermediate measures. The Bayesian tempering setting provides us with a condition that these parameters must be satisfy ( n * n=1 α −1 n = 1) to bridge the prior and posterior. We encode this condition for the termination of EKI-DMC together with our new method for choosing α n .
We applied EKI-DMC for the solution of two PDE-constrained inverse problems: (i) electrical impedance tomography with the complete electrode model and (ii) estimating permeability of a com-posite preform from pressure collected during resin injection. In these problems, the unknown physical property was parameterised via different maps that enabled us, via the EKI framework applied to the appropriate parameters, to characterise smooth functions and piece-wise smooth physical properties. Wittle-Matern fields were at the core of these parameterisations which included the intrinsic lengthscales as inputs that we estimate within EKI. For the piece-wise smooth, we use a truncated WF field (the level-set function) to characterise discontinuities between different regions. Variability on each region is also incorporated via WM fields. Our experimental settings also explored WM parameterisation of heterogenous lengthscales for the level-set function.
As with any other EKI method, our results show that the performance of EKI-DMC relies on reasonable choices of ensemble size, J. For sufficiently large choices of J, our experiments show that EKI-DMC is quite robust and capable of producing accurate identification of physical properties parameterised in the appropriate manner. We performed a performance comparison between EKI-DMC and the EKI-LM approach of [2]. In most cases EKI-DMC outperforms the accuracy of EKI-LM in terms of error w.r.t the truth and data misfit. Even when the performance of both algorithms was comparable (for certain choices of tuning parameters in EKI-LM), the computational cost of EKI-DMC was substantially smaller. The advantage in terms of computational cost was even more outstanding in the case of more challenging settings where we inferred unknowns parameterised via highly-complex mappings. In these settings, not only larger ensembles abut also more iterations are required to achieve convergence of EKI. EKI-DMC was quite robust and managed to produce accurate and stable results with reasonable computational resources.
While our numerical results are focused on the comparison against EKI-LM, we recognise that other choices of α n , including those discussed in Section 2, could also perform in a robust and accurate fashion. However, most of these approaches rely on the choice of additional tuning parameters. Of course, any of these other methodologies may display a comparable performance to EKI-DMC when those tuning parameters are selected optimally. Such a an optimal selection, however, can only be found after a careful numerical investigation on the given problem-specific setting. We reiterate, that the aim of EKI-DMC is to provide EKI with a robust self-tuning regularisation strategy that can be used in practical and operational identification problems for which a optimal choices of tuning parameters may not be computationally feasible. Our results show that this aim has been successfully achieved via the proposed framework.

Remark Appendix B.2 (Squared-root EKI).
It is worth noticing that other approaches can be used to approximate N (m n+1 ,C n+1 ) above. This includes the so-called ensemble square-root formulations [67] in which the particles are cleverly updated so that their sample mean and covariance coincide (exactly) withm n+1 andC n+1 . While these has been shown to be beneficial for very small samples (i.e. < 50), we note that N (m n+1 ,C n+1 ) does not, in general, coincides with ν n+1 = N (m n+1 , C n+1 ) (unless G is linear). Remark Appendix B.3 (EKI as a derivative-free approximation of LM). For sufficiently large J, the mean of the ensemble u n+1 approximatesm n+1 and so m n+1 which is, in turn, the iteration of the LM scheme in (5) (see Remark Appendix B.3). Therefore, we can interpret EKI as a derivativefree approximation of the LM scheme constrained to the subspace generated by the initial ensemble {u (j) 0 } J j=1 . More specifically, the ensemble mean u n define by the recursive formula (4) satisfies the following subspace invariance property (see Theorem 2.1 in [1]) for all n ∈ N. We expect EKI to produce approximate solutions to (3) within the subspace defined above. While the numerical experiments from [2] provides evidence of such a claim, to the best of our knowledge, the convergence of EKI in this context is still an open problem.