FD-injection-based elastic waveﬁeld separation for open and closed conﬁgurations

An important step in the processing of seismic data that are recorded at the free surface is the isolation of the primary incident waveﬁeld from the total recorded waveﬁeld (which is contaminated with the immediate reﬂections off the free surface). We present a 3-D waveﬁeld reconstruction technique, based on numerical waveﬁeld injection along a closed boundary, that allows us to isolate this primary waveﬁeld from measurements at the free surface. The technique consists of injecting only the three-component particle velocity recordings acquired at the free surface into a numerical waveﬁeld simulation, and additionally requires information about the medium properties. The result of our proposed procedure is the separation of elastic waves into their ﬁrst-order incident and reﬂected constituents, even when the recording or injection surface has sharp corners. With the use of synthetic data it is shown that the method achieves close to numerically exact waveﬁeld separation, provided that the true elastic model in the interior is used. In practice, the parameters for a homogeneous elastic model can be determined efﬁciently from the surface data itself using an optimization scheme. Finally, the waveﬁeld separation technique is successfully applied to experimental data, with particle velocity recordings acquired along ﬁve faces of a cubic granite rock volume. In addition to characterizing materials in laboratories, the presented technique has applications in numerical modelling and in so-called immersive experimentation, where the incident ﬁeld is required to immerse an elastic object in an arbitrary larger, virtual elastic environment. 20 kHz central frequency along the y -direction. The LDV used to the waveﬁeld individual cm ﬁve 8 and the top 4 corners of the rock. The bottom face, is the XY slice at the of the -axis with the for the We made 100-ms-long recordings with a 250 kHz sampling rate resulting in a temporal resolution of (cid:2) t LDV = 4 μ s. Additionally, the data were bandpass ﬁltered between 2 and 50 kHz. The signal-to-noise ratio of the


Elastic wavefield separation for open and closed configurations
1647 separate wavefield data recorded by densely deployed receivers at the Earth's free surface into its up-and downgoing components. Similarly, Van Renterghem et al. (2018) combine translational measurements and spatial wavefield gradient estimates to separate multicomponent land and seabed seismic data. Sollberger et al. (2017) utilize polarization analysis of seismic, 6-component data (3C particle velocity and 3C rotation data), recorded at the surface to identify wave type and propagation direction. Most of the methods described require transformations from the temporal to the frequency domain (Robertsson & Curtis 2002) or use temporal windowing of the data to isolate individual events (Sollberger et al. 2017).
Different wavefield separation methods are applicable when one has access to the domain interior, that is when recordings inside the target volume are available. One such class of methods is based on the principles of finite-difference (FD) wavefield injection introduced by Robertsson & Chapman (2000). In this paper, we present a modification of such a method. Wavefield separation enabled by FD-injection is based on a series of discoveries made by Blanch (2012), Ravasi & Curtis (2013), Vasconcelos (2013) and : when a surface of multicomponent particle velocity recordings is injected into seismic FD wave simulations, then the upgoing part of the wavefield in the recordings will propagate purely up from the injection surface in the FD simulation; simultaneously, the downgoing wavefield will propagate only in the downwards direction. Hence, FD-injection provides a method to decompose wavefields recorded at a single level into their up-and downgoing constituents, through the use of one FD simulation.
Utilizing this observation,  showed that by FD-injection of recorded multicomponent streamer data into an infinite water layer, without a sea surface, a separation of the wavefield is achieved. By keeping only the upgoing portion of the wavefield, the spatial and temporal bandwidth of the data is increased (i.e. is deghosted on the receiver side). Subsequently, Robertsson et al. (2015) showed that up-and downgoing P-and S-wave constituents can naturally be separated to numerical precision with isotropic elastic FD modelling. Börsing et al. (2017) implemented the method of multiple point sources (MPS;Mittet 1994;Masson et al. 2014), a wavefield injection scheme based on the Kirchhoff-Helmholtz integral, to separate a measured wavefield (recorded close to a planar free surface in a water tank) into its incident and reflected constituents. While the method of MPS is closely related to FD-injection, it differs in the computational implementation.
All of the mentioned surface-based wavefield separation methods fail when the recording surface has sharp corners. On the other hand, FD-injection based wavefield separation methods, applied to reconstruct the wavefield constituents in the interior of the modelling domain, handle sharp corners (and stair-cased surfaces) naturally. Therefore, following Thomsen et al. (2018), we leverage the interior, FD-injection based methods and show how to adapt them to the separation of surface elastic data. We present a novel wavefield separation method requiring only three-component (3C) particle velocity data recorded at the free surface and information about the medium properties. Our scheme is capable of retrieving the first-order incident and full reflected wavefields at a closed surface with sharp corners of 90 • , whilst correctly handling all angles of propagation as well as correctly decomposing surface waves. We show how such a scheme can be used to separate experimental and numerical, 3C free surface data. Our proposed scheme is particularly relevant for isolating the incident and reflected wavefield components in physical laboratory experiments, for example for implementing 3-D elastic immersive wave experimentation (Vasmel 2016). In such experiments, applying the three components of traction incident on the free surface of an elastic object as the components of point-force sources makes it possible to cancel any reflections of the wavefield at said boundary (Thomsen et al. 2019). Furthermore, the incident wavefield can be extrapolated through an arbitrary larger virtual numerical domain before being re-injected into the object, enabling a seamless interaction of physical and virtual wave propagation (Thomsen et al. 2019;Becker et al. 2018). Additionally, our proposed scheme can be adapted to provide a means of separating seismic data acquired at the Earth's surface, as well as for 3-D elastic numerical modelling.
Section 2 introduces FD-injection as a numerical method to reconstruct an elastic wavefield in a subvolume using only recordings surrounding an auxiliary surface enclosing the subvolume. Then, in anticipation of an eventual application to the free surface, we show how FD-injection can be adapted when some of the recordings are exactly on the auxiliary surface, without affecting the reconstruction accuracy. Particularly, we explain how the resulting approach naturally treats the sharp corners and edges of a rectangular injection surface, without introducing any distortions to the reconstructed wavefields. This is followed, in Section 3.1, by a brief summary of how FD-injection has been used previously for wavefield separation along open surfaces, leading into an analysis, in Section 3.2, of the requirements when both the recording and injection occur along a closed boundary. We also explain, why our method is particularly suitable for wavefield separation at a free surface boundary. Finally, after demonstrating the accuracy of our method using synthetic data in Section 4, we show its success separating real, laboratory particle velocity data into their incident and reflected components in Section 5. The laboratory data were acquired on the surface of a granite rock volume using a state-of-the-art robotized 3-D scanning Laser Doppler Vibrometer (LDV).

N U M E R I C A L WAV E F I E L D R E C O N S T RU C T I O N U S I N G F I N I T E -D I F F E R E N C E I N J E C T I O N
FD-injection was originally proposed by Alterman & Karal (1968) as a method to introduce analytical source wavefields in a finite-difference grid in an acoustically transparent manner and to avoid dynamic range issues due to the limited precision of floating-point computations at that time. Robertsson & Chapman (2000) realized that this approach could be adapted to yield an efficient method to recompute full-waveform seismograms for a suite of localized model alterations and coined the term FD-injection. Elastic wave propagation in space is governed by the velocity-stress system (Aki & Richards 2002). In the following, we use a conventional, second-order accurate [O(2,2)], staggered-grid leap-frog discretization to model elastic wave propagation (Virieux 1986;Graves 1996  Those quantities recorded along S are indicated by the light blue rectangle. Whereas, quantities recorded above and below, and subsequently averaged to their midpoint on S, are indicated by the red rectangles. (b) The first boundary FD inj1 runs in between the gridpoints and the arrows indicate the additional injections which must be performed during FD-injection. (c) The second FD-injection boundary (FD inj2 ) is indicated, as well as the associated FD-injection updates. (d) By averaging the updates in (b) and (c), we define a new injection surface S, running through the gridpoints. Any quantities recorded half a grid spacing away from and normal to the defined injection surface are averaged and injected into the associated quantities on the surface S (indicated by the dark, red arrows); simultaneously, recordings of those quantities located directly on the injection surface are spread equally across the associated quantities located half a grid spacing away and normal to the injection boundary (indicated by the light, blue arrows). quantities are given in Appendix A. The full elastic FD system is described by the three component velocity vector (v x , v y , v z ), the normal stresses (s xx , s yy , s zz ) and the shear stresses (s xy , s xz , s yz ). Fig. 1(a) illustrates a XY cross section of the staggered 3-D FD grid (see Fig. A1b).
FD-injection is based on the linearity of the wave equation and compact support of FD stencils (van Manen et al. 2020). It comprises of two steps: first, the wavefield propagating on a numerical model is recorded in a small band of gridpoints normal and directly next to a user-defined injection surface running in between the FD nodes. Figs 1(b) and 1(c) depict two separate rectangular FD-injection surfaces, FD inj1 and FD inj2 (indicated by the dashed black lines), at slightly different locations within the XY cross section (see Fig. 1(a)). Those quantities recorded during wavefield propagation are located a quarter of a grid spacing normal to the defined injection boundary. Subsequently, these recorded wavefield quantities are used in a second wavefield simulation. In this second simulation, the FD updates will be augmented by injections of the recorded quantities. The effect of these injections is, for example, that a wavefield is fully absorbed at the injection surface. It is important to note that the injections take place solely at the gridpoints whose stencil cross the injection surface. The corrections occur in addition to the regular FD updates and are indicated by the arrows in Figs 1(b) and (c). For instance, quantity s xx previously recorded outside of the surface FD inj1 (see Fig. 1b) is used to correct the affected quantity v x , located inside of FD inj1 . Simultaneously, v x , previously recorded inside, is used to correct s xx , s yy and s zz , the affected gridpoint located outside. These corrections create and maintain a discontinuity in the wavefield across the injection surface such that the wavefield is reconstructed on one side of the surface, while no wavefield is introduced on the other side.
Note that in conventional applications of FD-injection, only a single injection surface in between the grid nodes is required. We furthermore note that sharp corners, as depicted in Fig. 1, are handled in a numerically exact way.

Implementation of two averaged FD-injections
The described FD-injection methods are designed for wavefield reconstruction in the interior of the modelling domain. When considering their application at the boundary of the domain, e.g. an exterior free surface setting, we must consider two significant apparent drawbacks. First, the injection surfaces FD inj1 and FD inj2 , depicted in Fig. 1(b) and (c), run in between the nodes, and therefore in between the locations where the physical quantities are defined on the FD grid. Hence, the FD-injection surfaces are not colocated with the chosen boundaries (e.g. free surface) which are typically defined (and applied) at the gridpoints. Secondly, it would require wavefield recordings made at locations offset from the injection surface: for example one would need measurements made above and below the free surface, which is of course impossible in practice.
Here we show that, for a second order scheme, both drawbacks are efficiently circumvented by our proposal of averaging two FDinjections. As shown in Figs 1(b) and (c), we can position the two FD-injection boundaries to surround a single 'layer' of grid nodes. In Fig. 1(d), we illustrate how the combined effect is akin to defining a single boundary that runs through the FD quantities, that is following the black line denoted with S.
As an example, we now derive the necessary updates for the quantities v x (located on S) and v y (located to either side of S), for the top boundary in Fig. 1(a). We assume to have a source inside the domain, and start with zero initial conditions at t 0 . First, we focus on the quantity v x , located below the top boundary FD inj1 (see Fig. 1b) and inside the volume bound by FD inj1 . During the 3-D wavefield simulation, the regular FD update for the velocity node v where the portion of eq. (1) underlined by the dashed line corresponds to the part of the stencil crossing FD inj1 . Therefore, this quantity must be recorded during the initial simulation.
Subsequently, during the wavefield injection simulation, both the regular FD update of v n+1/2 x|i, j,k is computed, as defined by eq. (1), and the previously recorded quantity is used to correct v n+1/2 x|i, j,k as per eq. (2): The left arrow (←) indicates that the FD-injection update is applied in addition to the regular FD update using the recorded quantities indicated by a tilde (−s n xx|i−1/2, j,k ). This additional update, that is the wavefield injection, leads to the reconstruction of the originally present wavefield below FD inj1 [as demonstrated by Robertsson & Chapman (2000)]. Fig. 1(c) shows the location of the second FD-injection boundary, FD inj2 . As before, the parts of the stencil normal to and crossing the injection boundary must be modified during the wavefield injection step. This time, v n+1/2 x|i, j,k is located outside the volume bound by FD inj2 . The portion of eq. (1) representing the FD-injection update for v x , is underlined by a solid line and must be recorded during the initial simulation. Eq. (3) shows the FD-injection update which is achieved by subtracting the normal stress recorded just inside the injection boundary (s n xx|i+1/2, j,k ) from v x located just outside the boundary: Note here that, taken in isolation, the effect of this step is that the state below FD inj2 evaluates to the originally present wavefield.
Our proposal now is to combine the two FD-injection surfaces through averaging. The effective injection into v x may then be written as the combined effect of the two measured stress nodes. We can show this injection term from the combination of the two FD-injection terms as defined by eqs (2) and (3): If the FD updates for both FD-injections are averaged, we are left with one update term for v x , as displayed by eq. (4). The combined injection as per eq. (4) is indicated by the red (dark) arrows in Fig. 1(d). We note an interesting feature of eq. (4), which is that the average of the two s xx may equally be considered a second-order accurate interpolation of s xx towards the surface S.
Besides the modified update eq. for v x , we must also modify the update for v y . When utilizing the first FD-injection boundary FD inj1 , the wavefield above must be made equal to 0, that is only the wavefield propagating within the boundary FD inj1 should be reconstructed. This is achieved for node v n+1/2 y|i−1/2, j+1/2,k , located directly above the transparent surface FD inj1 in Fig. 1(a), by the update term: Similarly, the FD-injection update of v y when defining FD inj2 , depicted in Fig. 1(b), must ensure reconstruction of the wavefield within the volume: Both v n+1/2 y|i−1/2, j+1/2,k , located directly above the transparent surface FD inj1 in Fig. 1(a), and v n+1/2 y|i+1/2, j+1/2,k , located directly below the transparent surface FD inj2 in Fig. 1(b) are updated using the shear stresss n xy|i, j+1/2,k (see eqs 5 and 6). This quantity is located on S when averaging the two FD-injections. Consequently, the recorded shear stress is spread out equally across both v y velocity nodes normal to the injection surface S in the FD-injection step: These two FD-injections occur in addition to the respective regular FD update of v y , as per eq. (A2). The remaining updates to perform during wavefield injection can be derived similarly for all the other FD quantities associated with the defined volume boundary. Appendix B lists all the update terms for the chosen 3-D FD scheme, including the edges and corners of the closed injection boundary. The method now only requires quantities recorded at the boundary S (surrounded by the light blue box in Figs 1a and d) and those quantities that are averaged from points on both sides of the boundary (surrounded by the light red boxes in Figs 1a and d), that is with their effective midpoint on S. Figs 1(d) and A1(d) and (e) show that the injections along the edges and corners of the domain utilize purely recordings made right outside and along the boundary S, that is the updates associated to FD inj1 . All quantities are readily available in numerical simulations.
As required, by averaging two FD-injections a new injection surface S running through the grid nodes is defined. Furthermore, by utilizing FD-injection, results accurate to machine precision can be obtained numerically even in the presence of sharp corners and edges. We show below that for practical applications (i.e. surface wavefield separation), instead of the quantities averaged from points on both sides of the boundary, we can simply use quantities measured on the boundary and still obtain highly accurate results.

F D -I N J E C T I O N B A S E D WAV E F I E L D S E PA R AT I O N
We shortly review how FD-injection can be used along flat open surfaces to achieve wavefield separation. Subsequently, we analyze the modifications necessary for wavefield separation along a closed surface. The FD-injection method allows for a recomputation of wavefields within closed surfaces. In particular, the FD-injection method generates waves from the domain boundaries such that waves can correctly propagate into the domain. Conversely, waves approaching the boundaries vanish due to destructive interference with the FD-injection terms. It is assumed that the internal medium (bound by the injection surface) and the boundary conditions do not change. However, violating this assumption has interesting effects. For example, modifying the velocity model prior to the FD-injection simulation results in waves that propagate beyond the FD-injection boundaries (Robertsson & Chapman 2000).

FD-injection based wavefield separation for open surfaces
The effect of FD-injection can be further modified by recording and injecting along a flat and open surface of infinite length, as well as using a homogeneous velocity model during injection (e.g. Blanch 2012). The open surface may be considered as just one of the sides of a rectangular FD-injection surface along which waves are then generated during the injection process. Using the homogeneous velocity model ensures that no waves are reflected back towards the injection surface. Consequently, both the full up-and downgoing wavefields are reconstructed at and propagate away from the single open surface. Successful applications of this up/down wavefield separation technique have been shown in, for example  and Robertsson et al. (2015).

FD-injection based wavefield separation for closed surfaces
We now consider the state of a reference wavefield recorded along a predefined, closed boundary (denoted by the dashed line in Fig. 2a). The wavefield is created by a source inside the boundary and interacts with a scatterer located above and outside of the boundary. Consequently, the wavefield crosses the closed boundary three times. We apply the wavefield injection scheme as described in Section 2, using the previously recorded wavefield, but disregard the internal source in the modelling. Fig. 2(b) shows that the boundary sources will emit a wavefield of opposite polarity to the originally recorded wavefield along the boundary. Thus in effect, the outgoing wavefield due to the state inside the closed boundary is obtained, albeit with an inverted polarity. All higher-order interactions with the exterior domain are computed. In contrast, the wavefield propagating inside the boundary is cancelled out. However, Fig. 2(c) shows that if we leave out the internal source and make the exterior domain homogeneous, we obtain only the first-order outgoing wavefield with flipped polarity, as well as the complete ingoing wavefield with the correct polarity.
To separate the wavefield at the free surface of an elastic target, we require the injection surface to be co-located with said boundary condition. Hence, we are interested in a case similar to the one depicted in Fig. 2(c). Note, that the scattering interface now coincides with the recording/injection surface as displayed in Fig. 2(d). A source placed inside the domain causes waves to propagate towards and interact with the free surface boundary (events marked 1-3) where the superimposed in-and outgoing waves are recorded. Outgoing waves are defined as those waves incident on the free surface boundary, and ingoing waves as those reflected at the free surface. In the method that we propose, these recordings are then injected along a boundary with the same dimensions as the recording surface during the wavefield reconstruction step. However, rather than having a free surface in the second (i.e. injection) model, we choose an infinite homogeneous domain for the exterior medium (i.e. such that no scattering takes place at the injection surface, as per Fig. 2e). The outgoing waves can then be recorded very close to the injection boundary (e.g. one gridpoint away for a second-order scheme) and extrapolated back to the location of the surface, if required, to find the wavefield incident on the free surface. Note that any outgoing events which have interacted with the free surface more than once (see event 2 in Figs 2d-e) cannot be retrieved as they destructively interfere with the ingoing waves injected on the opposite side of the boundary. This important observation was also made by, for example Thomsen et al. (2018).
Having the free surface as the governing boundary condition provides additional relaxation in terms of which quantities need to be recorded along the surface S. A free surface defines a traction-free boundary condition, that is s ij n j = 0. Therefore, all stresses normal to the Omitting the source and only injecting the recordings made along the boundary, whilst preserving the same model characteristics, recreates only the wavefield outside the transparent boundary. (c) Omitting the source whilst injecting the recordings along the boundary into a homogeneous model leads to a separation of the wavefield along the transparent boundary. (d) A source is fired inside a domain bounded by a free surface, with waves recorded at the free surface. Three separate events are depicted. The first (event 1), propagates towards and reflects off the free surface. It then propagates to the opposite side where it reflects off the free surface again (event 2). A third event is scattered internally before interacting with the free surface boundary (event 3). (e) When injected along a transparent boundary, all events recorded at the free surface are separated into their outgoing (red wave, with opposite polarity) and ingoing (black wave, with correct polarity) parts. Only outgoing events, for which the wavefield interacted with the free surface only once are retrieved. On the other hand, all ingoing events are retrieved.
respective interface of the free surface are zero. This leaves us with needing to record only the three component particle velocity vector at the free surface and injecting the recorded quantities into the associated stresses during wavefield injection (see Appendix B2).

S Y N T H E T I C WAV E F I E L D S E PA R AT I O N F O R A F R E E S U R FA C E W I T H S H A R P C O R N E R S
We now demonstrate our method of wavefield separation using synthetically modelled recordings of the 3C particle velocity vector at the free surface. It is inherently difficult to output accurate free surface data using staggered-grid finite-difference modelling as the modeled wavefield constituents are staggered and do not coincide in space. Therefore, we utilized the spectral element modelling package Salvus (Afanasiev et al. 2018) to generate and record the 3C particle velocity vector at the free surface of a homogeneous cube. The synthetic data generated can be understood as a substitute for a real experiment, where only access to the free surface of the solid target is available. The dimensions of the cube run from 0 to 0.5 m in the x, y and z directions and it is characterized by a density of ρ = 2600 kg m -3 , P-wave velocity of v p = 3200 m s -1 and S-wave velocity of v s = 1200 m s -1 . We excited a vector point source in the x-direction inside the cube, located at (x src , y src , z src ) = (0.3, 0.26, 0.29). The source time function was a Ricker wavelet with central frequency f c = 15 kHz. The particle velocity components v x , v y and v z of the wavefield propagating at the free surface were then recorded at equally spaced gridpoints along the faces, edges and corners of the cube. Every recording node in the Salvus simulation was mapped onto a wavefield injection point when subsequently applying our FD-based wavefield injection scheme in MATLAB. To achieve an accurate FD modelling result for our O(2,2) implementation, the recording/injection nodes were spaced at x < λ/10, with λ = v s /(2 · f c ) = 4 mm being the minimum wavelength of the propagating wavefield. Furthermore, the sampling interval t, adhered to the CFL stability criterion, t < x/(v max √ 3). Consequently, we set x = 4 mm and t = 0.5 μs. Additionally, Salvus was used to create a reference solution of just the wavefield incident on the free surface by simulating a wavefield propagating on a larger domain without a free surface present.
Time gathers of the velocity components v x , v y and v z , recorded along the bottom face of the XY cut plane of the free surface are displayed in Figs 3(a), (c) and (e), respectively. The first wave recorded at the free surface can be identified as the P wave. The second event recorded is the S wave, which can be distinguished from the P wave due to its larger curvature. The events crossing through the gather from either side are reflections of the wavefield from the adjacent faces, edges and corners of the cube. Furthermore, higher-order reflections, that is reflections of the propagating wave with more than one face, are visible at later points in time. We injected these recordings in our finite-difference based wavefield separation scheme. Snapshots of the resulting wave propagation are displayed in Fig. 4. They show that the reconstructed wavefield is not affected by the sharp edges and corners of the model. To retrieve the wavefield incident on the free surface, we record the wavefield propagating to the outside of the injection boundary one gridpoint away. The time gathers in Figs 3(j), (k) and (l) show the retrieved incident wavefield along the bottom face of the XY cut plane displayed in the left-hand column of Fig. 4. The result can be compared to the reference incident wavefield displayed by the time gathers in Figs 3(g), (h) and (i). Figs 3(b), (d) and (f) show traces of the wavefield recorded at the free surface, the retrieved incident wavefield and the reference solution at a single receiver station of the time gathers. For all three components, the retrieved incident wavefield matches the reference solution. As shown in the schematic in Fig. 2(e), we are able to retrieve those events of the incident wavefield that interacted with the free surface boundary once. Taking the divergence and curl of the propagating wavefield during injection makes it possible to further separate the incident P and S waves (Sun et al. 2004). We note that according to our method, recordings of the particle velocities are expected in a staggered grid, that is both on, or directly inside and outside of the boundary (even if this leads to an obvious contradiction for a free surface condition). Even so, in practice all the recordings of the particle velocity were made exactly on the free surface and the wavefield injection scheme still leads to a highly accurate reconstruction of the wavefield as shown by Fig. 3. To achieve this accuracy during the wavefield reconstruction step, knowledge of the model parameters in the interior of the closed boundary of the target was used. Any mismatch of the model parameters defined during wavefield injection to those characterizing the real model cause energy to 'leak' into the retrieved incident wavefield gathers. This is best understood by further analysing Fig. 2(e). If, for example the velocities used during wavefield injection are erroneous, the injected, ingoing (reflected) wave for event 1 (black wiggle line) will not reach the opposite side of the injection boundary at the right time. Consequently, this event will not properly destructively interfere with the outgoing injected wavefield of event 2 (red wiggle line) and energy will 'leak' to the outside of the injection boundary. In Appendix C, we further explain how this observation can be utilized to constrain the unknown model parameters prior to wavefield injection using a minimization scheme. In the test case presented here, however, we required just three elastic parameters to describe the entire homogeneous model during wavefield injection. The left-hand column shows wavefield propagation in the XY cut plane, the middle column propagation in the XZ plane and the right column propagation in the YZ plane. Stepping forward in time, one can observe the incident wavefield propagating away from the injection boundary and the reflected wavefield propagating towards the inside of the boundary. Note that no artifacts are caused by the sharp edges of the cubic injection boundary. See Supporting Information Movie S1 for an animation of the snapshots presented here.

E X P E R I M E N TA L WAV E F I E L D S E PA R AT I O N F O R A F R E E S U R FA C E W I T H S H A R P C O R N E R S
After demonstrating the accuracy of our wavefield separation scheme in a numerical study, we applied our method to real x, y and z-component particle velocity data acquired along the free surface of a solid target. Using a Polytec PSV-500-3-D Scanning Vibrometer (Polytec GmbH 2019), we measured the wavefield propagating along five faces of a granite rock. Fig. 5(a) depicts the three-headed Laser Doppler Vibrometer (LDV) positioned to scan one side of the granite rock. The rock itself is assumed to be homogeneous and is 0.574 m long, high and wide.
A 20-cm-deep and 3-cm-wide hole was drilled into the base of the rock as indicated in Fig. 5(a) into which we placed a multi-axis PICA 143.01 shear actuator (PI Ceramic GmbH 2018) (see Fig. 5b), embedded in a plastic holder (see Fig. 5c), to act as a source. Located at (x src , y src , z src ) = (0.255, 0.207, 0.18) m, the source was amplified using a PiezoDrive PD32 amplifier (PiezoDrive 2019) and used to excite a Ricker wavelet of 20 kHz central frequency along the y-direction. The LDV was then used to measure the wavefield at individual scan points, spaced 2 cm apart, on five faces, 8 edges and the top 4 corners of the rock. The bottom face, that is the XY slice at the origin of the z-axis with the hole for the source, was not scanned. We made 100-ms-long recordings with a 250 kHz sampling rate resulting in a temporal resolution of t LDV = 4 μs. Additionally, the data were bandpass filtered between 2 and 50 kHz. The signal-to-noise ratio of the measurements made by the LDV was increased by placing reflective tape along the surfaces of the rock at each scanning point (see Fig. 5a) and averaging each measurement 200 times. Figs 5(d)-(g) show snapshots of the LDV data recorded. The source located inside the granite rock is activated at t 0 = 20.00 ms after which, at t 1 =20.15 ms one observes the first arriving P wave on the front face of the rock. This event is directly followed by the secondary S wave at t 2 =20.18 ms. Snapshots at times t 3 =20.25 ms and t 4 =20.28 ms show how the wavefield arrives at the top face. Furthermore, the wavefields already propagating along the rocks faces interact with the corners and edges of the rock.
While the real data was recorded using the LDV, the wavefield injection was performed numerically with MATLAB. However, different to the purely synthetic free surface data set, as an additional step before wavefield injection, we had to interpolate the LDV data both in space and time to satisfy the stability criteria for numerical modelling. The LDV recordings were acquired at a spacing of 2 cm and interpolated linearly to a spatial resolution of x =4.1 mm. The LDV data set was then further interpolated in time to t = 0.5 μs.
Figs 6(d), 7(d) and 8(d) show three examples of time gathers of the interpolated v x , v y and v z component LDV data recorded along three separate faces of the rock. We injected the recordings according to our wavefield injection scheme and recorded the separated incident wavefield one gridpoint adjacent to the injection boundary. As in the numerical study, the individual velocity recordings were treated as if they were recorded on a staggered grid. The data were injected into a homogeneous model characterized by a P-wave velocity of v p = 3430 m s -1 , S-wave velocity of v s = 2170 m s -1 and density of ρ = 2600 kg m -3 . A rough estimate for these values was first obtained using transmission experiments and subsequently further optimized using a minimization scheme based on 2-D wavefield injection (see Appendix C).
Snapshots of the resulting wavefield reconstruction are depicted as XY (Figs 6a-c), XZ (Figs 7a-c) and YZ (Figs 8a-c) slices. For each slice, the snapshots depict the v x , v y and v z component. From the individual snapshots, we observe spherical radiation patterns of the wavefield emitted from the injection boundary. Furthermore, at, for example the right boundary in Fig. 6(b), 7(b) and 8(b), no significant artefacts can be observed at the corners (edges of the 3-D injection boundary) during wavefield injection. For each slice, we depict a time gather along one face of the injection boundary, as well as a trace comparison between the recorded free surface wavefield and the retrieved incident wavefield. We focus on comparing the early arrivals, that is the incident P and S waves, of the retrieved incident wavefield gathers to the free surface recordings and observe comparable characteristics. Specifically, the traces of the retrieved incident wavefield are similar in phase and smaller in amplitude, which is expected. As previously stated, any mismatches in the model characteristics, for example velocities between the homogeneous injection domain and physical granite rock in the lab, cause energy to leak to the outside of the injection boundary. This can be observed in later times in the gathers and traces. Note that no wavefield injection currently takes place along the face of the rock with the hole as can be observed in, for example Figs 7(a)-(c) along the open, left boundary. Hence, we can see artefacts appearing during injection along the adjacent faces. To obtain a better understanding of the accuracy of our method when using experimental data, we used Salvus to simulate the wavefield at the free surface using a model matching the physical experiment as closely as possible. This included the presence of a cylindrical hole for the source and the model characteristics defined during wavefield injection of the LDV free surface data. Instead of interpolating the recordings as in the laboratory experiment, we record the synthetic wavefield at all the required positions for the FD-injection simulation.
To simulate the physical actuator, we placed two vector point sources acting in phase in the interior of the rock at the edges of the 3-cmwide hole and excited a Ricker wavelet of 20 kHz central frequency. The synthetic data were then injected according to our wavefield separation (a) (d) (a) (d) We compare the synthetic reference wavefield recorded at the free surface (Figs 6j, 7j and 8j) as well as the associated retrieved incident wavefield (Figs 6k, 7k and 8k) to the free surface LDV time gathers (Figs 6d, 7d and 8d) and retrieved incident wavefield (Figs 6e, 7e and 8e). Overall, the LDV data appears to have a lower bandwidth. Even though the piezoelectric actuator is excited at 20 kHz, factors such as the source mounting or electrical components could cause the observed shift in frequency and introduce ringing to the signal, as could attenuation during propagation through the rock volume towards the free surface. Nevertheless, when analysing time gathers (e.g. Figs 6d and j), comparable radiation patterns can be observed, such as the polarity changes of the wavefield from the left to the right. Furthermore, time gathers in Figs 7(d) and (j) depict a similar sequence of events: an off-centred, S-wave arrival, as well as the dominant events crossing from the left and right caused by the wavefield interacting with the adjacent faces and edges of the rock domain. The effects of not being able to record and subsequently inject along one face of the rock are also visible in the snapshots and time gathers of both Figs 7 and 8. Higher order outgoing events, such as the event crossing from left to right in Fig. 7(d), are not accurately suppressed during the wavefield injection as indicated by the schematic displayed in Fig. 2(e).
Finally, comparing the traces of the retrieved synthetic incident wavefield to the synthetic free surface data (Figs 6l, 7l and 8l) shows a similar reduction in amplitude of the incident wavefield as observed in the real data example (Figs 6f, 7f and 8f). Additionally, the shapes of the traces are comparable between the synthetic and real data examples for the early arrivals, i.e. the first order incident wavefield.
Overall, the synthetic data results support the real data results of our wavefield injection scheme. We conclude that an accurate reconstruction of the incident wavefield has been achieved. Artefacts observed can be attributed to not including the hole in the homogeneous model during injection, and the missing wavefield injection along one face of the model. Hence, the synthetic results allow for a qualitative, yet detailed, comparison with the real data example.

D I S C U S S I O N
The fact that only the first-order outgoing wavefield is obtained during wavefield separation distinguishes this work from other approaches. The retrieval of the wavefield incident at the free surface of an elastic target is an essential component of implementing elastic immersive  Vasmel et al. (2013) first described the concept of immersing a physical laboratory experiment in a virtual, numerical environment. By controlling the physical boundary condition of the experimental target (van Manen et al. 2007;Broggini et al. 2017), it is possible to introduce seamless interactions between wavefields propagating in the physical and virtual domains (Becker et al. 2018). As described in Section 3.2, our method of wavefield reconstruction along a closed boundary will only reconstruct the first order interactions of the incident wavefield with the free surface. Thomsen et al. (2019) numerically show that the traction associated to exactly this constituent of the wavefield must be applied as a boundary condition to effectively cancel any unwanted boundary reflections at the free surface of an elastic target. This is in fact the application for which our method of wavefield separation was developed.
Nonetheless, our scheme relies on the interior of the model to be characterized correctly. In Appendix C, we present a method to constrain the unknown model parameters of the physical target to ensure that the parameters used during wavefield injection are as close to the real model as possible. However, this assumes a homogeneous interior of the model. This is obviously rarely a given and any recordings made at the free surface caused by internal reflections would lead to erroneous event reconstruction if we run the wavefield separation scheme with an erroneous interior model. In particular, the separated ingoing component of the wavefield, propagating into the interior of the injection domain, will not destructively interfere with the correctly reconstructed outgoing waves, as they do not propagate along the same trajectories as in the heterogeneous model they were recorded on. Similar artefacts could occur if the wavefield propagating through the rock volume in the laboratory is subject to attenuation which is not accounted for in the used model during wavefield injection.
The use of internal absorbing boundaries during wavefield injection would however completely alleviate the requirement to have prior information about the possibly, non-homogeneous interior of the model. By placing absorbing boundaries in the interior of the injection domain, the reconstructed, reflected waves, propagating into the interior after injection (event 2. in Fig. 2e), would not reach the opposite side of the injection boundary. The problem of leaking events caused by interior inhomogeneities would thereby be eliminated. Even so, these events do not destructively interfere with outgoing events caused by higher-order interactions with the free surface. Hence, the wavefield injection scheme will reconstruct all waves incident on the free surface of a closed boundary, that is including the incident part of higher-order interactions. Such an internal absorbing layer could perhaps be implemented by an adaptation of the methods utilized in cloaking and holography as presented by van Manen et al. (2007). This is the subject of current research. It is important to note that the elastic parameters at the receiver positions would still be required. However, geophysical methods can be used to estimate these material parameters, (e.g. via travel time or surface wave tomography and the use of spatial gradient data Socco et al. 2010;Robertsson & Curtis 2002). Furthermore, initial testing has shown that the medium parameters can also vary laterally along the free surface.
The multicomponent particle velocity recordings utilized by the proposed method must be sampled to at least the Nyquist criterion in both space and time to allow for an accurate interpolation to the FD grid. We note that interpolation methods, such as sinc-interpolation, could reduce possible interpolation errors. However, in the real data example presented, the free surface data of f c = 20 kHz dominant frequency was acquired at a spatial sampling of 2 cm using the LDV. Thus, the dominant wavelength λ dom = v s /f c = 2170 m s /20 kHz = 10.85 cm is significantly oversampled with respect to Nyquist and therefore the linear interpolation method used does not pose a problem.
We note that our wavefield separation method can be implemented for any closed boundary, as well as for open surfaces. As mentioned in Section 3.1, utilizing FD-injection for wavefield reconstruction along a flat, open injection boundary results in a reconstruction of the full incident and reflected wavefields. Both our numerical and experimental results show that we require only recordings made at the free surface of the target. This novelty expands on the findings by Robertsson & Chapman (2000), who required more than one recording datum. Therefore, our method is also of relevance for wavefield separation of, for example 3C land seismic data. In addition, the method presented could be used to more accurately treat corners commonly emerging when approximating irregular shaped free surface topographies in FD schemes. Such use cases are subject of future research.
Finally, we remark that our implementation of averaging two FD based wavefield injections can be interpreted as as a numerical representation of multiple point-source (MPS) injection along a closed boundary. Numerical implementations of MPS based wavefield injection can be found in for example Mittet (1994), Vasmel & Robertsson (2016) and Aaker et al. (2020). However, MPS does not yield machine precision accuracy at corner points in elastic FD simulations, whereas our implementation can handle sharp corners in elastic media up to machine precision, at least for the second-order accurate FD scheme used here. Hence, our implementation is also applicable for local, 3-D, elastic numerical modelling problems.

C O N C L U S I O N
We have presented a method to retrieve the incident and reflected wavefields at the free surface of an elastic target by injecting measured three-component surface particle velocity data in a modified FD wavefield injection scheme. For an open (e.g. horizontal) surface, the method only requires knowledge of the model characteristics directly at the receiver positions, while for a closed surface the properties in the interior of the target bound by the free surface are also required. In the latter case, only the incident part of the first-order interactions with the free surface are recovered while also correctly separating the wavefield in the sharp corners and edges of the injection surface without visible distortions. Furthermore, the method accounts for all angles of propagation. In an example, we successfully separated a wavefield experimentally recorded along five sides of a cubic granite rock volume. The ability to isolate the wavefield incident on the free surface of an elastic target can be utilized in the development of experimental focusing and imaging techniques. Specifically, retrieval of the incident part of the first-order interactions of the wavefield with the free surface is an important requirement for realizing 3-D elastic immersive wave experimentation.

A C K N O W L E D G E M E N T S
The authors would like to thank two anonymous reviewers for their valuable feedback. Furthermore, the first author thanks Nele Börsing, Theodor Becker, Patrick Elison and Xun Li for insightful discussions and suggestions. We also thank the Mondaic team for their support with Salvus. This research was supported by the Horizon 2020 ERC grant 694407 -MATRIX.

DATA AVA I L A B I L I T Y
Data associated with this research are available and can be obtained by contacting the corresponding author.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Movie S1. The animation displays wavefield propagation during wavefield injection along the black square. Displayed are the parameters v x , v y and v z (first to third rows) for several cut planes though the injection cube (X Y , X Z and Y Z from left to right). Starting from 0.10 to 0.50 ms, one can observe the incident wavefield propagating away from the injection boundary and the reflected wavefield propagating towards the inside of the boundary. Movie S2. A 3-D animation of the data acquired with the LDV on the granite rock. The time trace at the bottom displays the three components of the measured velocity at the indicated scan point. Visible are the arrival of the primary and secondary wave, followed by their respective reflections off the edges and corners of the rock. Movie S3. Animations of the propagating wavefield during wavefield injection at the black square using the data acquired with the LDV along five faces of the granite rock. Displayed are the parameters v x , v y and v z (first to third rows) for several cut planes though the injection cube (X Y , X Z and Y Z from left to right). One can observe the incident wavefield propagating away from the injection boundary and the reflected wavefield propagating towards the inside of the boundary. Clearly visible is the effect of the missing wavefield injection along one face of the model where no data were acquired with the LDV. Movie S4. Animations of the propagating wavefield during wavefield injection at the black square using synthetically modelled data. Displayed are the parameters v x , v y and v z (first to third rows) for several cut planes though the injection cube (X Y , X Z and Y Z from left to right). One can observe the incident wavefield propagating away from the injection boundary and the reflected wavefield propagating towards the inside of the boundary. Compared to Movie S3, the animations show exactly the same sequence of events propagating away from and into the square injection boundary. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

1661
The individual particle velocity and stress constituents must be recorded and injected along the defined boundary S in the following way: (i) For the six faces of the defined injection boundary, any quantities recorded just inside or outside of the defined boundary S are averaged and injected into the quantity located on the boundary S; simultaneously, the recordings of the quantities located on the faces of the injection boundary S, are spread evenly across the associated gridpoints located just inside or outside of the boundaries faces (see Fig. A1c).
(ii) Quantities located on the 12 edges or 6 corners of the boundary are recorded and spread out equally to the associated nodes purely outside of the injection boundary (see Figs A1d and e). Consequently, only those quantities recorded outside of the edges and corners are injected into the associated nodes on the edges or corners (see Figs A1d-e).
We define all the injection updates required for the XY-plane as a blueprint for the implementation of the other two cut planes. The sign convention during the wavefield injection is determined by the normal vector n = [n 1 , n 2 , n 3 ] at the respective face, edge or corner of the injection boundary.   Figs 2(d) and (e) displays how our wavefield separation scheme can be used to retrieve the wavefield incident on the free surface. Due to the closed injection boundary, we retrieve exactly the incident wavefield which interacted with the free surface once (events 1 and 3). Any higher order interactions with the free surface are not part of the retrieved incident wavefield as they destructively interfere with the ingoing waves injected on the opposite sides of the injection boundary. This observation holds only if the characteristics of the interior of the injection volume are known. In the simplest case the interior of the target is homogeneous. An erroneous estimation of the physical characteristics of the domain, such as the P-and S-wave velocities, leads to energy leaking into the outgoing wavefield (event 2) as the timing between the ingoing wavefield at event 1 propagating through the domain does not perfectly interfere with the injected outgoing wave of opposite polarity at event 2. Consequently, it is possible to vary the physical parameters of the model during wavefield injection, aiming to minimize the total energy of the recorded incident wavefield gathers. We find that the energy is minimized when we only retrieve the first order interactions and no energy leaks outside of the injection domain. The same holds if the interior of the model is inhomogeneous. In this case only the correct interior model leads to a minimization of the residual energy.
We used this observation to optimize the unknown model parameters v p and v s prior to wavefield injection of the LDV free surface data in Section 5. First, the data were converted from 3-D to 2-D along four faces of the rock using amplitude pre-processing (Wapenaar et al. 1992) to facilitate a much faster retrieval of the incident wavefield. The resulting 2-D data were then injected along a closed 2-D injection boundary (XY cross section) according to our O(2,2) FD-injection scheme. Next, we retrieved the incident wavefield for several v p , v s pairs. By minimizing the total energy of the resulting time gathers, we found the optimal velocity pair which was then used in the full scale, 3-D wavefield injection.