Target-enclosing inversion using an interferometric objective function

Full waveform inversion is a high-resolution subsurface imaging technique, in which full seismic waveforms are used to infer subsurface physical properties. We present a novel, target-enclosing, full-waveform inversion framework based on an interferometric objective function. This objective function exploits the equivalence between the convolution and correlation representation formulas, using data from a closed boundary around the target area of interest. Because such equivalence is violated when the knowledge of the enclosed medium is incorrect, we propose to minimize the mismatch between the wavefields independently reconstructed by the two representation formulas. The proposed method requires only kinematic knowledge of the subsurface model, specifically the overburden for redatuming, and does not require prior knowledge of the model below the target area. In this sense it is truly local: sensitive only to the medium parameters within the chosen target, with no assumptions about the medium or scattering regime outside the target. We present the theoretical framework and derive the gradient of the new objective function via the adjoint-state method and apply it to a synthetic example with exactly redatumed wavefields. A comparison with FWI of surface data and target-oriented FWI based on the convolution representation theorem only shows the superiority of our method both in terms of the quality of target recovery and reduction in computational cost.

hand, data-driven approaches that rely for example on Marchenko redatuming (e.g., van der Neut et al. 2015) only require a kinematically accurate, smooth model of the overburden. Moreover, despite early successful results, all of the aforementioned methods apart from that of Cui et al. (2020) are not fully target-oriented as they require estimating model parameters, to a relatively high degree of accuracy, in an infinite half-space from the datum of interest in the subsurface.
We propose a new full waveform inversion (FWI) framework to invert for the velocity model locally from full waveforms redatumed to an enclosing boundary of the subsurface subdomain of interest based on the objective function initially proposed by Vasconcelos et al. (2016). This objective function is conceptually different from both the traditional FWI misfit function and those used in previous approaches to target-oriented inversion: instead of relying on the equivalence of the source and receiver wavefields at the physical (or virtual) receiver locations, it builds on the equivalence of the wavefields reconstructed from the boundary data by the convolution and correlation representation formulas anywhere in the local subdomain. Examples of such representation formulas can be found, in e.g. classical texts such as Morse (1953) (acoustic) and Aki & Richards (2002) (elastic).
The idea of the proposed method stems from the fact that in a lossless medium, both the convolution and correlation representation formulas can be used to reconstruct the true wavefields in an open subdomain from wavefield measurements on the subdomain boundary, if the Green's function, i.e. model parameters within the subdomain are exactly known. Moreover, the wavefields reconstructed by the convolution and correlation representation theorems are the same up to the time direction. If the Green's function, i.e. the model parameters are incorrectly known, then these fields have similar features but differ in details. Our inversion scheme is therefore driven by the mismatch of the wavefields reconstructed by the convolution and correlation representation formulas when the model is incorrect.
In the ideal situation, when the boundary data are exact, the mismatch between the convolution and correlation wavefields is nil when all the features of the true model are recovered by the inversion.
Thus, the proposed method represents an ideal platform for local refinement of the subsurface model parameter distribution.
Moreover, as the boundary data contain all of the waves entering, leaving and re-entering the local domain, all of these wave phases are utilized by the inversion. Also, since the convolution and correlation representation formulas reconstruct the wavefields only within the injection boundary, the method is insensitive to model inaccuracies outside of the injection surface. Therefore, the proposed inversion method does not require us to invert for an infinite half-space below the top boundary. The modelled wavefields include the waves entering the local subdomain from below without the necessity to model wave propagation in the underburden at each step of the inversion. In other words, the local domain forward modelling process accounts for all events coming from outside of the local subdomain.
The proposed method requires the knowledge of the wavefields on a closed boundary completely surrounding the local subdomain. In practice such data can either be measured by borehole receivers located underground, or, more often, they are obtained by redatuming methods, e.g., Marchenko redatuming (Wapenaar et al. 2014;van der Neut et al. 2015;Ravasi 2017;Vargas et al. 2021). Within the context of local waveform inversion based on redatumed wavefields and the convolution representation formula, Cui et al. (2020) quantitatively assess the retrieval accuracy of the Marchenko method.
They show that inversion with redatumed boundary data is almost as accurate as the inversion with the exact boundary data, although, admittedly, they use a starting velocity model of a very high quality both for redatuming and local inversion. They also show that a local refinement of the subsurface model with the redatumed wavefields is comparable in resolution to the refinement obtained by full waveform inversion of surface data.
The main contribution of this work is to formulate a new full waveform inversion method based on the interferometric objective function of Vasconcelos et al. (2016) and investigate its properties. We formulate the interferometric objective function and the partial differential equation (PDE) constraints stemming from the convolution and correlation representation integrals, using the constant-density vector-acoustic formulation for pressure and particle displacement (Fleury & Vasconcelos 2013;Zheglova & Malcolm 2020), where the inverted model parameter is squared slowness. We derive the gradient of the interferometric objective function with respect to the model parameter by the adjoint state method. Then we implement an interferometric full waveform inversion (IFWI) method, where the model updates are produced with the help of an L-BFGS optimization engine. We investigate the properties of the interferometric objective function, its gradient and the new inversion method on a series of stylized synthetic examples, using exact (i.e., directly modelled) boundary wavefields. This choice allows us to focus on the features and effectiveness of the proposed approach in the ideal situa-Target-enclosing inversion using an interferometric objective function 5 tion. We show that with exact data, the proposed method is able to better recover the low wavenumber components in the model compared to the conventional surface FWI method and the local FWI method based on the convolution representation formula (Cui et al. 2020). Then we estimate with proxy examples the potential influence of using redatumed data. In particular, we investigate the influence on the inversion performance of kinematically incorrect boundary data, and the influence of missing data on the inversion result. Finally, we discuss potential outcomes and issues related to the implementation of the method with redatumed wavefields.

Preliminaries
In this section we review the vector-acoustic convolution and correlation representation formulas that will be used in the formulation of the interferometric full waveform inversion (IFWI) objective function. As shown by , equations (63) and (67), the vector-acoustic field w = where p is pressure and v is particle velocity, in an open subdomain D of the subsurface can be exactly reconstructed in forward or reverse time order from measurements on the boundary ∂D, if the Green's function G D in D is known exactly, and the physical source of the wavefield is located outside of the subdomain D. Denoting the field reconstructed forward in time as the convolution field, w conv , and the field reconstructed in reverse time order as the correlation field, w corr , Wapenaar's equations are written in the frequency domain as: where x is a point in D; x r is a point on the boundary ∂D; x s is physical source of the field, located outside of D ∪ ∂D; ω is angular frequency; N r is the matrix containing the outward normal components n to ∂D: Figure 1. Schematic of wavefield extrapolation: (a) convolution, (b) correlation. The wavefield from the source x s is propagated to a point x r ∈ ∂D, from which it is reconstructed at x ∈ D using equations (1) and (2).
K is the diagonal matrix encoding the change of sign in v due to time-reversal: Here, * is the complex conjugation in the frequency domain corresponding to time-reversal in the time domain, whilst 0 and I are the zero and identity matrices, respectively.
Schematically this process is shown in Figure 1. The wavefield from the source x s / ∈ D ∪ ∂D propagates to points x r ∈ ∂D, from which it is reconstructed at x ∈ D using equations (1) and (2).
Assume that w is propagated to ∂D without error, so that w(x r ) is exact. If the Green's function G D is also known exactly, i.e. it is equal to the true Green's function G true then the convolution and correlation representation formulas (1) and (2) reconstruct the same wavefield, which is also the true wavefield w. If G D is not known exactly, i.e. G D = G true , then both the convolution and correlation wavefields differ from the true field w and from each other. Writing this as equations, we have that except for pathological cases of symmetry that that are extremely unlikely to occur in real life situa- tions. This is also demonstrated numerically by Vasconcelos et al. (2016). In practice, if the fields are redatumed to x r using e.g. Marchenko redatuming, they may contain errors such as missing and nonphysical events (van der Neut et al. 2015), as well as errors due to incorrect kinematics of the macro velocity model used for redatuming. In this case, equations (5) and (6) may hold only approximately.
Target-enclosing inversion using an interferometric objective function 7

IFWI problem formulation
The above set-up leads to an FWI problem formulation with an objective function being the L 2 norm of the difference between w corr and w conv , where the convolution and correlation fields satisfy constraints (1), (2).

PDE constraints
We set up the inverse problem in the constant density acoustic formulation for pressure p and scaled Integral equations (1), (2) can be expressed as partial differential equations for pressure and scaled displacement using the forward and adjoint vector-acoustic differential operators (Zheglova & Malcolm 2020) and convolution and correlation areal sources along the boundary ∂D: where: corr} are the convolution and correlation pressure and particle displacement vector-acoustic wavefields; L and L † are the vector-acoustic differential operator and its adjoint: s conv and s corr are the convolution and correlation sources given at each x r ∈ ∂D as δ ∂D (x, x r ) denotes the areal source distributed along the boundary ∂D.
As per the time-domain adjoint-state approach, Equation (7) is solved forward in time, while equation (8) is solved backward in time.

Objective function
Based on the reciprocity relations we discuss above, the objective function we use is where ||x|| 2 2 = i x 2 i is the squared Euclidean l 2 norm and Λ is the data weighting operator. In this study we use which samples only pressure. Note that here, unlike the representations in Eqs. 1 and 2, we use an explicit time-domain representation for the sake of consistency with our time-domain implementation of the adjoint-state vector-acoustic wave equations.
Equation (9) together with equations (7) and (8) constitute the PDE constrained optimization problem that we solve for the model parameter m.
We remark that the convolution and correlation differential equations (7) and (8) together constitute the forward problem that is exactly satisfied by the convolution and correlation fields. The (redatumed) data w(x r , x s , t) at the boundary ∂D constitute the observed data. The convolution and correlation sources are fully determined by the redatumed data and the geometry of ∂D. Unlike in conventional surface FWI and local FWI, we do not directly compare the observed and modelled fields along a line of receivers. Rather, the modelled data are the pressures p corr = Λw corr and p conv = Λw conv , whose difference must be zero everywhere in the local domain. This process can also be expressed as minimization of Λ −Λ w corr w conv T 2 2 . The implications of such an objective function are far wider than it may initially sound: by imposing equivalence between two wavefields, neither of them must be known directly at the location where the objective function in evaluated, allowing us to consider a grid within the entire domain of interest. When only one representation theorem is used (e.g. Cui et al. 2020) the objective function can also be theoretically evaluated at any point in the grid; however, in practical application this would require being able to access the wavefield at such location by means of e.g., Marchenko redatuming, leading to a radical increase in the computational cost associated with the data preparation prior to inversion.

IFWI gradient
The gradient of the interferometric objective function with respect to the model parameter m is derived by the adjoint state method (Plessix 2006;Fichtner 2011). We derive it in Appendix A, and show here only the result. In the general form, the gradient of the interferometric objective function with Target-enclosing inversion using an interferometric objective function 9 respect to the squared slowness m is given by  are the adjoint convolution and correlations wavefields that satisfy the adjoint convolution and correlation equations: Equations (11) and (12) are solved respectively backward and forward in time.
With the choice of Λ as in (10), the gradient becomes Since the residual w corr − w conv is calculated over the whole D, the convolution and correlation adjoint sources, i.e. right-hand sides of equations (11) and (12)  We implement the iterative inversion by combining the IFWI objective and gradient with the L-BFGS optimization method.

Convexity of the interferometric objective function
In this section, we investigate the convexity of the interferometric objective function in the target domain using exact data, and compare it to the surface FWI objective function. We consider an example, where we obtain the true model, shown in Figure 2, by applying a mask to a modified Marmousi model. We place six sources between 1 and 3 km in the horizontal direction at a depth of 0.1 km and use a Ricker wavelet with a peak frequency of 15 Hz as the source signature.
Because a kinematically accurate macro velocity model is usually needed for both full waveform inversion and redatuming, we are primarily interested in the response of the objective function to the lack of high wavenumber components in the model. Therefore, we apply 2D Gaussian smoothing to the target area, where we vary the standard deviation of the Gaussian filter from 5 to 100 m to obtain progressively smoother models of the target, and compute the objective function value for each of the smoothed models. Examples of such smoothed local models are shown in Figure 2 (b), (c). Figure 3 shows the interferometric objective function variation with the STD of the Gausian filter. the calculation of FWI and VAFWI objective function, we apply no smoothing to the overburden layer. Note that, apart from the difference in absolute values, the interferometric objective function is more convex than both the FWI and VAFWI objective functions, in addition, it has a wider basin of attraction near the global minimum. Based on this experiment, we can expect the IFWI to exhibit fast initial convergence with a subsequent slow-down near the minimum. We can also expect the IFWI to be more resistant to errors in the initial model, which we show to be case with exact data.

Detailed analysis of the interferometric gradient
In this section, we take a closer look at the gradient of the interferometric objective function. First, we observe that each of the two terms in the gradient (equation 13) is a zero-lag correlation of the forward field with the corresponding adjoint field. Each of the adjoint fields in turn consists of two parts: a part due to the convolution forward field, and a part due to the correlation forward field, as follows from the adjoint equations (11) and (12). We denote these adjoint field components as follows: (i) w conv † conv = Γ † Λ † ΛΓw conv is the part of the convolution adjoint field w conv † due to the convolution part of the adjoint source Λ † Λw conv (ii) w conv † corr = −Γ † Λ † ΛΓ † w corr is the part of the convolution adjoint field w conv † due to the correlation part of the adjoint source −Λ † Λw corr (iii) w corr † corr = ΓΛ † ΛΓ † w corr is the part of the correlation adjoint field w corr † due to the correlation part of the adjoint source Λ † Λw corr , and (iv) w corr † conv = ΓΛ † ΛΓw corr is the part of the correlation adjoint field w corr † due to the convolution part of the adjoint source −Λ † Λw conv , where Γ and Γ † are Green's functions for operators L and L † . Correspondingly, the four terms in the gradient are: ∂I ∂m ∂I ∂m We use two starting models to analyze the behaviour of the IFWI gradient, shown in Figure 2 (b) and (c). The model in Figure 2 (b) represents late inversion stages, while the model in Figure 2 (c) represents early inversion stages. The original field is generated by a Ricker wavelet source with a peak frequency of 15 Hz located at the middle top of the model (Figure 2) at z = 0.1 km, x = 2 km.
First, we look at the gradient obtained in the starting model in Figure 2  smooth and lacks any high-wavenumber details, so the convolution forward field does not undergo any scattering inside the local domain and is purely down-going: p conv ↓. It also generates a purely down-going adjoint field p conv † conv ↓ and a purely up-going field p corr † conv ↑. Zero-lag correlation of the fields traveling in the same direction generates low-wavenumber updates, whereas zero-lag correlation of the fields travelling in opposite directions generates low-wavenumber updates. The term ∂I ∂m (i) is a cross-correlation of two down going waves, and therefore updates only low wavenumbers, Figure   5 (a). The correlation field p corr has both transmitted and reflected waves from the true model that are respectively down-and up-going: p corr ↓ + ↑, and generates up-and down-going waves in the adjoint field components p conv † corr :↑ + ↓ and p corr † corr :↑ + ↓.

Target-enclosing inversion using an interferometric objective function 13
The second observation we make from Figures 4 (e), (f) and 5 (e), (f) is that the convolution and correlation components in the gradient contribute to different parts of the model update. The convolution term T 0 p conv p conv † dt, Figure 4 and 5 (c), updates the target area mostly above the bottom reflector of the target. The correlation term T 0 p corr p corr † dt, Figure 4 and 5 (f), updates the target area mostly from the top reflector to the bottom of the injection boundary ∂D.
Figure 6 (a) shows the gradient update, computed in the local subdomain by IFWI using 5 equally spaced sources located at z = 0.1 km and x = 1.2 to 2.8 km in the initial model in Figure 2 The update by the local FWI method of Cui et al. (2020) is shown in Figure 6 (b) for comparison. It is interesting to note here that the local FWI method of Cui et al. (2020) uses only the convolution field p conv and its adjoint computed from the misfit between the convolution field and redatumed data at the virtual receivers near the top boundary (green line in Figure 2). It is impossible to generate low wavenumber updates with only a purely convolution field and its adjoint, if there are no upgoing waves present in the forward convolution field. We observe that the update in Figure 2  We note also that even though the low wavenumber update appears possible with IFWI, a kinematically correct macro model is still important for redatuming, as an incorrect macro model of the target causes kinematic errors in the redatumed boundary data that propagate into the reconstruction.
We investigate this issue in the Examples section with a proxy example.

Computational cost
In this section we express the computational cost of the proposed method and compare it with that of the other FWI methods used in this paper as comparison. First, we consider the cost of IFWI. We use capital letters for full domain and surface quantities, and small letters for local subdomain quantities. In the implementation of IFWI, the computational cost is dominated by the forward and adjoint modelling. The cost of one constant-density vector-acoustic PDE solve in flops is proportional to n x n z n t in 2D and to n x n y n z n t in 3D, where n x , n y and n z are the dimensions of the local subdomain in grid-points and n t is the number of time steps. A total of four PDE solves on the local subdomain are required per source for one IFWI gradient evaluation: two forward solves, equations (7) and (8) order accurate discretization in time, and PML width of 20 points, C ≈ 100 for our current 2D implementation. This constant can be reduced at the expense of the code clarity and a small increase in memory use. Assuming that n x , n y , n z ∼ n, the cost of one iteration of IFWI is ∼ O(n 2 n t N s ) in 2D and ∼ O(n 3 n t N s ) in 3D.
The cost of redatuming, which only needs to be performed once, must be added to the cost of waveform inversion. Limiting ourselves here to data-driven redatuming, and more specifically Marchenko Cost (Red) = C red n red n iter 2[(2 n t log(n t ) N r + n t (4 N s (2 N r − 1))) + n t N r ]

Target-enclosing inversion using an interferometric objective function 17
where n red is the number of points in the subsurface where we need to redatum the fields n iter is the number of iterations of the Marchenko equations per one point N r is the number of surface receivers C red is an unknown proportionality constant.
For the terms in the square brackets, the first corresponds to forward and inverse FFTs, the second term is the discretized integral over the acquisition surface in complex-number arithmetic, and the third term is the application of muting. The factor of 2 in front of the square brackets arises from the fact that two multi-dimensional convolutions are required per iteration of redatuming. Finally, n red ∼ 4(n x + n z ) in 2D and n red ∼ 4(n x n z + n x n y + n y n z ) in 3D. This accounts for two horizontal and two vertical subdomain boundaries (four in 3D), to each of which pressure measurements need to be redatumed at two layers of points in order to obtain displacement components/pressure derivatives (Cui et al. 2020).
In order to compare the computational cost of local FWI and redatuming, we observe that the cost of redatuming is dominated by n t and N r . We assume N r to be in the order of O(N x ) in 2D and O(N x N y ) in 3D, where N x , N y represent the dimensions of the full domain in grid-points. Moreover, if n x , n y , n z ∼ n, and N x , N y ∼ N then redatuming the surface data to the whole local domain boundary has complexity ∼ O(nN n t N s ) in 2D, which is comparable to the order of complexity of one iteration of IFWI, and in practice is likely to have the cost of a few iterations of IFWI due to other multiplicative factors in equation 18 and the fact that N > n. In 3D, the cost of redatuming is O(n 2 N 2 n t n s ), which is larger by a factor of N than one iteration of FWI.
Conventional surface FWI requires two PDE solves on the full domain per source: one forward and one adjoint solve, and has the cost ∼ O(N 2 n t N s ) in 2D and ∼ O(N 3 n t N s ) in 3D per iteration.
If V f ull = N x (N y )N z and V local = n x (n y )n z denote the volume of the full model and the volume of the local subdomain in grid-nodes, and they are discretized using the same temporal grid, then the cost of IFWI roughly becomes: Therefore, the IFWI method is practical in 2D compared to the conventional surface FWI when the local domain occupies less than roughly 1/2 of the full model. Due to the cost of redatuming in 3D, 1/2 reduction in the domain size might not be sufficient to offset the redatuming cost. In local domain applications the goal is to make V local V f ull as small as possible. In the proposed implementation, since the side injection boundaries are included in the forward modelling process, the local subdomain size can be reduced both in width and height.
When compared to the local convolution-based FWI of Cui et al. (2020), IFWI is admittedly twice as expensive per gradient evaluation for the same local domain size, number of sources and discretization, since local FWI requires only two PDE solves per source per gradient evaluation. Local The cost estimates for the three inversion methods are summarized in Table 1.

Exact redatuming
In this section, we demonstrate the performance of the interferometric FWI on stylized examples with exact redatuming. We compare performance of our method to the performance of surface FWI and the local convolution-based FWI of (Cui et al. 2020) (local FWI). For the first example, the true velocity model is the middle part of the velocity model shown in Figure 2 (a) between x = 1 and x = 3 km to reduce the cost of full domain modelling and inversion.
For the second example, we add a reflector below the target. The true model with reflector is shown in Figure 7, we also use the middle part of it between x = 1 and x = 3 km.
We place five evenly spaced sources at depth z = 0.1 km from x = 0.2 to x = 1.8 km. For the conventional FWI, we use the receivers at a depth of z = 0.1 km from x = 0.1 to x = 1.9 km. As before, the local domain is marked by the red line, the white line denotes the injection boundary ∂D, and the green line denotes the receivers used in the local FWI method.
For the proposed IFWI method and the local FWI, we start with the initial model shown in Figure   2 (c). For surface FWI, we additionally apply smoothing by Gaussian filter with STD of 20 m and 80 m to the overburden and the bottom reflector.
For the local FWI we do not apply any preconditioning, while for FWI we suppress the receiver Target-enclosing inversion using an interferometric objective function 19 footprint during the inversion, since in practice, receivers are likely to be in the water. We note that suppressing the receiver footprint for local FWI does not change the result significantly. We stop all inversions after 200 iterations, regardless of whether any in-built stopping criterion of the L-BFGS solver was reached.
3.1.1 Example 1: exact redatuming, no bottom reflector    below, such misfit measure significantly improves convergence and resolution of that method, but also dramatically raises the redatuming cost compared to the proposed IFWI method.
In Figure 11 (d), we show the inversion result obtained by the method of Cui et al. (2020) where the misfit is also computed everywhere inside the local subdomain just as for our proposed IFWI method. Figure 12 shows the objective function and the RMS model error plot for this inversion as dash-dot lines. With the data measurements everywhere, the local FWI converges in about 50 iterations to a model comparable in resolution to the IFWI recovery achieved after about 150 iterations. This is not surprising, since with full data available everywhere in the local subdomain, the convolution-based inverse problem becomes the problem of solving an overdetermined linear system where Lw and s are known everywhere. So it is a linear inverse problem, while IFWI remains a non-linear problem. Stated differently, the local FWI compares the modelled convolution field to the "known true" field everywhere, while IFWI compares two modelled fields to each other without access to the "true field".
Clearly, the ability to evaluate the objective function at every point in the local subdomain has a substantial impact on image resolution. However, the need to have pressure measurements at every cost of redatuming for that method from O(nN n t N s ) to O(n 2 N n t n s ) in 2D and from O(n 2 N 2 n t n s ) to O(n 3 N 2 n t n s ) in 3D. While doing this is still possible, this extra redatuming cost is likely to offset in 2D and exceed in 3D the extra cost per iteration and slower convergence of the proposed IFWI method. On the other hand, in the IFWI method, the correlation forward and adjoint fields are obtained at the cost of two additional local domain modelings per iteration. Thus, the same result is obtained by IFWI at a similar or lower cost.

Inexact and missing data
3.2.1 Example 3: inexact redatuming the data at the bottom and side boundaries, while the top boundary data is still exact. We invert this data using the biased initial model. Figure 14 shows the pressure field on the boundary used for this example and pressure data error.
As expected, the error is the largest at the bottom boundary and increases with depth at the side boundaries. Figure 13 (b) shows the IFWI reconstruction after 109 iterations. After this the image is still reasonable but becomes slightly grainy. We observe that the inversion is able to reproduce most of the high wavenumber details in the model, although some of the low wavenumber bias propagates into the reconstruction, particularly, the low velocity in the right and left bottom corners.
While some insight can be gleaned from this proxy example, a more accurate assessment needs to be made with redatumed data.

Example 4: missing data
In (ii) Low wavenumber updates: whilst the ability of the interferometric objective function to retrieve low wavenumber components in a target area is quite notable, this heavily depends on kinematic fea-tures captured by transmissions in the data (mostly along the bottom boundary). Therefore, the applicability of the proposed method for the macro model update may be limited for two reasons: (i) an accurate macro velocity model is still required to obtain kinematically correct redatumed data at the boundaries, and (ii) the kinematic errors in the boundary data may well lead to incorrect retrieval of the low wavenumber components in the model by inversion. The proxy example in section 3.2.1 shows that if kinematic errors from the target area are present in the data at the side and bottom boundaries, they propagate into the reconstruction. Therefore, in its current implementation, the proposed method is more suitable to localized model refinement rather than to macro velocity model building. An alternative update strategy may be necessary to mitigate the incorrect low wavenumber components in the redatumed data. Furthermore, more research is needed to assess the influence of redatumed-waveform errors caused by the overburden versus those related to the target medium, as we expect these to affect our method in different ways.
(iii) Distribution and density of volume points: the proposed objective function can be evaluated at any point inside the enclosing boundary. This plays the counterpart of physical receivers in conventional and local FWI. However, whilst the distribution and number of receivers is fixed and dictated by physical (in FWI) or algorithmic (in local FWI) constraints, the choice of the grid of points where the cost function is evaluated is totally arbitrary in our case. Our current implementation relies on a dense grid where the spatial sampling equates to that of the FD grid used for modelling of the wavefields.
The use of coarser (or finer), and possibly iteration dependent grids e.g., to balance bandwidth-related resolution with memory usage), will be a subject of future studies. For example, we envisage the combining low-frequency, coarse-grid FWI for redatuming-ready background and overburden models, with wide-bandwidth, denser-grid IFWI target inversion.
(iv) Cost function evaluation outside of the enclosing boundary: an additional feature of the proposed objective function is that it can also be evaluated outside of the enclosing boundary. In this case, it is not only the difference between the convolution and correlation representation theorems that will be zero in the presence of a correct model, but also each individual term alone must vanish as a result of writing the convolution and correlation representation formulas when both sources are outside of the boundary. This property is however not satisfied when the Green's function G D is modelled in the incorrect medium. Future research will investigate the benefit of enriching the interferometric objective function with grid points outside of the boundary.
(v) Multi-parameter and elastic inversion: whilst the numerical example presented in this paper targets only one parameter (i.e., velocity), the proposed framework is suited to multi-parameter inversion, namely (visco)elasticity and density, in that it relies on representation theorems that contain and reconstruct both pressure and particle velocity recordings -or to include particle velocity and stress Target-enclosing inversion using an interferometric objective function 27 fields in the more general cases. Moreover, extension to elastic media is also straightforward using the (visco)elastic counterparts of the convolution and correlation representation theorems (Aki & Richards 2002) -with the inclusion of viscoelasticity being also possible, granted the representation terms then include the appropriate volume terms. Moreover, given the much-increased computational costs associated with elastic media, our target-oriented approach may prove an essential tool to achieve detailed inverted models of target reservoirs at depth.

CONCLUSIONS
We propose a novel, target-enclosing, waveform inversion method based on an interferometric objective function that minimizes the difference between wavefields reconstructed from the boundary data by convolution and correlation representation formulas. The representation formulas are reformulated as PDEs constraints for the inversion. The method is formulated using full vector-acoustic boundary data consisting of pressure and scaled particle displacement. We derive the gradient of the objective with respect to the model parameter and implement the inversion using L-BFGS optimization engine. The method is shown to be very suitable for high-resolution model refinement, provided that an accurate macro velocity model is available.
In this initial work, we test the new inversion algorithm on stylized examples. We find that the objective function is also able to recover the low wavenumbers in the model better than other waveform inversion techniques when the boundary data is exact, i.e. contains the correct information about the macro velocity model. This seems to affect the resolution. In practice, the quality of the reconstruction also depends on the accuracy of redatuming, which in turn relies on the kinematic velocity model. Therefore, the ability of IFWI to recover the long wavenumber components is likely to be limited in practice. From the numerical examples, a key feature of the method appears to be evaluation of the objective function everywhere in the local subdomain. This is achieved at the expense of two extra local domain PDE solves per iteration, compared to the closest competitor method: the local convolution-based FWI method of (Cui et al. 2020). Considering that comparable resolution for the local FWI method is likely to require extra redatuming effort, our proposed method appears to be at least competitive in cost.
Further assessment of the robustness of the proposed technique with respect to redatuming is necessary. Future work will be focused on combining the proposed objective function with boundary data redatumed by means of the Marchenko method both on synthetic and field data.

APPENDIX A: IFWI GRADIENT DERIVATION
We begin by defining the convolution and correlation forward modelling operators: where M and D are the model and data spaces respectively, and the mapping is computed by solving equations (7) and (8) .