Internal absorbing boundary conditions for closed-aperture wavefield decomposition in solid media with unknown interiors

Physical modeling in a laboratory often involves the use of small-scale (analog) models of natural or synthetic materials and constitutes a primary tool for the study of mechanical wave propagation in a variety of settings. Until the late 1970s, physical modeling was used as a proxy for seismic field experiments to study how recorded seismograms can be interpreted for subsurface geological structures (Angona, 1960; Brien and Symes, 1971; Hilterman, 1970; Oliver et al., 1954). In the 1980s and 1990s, the advent of full wavefield numerical modeling enabled simulations [mainly two-dimensional (2D) acoustic and elastic] to be run on the same scales as seismic field experiments (Levander, 1990; Marfurt, 1984; Robertsson et al., 1994; Virieux, 1986), and as a general trend, physical modeling fell out of favor despite the fact that numerical modeling requires strong assumptions about the physics of wave propagation that may not hold in the natural world (Mo et al., 2015).

Since the 2000s there has been a resurgence of physical modeling research and development due to a wide range of interest in, amongst others, material anisotropy (Qi et al., 2015; Tsvankin, 2012; Yang et al., 2016), fracture characterization (King, 2002; Nosjean et al., 2020; Pinto et al., 2018), noise interferometry (Curtis et al., 2006; Snieder et al., 2002; Wapenaar et al., 2010), and wave-based imaging, tomography, and inversion applications (Brenders and Pratt, 2007; Duan et al., 2017; Fichtner, 2010). These wave phenomena and associated applications can be challenging to study using numerical modeling. For example, attenuation and dispersion in natural materials are frequency-dependent wave phenomena (Jakobsen and Chapman, 2009), but using a fully linear (numerical) model commonly involves the assumption that the medium has an equal response at different frequencies. Such simplifications of real-world wave interactions with structures encountered in materials in the Earth's subsurface and other environments can cause problems when trying to replicate mimetically the data observed in the field or laboratory. As a result, numerous new laboratories have been developed for physical experimentation over the past two decades (Arthur et al., 2012; Bretaudeau et al., 2011; Mikesell and van Wijk, 2011; Sivaji et al., 2002). A new wave laboratory at ETH Zürich adopts a radically different approach aimed at overcoming such issues by taking full control of the physical boundary conditions imposed on wavefields, paving the way for complex wave experimentation in real three-dimensional (3D) volumes (Becker, 2019; Börsing, 2019).

The internal structure of 3D solid experimental (rock) volumes can of course be probed using elastic waves, but the limited accessibility of solids usually only allows wavefield measurements to be made at their traction-free surface.¹ Processing and interpreting such recordings is challenging because the data necessarily involve a superposition of (1) outgoing waves (incident on the free surface from the interior) and (2) their immediate reflection and mode-conversion at the free surface (i.e., ingoing waves), as shown in Fig. 1(a). The same challenge exists when processing seismic data acquired at the Earth's surface, in which context the solution usually involves decomposing the data into upgoing (incident) and downgoing (reflected) constituents. Best practice in seismic data processing suggests that decomposing free-surface recordings is also an essential step in the laboratory, and methods for doing so are discussed in detail below. The free surface also gives rise to so-called surface-related multiples (copies of the primary scattered wavefield) as their surface reflections interact with the interior for a second, third, fourth, … time. Eliminating such free-surface-related multiples and replacing them with more meaningful interactions from a virtual exterior that better represents the context of the sample embedded in natural or artificial worlds is also a key driver for the ETH laboratory, as discussed below.

Elastic wavefields can be decomposed using numerical wavefield injection in time-domain finite-difference (FD) modeling (Amundsen and Robertsson, 2014; Ravasi and Curtis, 2013; Robertsson et al., 2015; Vasconcelos, 2013). The experimental data are first recorded along an array of densely spaced receivers, and the data are subsequently injected into a FD simulation as the signatures of a set of densely distributed multi-component point sources. These sources are arranged in the same geometry as the receiver array deployed in the physical experiment and embedded in an extended numerical model with the same material properties as the true medium inside the wavefield injection surface. Figure 1 illustrates this physical data acquisition and the associated wavefield injection in two dimensions [Figs. 1(a) and 1(b)]. The FD propagator, running alongside the wavefield injection, radiates the out- and in-going parts of the data [rays 1 and 2 in Fig. 1(a)] away from the injection surface $Ssep$ in their correct propagation directions [rays 1 and 2 in Fig. 1(b)], thereby realizing the separation (Robertsson et al., 2015).

This methodology works very well for data that are acquired on a single side of the experimental domain (Vasmel and Robertsson, 2016). However, for a closed-aperture recording geometry, it was found that this injection-based wavefield decomposition does not isolate all outgoing energy present in the recorded data (Thomsen et al., 2021): only primary outgoing waves (wave energy that did not previously interact with the free surface) are isolated. Furthermore, a practical limitation exists with this approach: the elastic model used in the FD simulation has to match the (unknown) physical experimental domain inside the wavefield injection surface $Ssep$ exactly (Thomsen et al., 2021). As explained below, these limitations can be overcome using a straightforward strategy: prevent the synthetically reconstructed ingoing waves from propagating across the interior such that they do not reach the other sides of the injection surface.

We propose an internal absorbing boundary condition (IABC) that cancels the reconstructed ingoing waves [ray path 2 in Fig. 1(b)] during the FD simulation, affecting the 3D closed-aperture wavefield separation such that all orders of the outgoing part of free-surface interactions are retrieved correctly from the recorded data. This IABC is developed using a combination of active and passive boundaries. The active boundary cancels the undesired ingoing waves using a prediction-annihilation scheme, closely related to immersive boundary conditions (Broggini et al., 2017; Li et al., 2022; van Manen et al., 2015). As the active boundary is transparent for propagating (and evanescent) waves, it can be complemented by a passive boundary inside, which is used for wavefield extrapolation and within which we rely on damping profiles made from convolutional perfectly matched layers (C-PMLs) (see Martin and Komatitsch, 2009; Roden and Gedney, 2000). To our knowledge, the resulting injection-absorption-based separation scheme is the only scheme that enables all orders of outgoing waves to be separated for closed-aperture configurations. As discussed below, using passive absorbing boundaries only cannot effectively absorb ingoing wave energy due to the limited number of FD grid points around the edges of the absorbing region. The IABC is the key enabler of so-called elastic immersive wave experimentation currently under development at our ETH laboratory, by which a physical wave experiment can be immersed into a virtual environment (Thomsen et al., 2019). The IABC allows the physical experiment to involve arbitrary medium with unknown physics in the interior, which is desired in general laboratory wave applications.

In Sec. II, we introduce the elastodynamic governing equations and recapitulate the theory for numerical wavefield injection of free-surface data. Next, we present the theory for the new elastic IABC and discuss the experimental setup in detail. In Sec. III, we present wavefield separation results for 3D synthetic free-surface data with the internal absorbing boundary activated. Since the results are complicated, even for the simplest case of a fully homogeneous 3D elastic medium, we employ compressional/shear (P/S) separation and ray tracing to explain many of the features. We then present wavefield separation results for the recorded experimental free-surface data, obtained using the same methodology. In Sec. IV, we discuss the difficulty of building such an IABC and the application of the IABC in elastic immersive wave experimentation before summarizing the conclusions.

Mechanical wave energy propagating through a solid medium carries information about material properties within. For an elastic experimental domain, wave propagation can be described by two coupled first-order differential equations relating the tensorial stress wavefield $τij(x,t)$ and particle velocity wavefield $vi(x,t)$⁠. The equation of motion follows (Pujol, 2003):

where x represents an arbitrary location in a Cartesian coordinate system, t is time, $ρ(x)$ is mass density, and $fi(x,t)$ denotes a distribution of body force density sources. Einstein's summation convention applies to repeated subscripts (here, for i, j representing spatial axes x, y, z). The additional constitutive relation between the stress and particle velocity wavefields is

where c_ijkl is a fourth-rank stiffness tensor containing the elastic parameters of the experimental domain, and $hkl(x,t)$ denotes a distribution of deformation rate density sources.

In a physical modeling experiment [i.e., Fig. 1(a)], particle velocities $vi(xF,t)$ are recorded at the free surface $Sfree$ (⁠$xF∈Sfree$⁠), which is also the boundary of the experimental domain. The recorded data $vi(xF,t)$ can be injected into a numerical simulation carried out in a time-space (⁠$t−x$⁠) domain in which Eqs. (1) and (2) are discretized in space and solved alternately every half time step (Virieux, 1986). The injection relies on using the following source term in the simulation (which is a function of the recorded free-surface data):

where the injection surface $Ssep$ [Fig. 1(b)] has the same geometry as the free surface $Sfree$ enclosing the experimental domain [Fig. 1(a)], $xS$ denotes a wavefield injection point on $Ssep$⁠, h_ij denotes the deformation rate density source [Eq. (2)], and $δ(·)$ denotes the Dirac distribution. Note that no distribution of body force sources is involved in this source injection, since their role is to inject signatures of tractions at the free surface that are known to be zero. The injection surface $Ssep$ is a transparent surface in the numerical simulation, but its counterpart $Sfree$ is a reflecting boundary in the physical experiment. Wavefield injection requires the correct model parameters (i.e., $cijkl,ρ$⁠) to be set at the injection locations (⁠$xF$ on $Ssep$⁠), matching the physical parameters at the recording locations (⁠$xS$ on $Sfree$⁠) in the physical experiment. This wavefield injection radiates the outgoing (incident) constituents of the recorded data outward and their ingoing (reflected) constituents inward away from the injection surface $Ssep$ (the former with opposite polarity), resulting in the incident wavefield being retrieved outside the surface $Ssep$ with opposite polarity. However, wavefield separation by injection along a closed surface is limited by an important observation: only the primary incident wavefield can be reconstructed and propagated away from the wavefield injection surface $Ssep$ (Thomsen et al., 2021). This can be understood from the fact that the separated ingoing waves that propagate across the interior of the numerical domain destructively interfere with the separated outgoing waves, injected with opposite polarity, along the other sides of the numerical wavefield injection surface $Ssep$⁠, thereby preventing the reconstruction of the higher-order outgoing waves at the surface $Ssep$ [e.g., ray path 3 in Fig. 1(b)].

The domain enclosed by the injection surface $Ssep$ in the numerical simulation has exactly the same size as the experimental domain enclosed by the free surface $Sfree$⁠. In the numerical domain, the part of the model outside of the injection surface must be kept homogeneous with elastic parameters equal to those of the physical background medium to avoid back-scattering of the separated outgoing energy. The physical background medium is the embedding medium surrounding the heterogeneities in the experimental volume (Fokkema and van den Berg, 2013). When the elastic parameters near the free surface are constant, the background parameters are set to these constant near-surface parameters. When the near-surface elastic parameters vary in space, the background parameters are set to their average values, and in this case, one needs to interpolate the model around the injection locations such that the model is smooth from the laterally varying near surface to a constant interior and exterior. This helps to avoid unwanted scattering at the injection surface.

In this paper, we employ a FD simulation that is second-order accurate in time and space (Virieux, 1986). This FD method is implemented on staggered velocity and stress grids, and the wavefield injection along surface $Ssep$ is implemented with a method of multiple-point sources (MPS) (Li et al., 2022; Thomsen et al., 2021). In the FD simulation, the grid spacing is chosen to be 25 points (⁠$25·Δx$⁠) per dominant S wavelength (>20 is a rule of thumb), and the time step $Δt$ satisfies the Courant–Friedrichs–Lewy (CFL) criterion (Igel, 2017). In the numerical simulation, a radiation boundary condition is implemented around the simulation domain using C-PMLs to absorb the injected, separated outgoing waves (Komatitsch and Martin, 2007).

As explained, the wavefield injection method alone cannot be used to isolate the full outgoing (incident) wavefield from the recorded free-surface data because the reconstructed ingoing waves propagate across the interior of the FD simulation and destructively interfere with the higher-order outgoing waves on the other sides. In principle, a straightforward solution is to prevent the propagation of the ingoing waves across the interior after they are produced at the injection surface $Ssep$ by using an internal absorbing boundary.

A first attempt to implement such an absorbing boundary could be based on so-called immersive boundary conditions (IBCs) by which (numerical) wavefields can be actively canceled at the edge or in the interior of a computational domain by predicting the fields arriving there ahead of time (van Manen et al., 2015). Figure 2 shows the geometry of an internal absorbing boundary that is deployed inside of a wavefield injection surface $Ssep$⁠. To cancel the undesired ingoing waves, we place an array of numerical sources on the emitting surface $Semt$⁠, whose signatures are calculated by extrapolating observed wavefield values from points in between the surfaces $Semt$ and $Ssep$⁠. These points form the (passive) recording surface $Srec$⁠, as shown in Fig. 2, and allow us to observe the reconstructed ingoing waves and extrapolate them to the emitting surface $Semt$ at each time step of the FD simulation.

However, as in IBCs, the wavefield extrapolation relies on computing the elastic Kirchhoff–Helmholtz integrals, which involve connecting all FD grid points on the recording surface $Srec$ to all FD grid points on the emitting surface $Semt$ (which is the surface that actually cancels the wavefield). Such an extrapolation approach not only cancels the ingoing wavefield impinging on the emitting surface from a particular side but also reproduces the same time-evolved wavefield on the other side (as if the wavefield had actually transmitted through the region enclosed by $Semt$⁠). Hence, this naive approach of using the all-to-all elastic Kirchhoff–Helmholtz integrals does not achieve any (internal) absorbing effect for ingoing waves.

To be able to realize an absorbing boundary condition, we still extrapolate the wavefield measured at the recording surface $Srec$ (see Fig. 2) to obtain the wavefield quantities needed for building the active boundary condition at the emitting surface $Semt$⁠. However, this extrapolation is slightly different from that in IBC theory (Li et al., 2022) as we only extrapolate wavefields from the recording surface $Srec$ to the nearest face of the emitting surface $Semt$⁠. To this end, we first introduce a scaling factor α depending on the locations of the sources on the emitting surface $Semt$ and the receivers on the recording surface $Srec$⁠,

where the operator $·$ represents the dot product of two vectors: g (the normal of $Semt$⁠) and $xemt−xrec$ (a vector pointing from a recording point on the recording surface $Srec$ to a source point on the emitting surface $Semt$⁠). Incorporating this scaling factor α into the extrapolation integrals allows only the ingoing wavefield on the same side to be extrapolated,² while the transmitted waves that go through the region enclosed by $Semt$ are not extrapolated to the other sides of $Semt$⁠, which gives

and

where the time variable t is omitted, m_l is the normal vector component of the recording surface $Srec$⁠, the Green's functions (GFs) Γ are pre-computed prior to performing the FD simulation (for wavefield separation) and are associated with a homogeneous free-space³ model whose elastic parameters are equal to those of the numerical medium between the emitting and recording surfaces when performing the separation. Without inserting this factor α in the wavefield extrapolation integrals, or when setting α to a constant value of one, Eqs. (5) and (6) will act exactly as the extrapolation integrals used in IBCs. Figure 3 illustrates the principle of defining the factor α at the upper left corner of the recording and emitting surfaces (in 2D), and the situation for corners, edges, and faces in 3D follows the same principle.

The GFs used in the wavefield extrapolation integrals [Eqs. (5) and (6)] are computed in a homogeneous medium using the same time-domain FD method as used for wavefield separation. To apply them in conjunction with the FD simulation, Eqs. (5) and (6) need to be discretized in space and time [refer to Li et al. (2022) for details]. In a 3D simulation, computing the extrapolation integrals at each time step (of the simulation) is computationally intensive and highly limited by the memory resources available for storing the pre-computed time-domain GFs, which connect a massive number of grid points across the recording surface $Srec$ (i.e., one receiver at every FD grid point) to a massive number of grid points across the emitting surface $Semt$ (i.e., one source at every FD grid point). However, the definition of the factor α in Eq. (4) means that ∼5/6 (for a cube in 3D) of these GFs do not need to be computed or stored for evaluating Eqs. (5) and (6) during the FD simulation. Also, since the GFs are associated with a homogeneous model, translational invariance can be exploited to further compress them, saving both computational cost and memory, as described in detail in  Appendix A.

This active, absorbing boundary provides an effective solution to the problem caused by the ingoing waves in the wavefield separation for almost all of the propagating waves. After terminating the propagation of these waves in the FD simulation, all the incident (outgoing) parts of the second- and higher-order interactions with the free surface will no longer be canceled by destructive interference and will be radiated away from the wavefield injection surface $Ssep$⁠. Hence, the outgoing part of all orders of free-surface interactions can be retrieved outside the surface $Ssep$⁠. Also, accurate estimation of the physical model parameters in the deeper interior of the experimental domain is no longer required for the wavefield separation because the ingoing waves can no longer propagate inside the closed injection surface $Ssep$⁠, except in the small region between $Ssep$ and $Semt$⁠, which are typically chosen to be only two (full) FD grid points apart (⁠$2·Δx$⁠). Note that the elastic parameters at the surface of the experimental volume are still needed because the numerical elastic parameters at the injection surface $Ssep$ in the FD simulation must be identical to physical elastic parameters there.

In practice, we found that the scheme shown in Fig. 3 for defining the factor α is problematic for wave energy that travels toward the sharp corners and edges (in 3D) of the emitting surface $Semt$ because the wave path associated with a corner or edge (e.g., the black arrow in Fig. 3) cannot be clearly judged to be impinging (α = 1) or transmissive (α = 0). In practice, we set α = 1 for any wave energy propagating toward the corners and edges of the emitting surface $Semt$⁠, which results in minor artefacts in the FD simulation (we also tried α = 0 for this case, which resulted in more artefacts). These artefacts consist of diffracted waves originating at the corners and edges of the emitting surface $Semt$⁠, due to (1) the sharp transition between canceling impinging waves (α = 1) and not producing transmitted waves (α = 0) on the emitting surface $Semt$ and (2) an abrupt end to each side of the recording surface $Srec$ (namely the limited aperture of the recorded wavefields for wavefield extrapolation).

To mitigate such artefacts, we add C-PML layers inside the emitting surface $Semt$ such that the artefacts can be further attenuated when they propagate inside. These layers consist of the damping profiles (e.g., d_x and d_y shown in Fig. 2) made from C-PMLs (Komatitsch and Martin, 2007). Compared to the classical C-PML profiles used to attenuate outgoing waves and providing the effect of external absorbing boundary conditions, the C-PML profiles incorporated inside of the emitting surface $Semt$ are spatially reflected through each normal direction of $Semt$ such that the incorporated absorbing layers can attenuate ingoing waves. The modified C-PML profiles are shown in  Appendix B. Figure 2 shows the geometry of the whole IABC, which includes the active boundary on the surface $Semt$ and the passive boundaries inside $Semt$⁠, both of which are used jointly for the wavefield separation of the synthetic and laboratory free-surface data. In Sec. II D, we briefly introduce the setup for 3D elastic wave experimentation, before presenting the synthetic and real data results.

The experimental setup used to study 3D elastic wave propagation in solids is shown in Fig. 4(a) (Masfara et al., 2020; Thomsen et al., 2021, 2019). A scanning laser Doppler vibrometer (LDV) is used to acquire the three-component particle velocity vector across all accessible surfaces of a target object (Blum et al., 2010; Pouet et al., 2012; Zaccherini et al., 2020), in this case a cubic granitic block. The LDV system provides high-fidelity wavefield measurements with a broadband frequency response; in particular, the measurement is non-contact, avoiding receiver coupling issues (Scruby and Drain, 1990). A single, highly repeatable piezoelectric source is placed inside the rock and excites elastic wave energy in the experiments. To deploy this source, a vertical borehole was drilled in the bottom face of the rock and subsequently refilled using segments of the core and epoxy resin (Sikadur-42 HE). For the LDV measurements, patches of reflective tape of 1 cm in diameter, as shown in Fig. 4(a), upper inset, are used on the surface to enhance the reflection of laser light, which increases the signal-to-noise ratio (Gasparetti and Revel, 1999; Hasanian and Lissenden, 2017). The patches are evenly spaced 2 cm apart.

The excellent repeatability of the piezoelectric source allows the physical experiment to be repeated many times, setting up the same propagating wavefield in the rock over a long period of time (e.g., days). This allows the LDV system to acquire the particle velocities at the densely sampled set of points with high data quality: for each measurement point, the same physical experiment is repeated 400 times, and the recorded signals are averaged to improve the signal-to-noise ratio. The time interval $tLDV$ between two consecutive experiments is set to $tLDV=60 ms$⁠, which is long enough for the wave energy to be dissipated to the background noise level at the end of each individual measurement. The tri-axial LDV head can be programmatically moved by the robotic arm for multi-position acquisitions such that all the measurement points on the rock's surface $Sfree$ can be reached. The fully automated data acquisition of 5020 scan points (i.e., 15 060 seismograms) takes around 64 h. Figures 4(b)–4(d) show the recorded wavefield at the surface of the rock, and Fig. 4(e) shows one example of the (averaged) recorded particle velocities at a single measurement point on the rock.

In Table I, we summarize the configuration of the experimental domain, the medium properties of the granite rock, and the parameters of the LDV-based wavefield measurement in the laboratory. In the wave propagation experiments, the piezoelectric source placed inside the rock volume exerts an upward-pointing force (along the z direction) with the source signature corresponding to a Ricker wavelet with peak frequency f_p = 20 kHz. A bandpass filter between 2 and 45 kHz was applied to each trace of the recorded seismic signals.

We assume that the granite rock is an isotropic, elastic medium (which is a good approximation everywhere except at the refilled hole). In this case, the stiffness tensor c_ijkl in Eq. (2) can be replaced by the Lamé parameters λ and μ via $cijkl(x)=λ(x)δijδkl+μ(x)(δikδjl+δjkδil)$⁠, where δ is Kronecker delta. The Lamé parameters λ and μ are related to the compressional and shear wave velocities V_P and V_S (in Table I) via $VP(x)=[λ(x)+2μ(x)]/ρ(x)$ and $VS(x)=μ(x)/ρ(x)$⁠. Note that the interior numerical model used for wavefield separation in practice also does not include any inhomogeneities because the exact shape and properties of such inhomogeneities are generally unknown. Instead, the wave velocity and density models used in the wavefield separation are all kept homogeneous with constant values equal to those of the granite medium.⁴

We first show the wavefield injection of the free-surface data obtained from a synthetic wave propagation experiment, without using the IABC. This synthetic experiment is simulated with the spectral element modeling (SEM) software SALVUS (Afanasiev et al., 2019) and resembles the physical experiment carried out in the laboratory (Fig. 4). The model used in the synthetic experiment is set to be homogeneous with the same dimensions as the rock volume, and the elastic parameters of the model are set equal to those of the granite rock (Table I). In the synthetic experiment, a body force is excited along the z direction (i.e., an f_z source) at the same location as the piezoelectric source in the physical experiment, and the source signature corresponds to a Ricker wavelet with peak frequency f_p = 20 kHz. Figure 5 shows the resulting elastic wavefield in the synthetic experiment. The f_z source causes both P and S waves that propagate with different velocities; these waves are reflected at the free surface of the simulation domain, generating (secondary) reflected and mode-converted waves.

In this 3D numerical experiment, we record the modeled three-component particle velocity vector seismograms at a dense set of points across the six planar faces of the free surface that also bounds the simulation domain. For wavefield decomposition, the recorded data are injected into an enlarged 3D numerical domain in which a time-domain FD propagator is running simultaneously with the wavefield injection as described in detail in Sec. II B. The model used in this FD simulation is fully homogeneous, with the elastic parameters exactly matching those used in the synthetic experiment. The parameters of the FD model are given in Table I.

Figure 6 shows the wavefield injection of the synthetic free-surface data. The corresponding movie is Mm. 1. The wavefield injection on transparent surface $Ssep$ radiates the incident wavefield constituents of the recorded data outward from $Ssep$⁠, while the reflected, ingoing wavefield constituents are reconstructed inside of $Ssep$⁠. As discussed in Sec. II B, the reconstructed outgoing waves that propagate outside of the injection surface $Ssep$ consist of only the first-order constituents of the recorded free-surface data, including only the primary P and S components (first P and first S) (with the former being separated earlier and propagating faster than the latter).

As explained previously, the wavefield injection on the closed aperture $Ssep$ cannot be used to isolate all the outgoing constituents from the recorded free-surface data. Instead, Fig. 6 only demonstrates how to retrieve the first-order incident wavefield constituents from the data. In  Appendix C, we show how this first-order wavefield retrieval relies on using the true elastic model, which has to exactly match the physical experimental domain inside the injection surface $Ssep$ in the FD simulation. However, in almost all cases of interest in laboratory experimentation, this is not realizable.

We now demonstrate the IABC methodology using the synthetic free-surface data generated on the homogeneous interior domain. Figure 7 shows the wavefield separation of the homogeneous free-surface data when using the IABC inside of the injection surface $Ssep$⁠. The corresponding movie is Mm. 2. In this case, the ingoing waves reconstructed on the injection surface $Ssep$ are almost completely annihilated as they cross the emitting surface $Semt$ propagating inward, and as expected, the volume enclosed by $Semt$ shows almost no wave energy compared to the exterior. Any remaining ingoing wave energy is further attenuated by the C-PMLs inside. Simultaneously with the IABC absorbing the unwanted ingoing waves, all orders of outgoing waves are reconstructed and successfully radiated away from the injection surface $Ssep$ (compare to Fig. 6, where the IABC is not used).

Figures 8 and 9 show the separated P and S wavefield constituents by taking the divergence-free and curl-free parts of the data (Robertsson et al., 2015; Yan and Sava, 2008). Superimposed, we show wave fronts computed from ray paths of the primary and secondary wave energy that propagates in the corresponding homogeneous synthetic experiment (Fig. 5, i.e., with the free surface present), accounting for the exact location of the interior source with respect to the six free surfaces in the 3D volume (Virieux and Farra, 1991). Note that these wave fronts are calculated only from the geometry of the experiment and thus only represent the kinematics of the wavefields (e.g., arrival times) without containing information about the wave dynamics (e.g., amplitudes).

Overall, there is a very good agreement between the wavefield separation results and the ray tracing. However, some marked wave fronts in Fig. 8 do not appear to correspond to matching separated arrivals outside the injection surface $Ssep$⁠; for example, on the x-y plane (first column), the separated primary P wavefield is barely visible. This is because in the synthetic experiment, a directional body force source in the z direction is used ( f_z) [resembling the piezoelectric source in the laboratory experiment], and this is different from the omnidirectional source assumed when calculating ray paths using the ray-tracing technique. Also, the wavefield quantity shown (vertical particle velocity v_z) is directional, giving rise to the characteristic radiation patterns of the separated wavefield (e.g., the separated primary S wave) as observed in Fig. 9.

Next, the wavefield separation methods demonstrated in Sec. III A are applied to the free-surface data acquired in the laboratory [shown in Figs. 4(b)–4(d)]. For densely spaced wavefield injection along the surface $Ssep$ in the FD simulation, the experimental data acquired at a spatial sampling of 2 cm (i.e., the interval between the reflective patches) are interpolated to the FD grid sampling $Δx=0.50$ cm and also to the FD time sampling $Δt$⁠. We assume that the near-surface elastic parameters of the rock are constant and the same as those of the background medium. To obtain the required P- and S-wave velocities corresponding to the background rock medium, the wavefield injection method without an IABC can be used as explained in detail in  Appendix D. Note that in our laboratory, data cannot be acquired over a small area of the rock's surface due to its contact with the posts of the stand on which it rests. We chose not to carry out any numerical wavefield injection for these missing-data regions in the FD simulations. See  Appendix E for a detailed illustration of the influence due to this choice.

Figure 10 shows the resulting wavefield injection of the experimental data along the transparent, closed surface $Ssep$⁠, without the internal absorbing boundary. The corresponding movie is Mm. 3. As discussed in Sec. II D, the employed FD model consists of a homogeneous domain, with the material properties estimated by the procedure described above due to the fact that the exact properties and dimensions of the refilled hole are unknown. Despite the fact that we have used a method to optimize the P and S wave velocity estimates, Fig. 10 still shows a lot of unexpected events emerging from the wavefield injection surface $Ssep$ at later times (especially at t = 0.32 and 0.38 ms), following the reconstruction of the primary outgoing P- and S-waves. These “leaked waves” are caused by the inaccurate elastic model used in the FD simulation, which cannot precisely match the physical model of the actual experimental domain in the laboratory. Also, some of the artefacts could stem from ignoring the attenuation of the elastic waves within the granite medium, which is difficult to account for in the FD model used for wavefield injection.

Finally, Fig. 11 shows the wavefield separation of the experimental data using the wavefield injection method in conjunction with the proposed IABC. The corresponding movie is Mm. 4. Again, the IABC enables the reconstruction of the higher-order outgoing waves at the injection surface $Ssep$⁠, compared to not using the IABC as in Fig. 10. Further, taking the wavefield separation of the synthetic counterpart as a reference (comparing to both the arrival times and radiation characteristics of the waves as shown in Fig. 7 and also in Fig. 6), we can clearly identify the primary outgoing P wave that is separated at an early time of the simulation (⁠$t=0.23$ ms) in Fig. 11, especially above the top face of the injection surface $Ssep$⁠. The wavefield injection across the bottom face of $Ssep$ also separates the primary P wave that has propagated through the refilled hole (i.e., hole-related P) in the physical experiment. This separated transmitted wave energy also contains shear modes (i.e., hole-related S) in the form of a ringing tail with significantly larger amplitude compared to other separated wavefield components. Our interpretation is that the ringing tail represents the imprint of the (complex) interaction of the primary outgoing elastic waves, as emitted by the source, with the epoxy in the refilled hole. The large amplitudes of the downgoing waves reconstructed across the bottom face are thus most likely due to the fact that the epoxy used to refill the hole is softer than the granite rock (the exact properties of the epoxy are unknown) and therefore that the piezoelectric source embedded in the epoxy emits more energy to its compliant back side (i.e., into the epoxy) than into the less compliant front side (i.e., into the rock). At later times $t= 0.32$ and 0.38 ms (third and fourth rows of Fig. 11), the primary outgoing S-wave energy (first S) is also reconstructed and separated, as can be identified from the shear-wave radiation pattern (due to the direction of the f_z source and the particular elastic wavefield quantity v_z shown), mirroring the synthetic results shown in Fig. 7. Also, in Fig. 11, secondary outgoing S waves (second S), instead of leaked artefacts can clearly be identified once they are radiated outside of the injection surface $Ssep$ at later times (⁠$t= 0.32$ and 0.38 ms).

We tried to simulate a hole-like anomaly in the synthetic experiments and assessed the influence of not knowing the anomaly on the wavefield separation. We found that such a small-sized elastic anomaly on its own is not enough to account for all the hole-related imprints as observed in Fig. 11, and, in particular, such a model does not explain the significantly higher amplitudes observed on the bottom face across the refilled hole. This leads us to believe that there could be some non-linearity involved—a hypothesis that has been further supported by the observation of strong harmonics in the spectra of individual measurements at the top surface of the refilled hole (at the bottom face of the rock cube) (Guyer and Johnson, 2009).

Different FD simulations have been presented to demonstrate the IABC methodology in wavefield separation, corresponding to synthetic or experimental free-surface data. Figure 12 further compares these separated outgoing wavefields as obtained along lines just above the top face of the wavefield injection surface $Ssep$⁠. Without using the IABC in the wavefield separation, the wave energy radiating outside the top surface $Ssep$ differs significantly from its synthetic counterpart, suggesting that the first-order wavefield separation can be problematic due to its requirement that we know the interior medium properties. The use of the IABC is key for the wavefield separation such that the “leaked” artefacts outside of the injection surface $Ssep$ are removed and that all higher-order incident wavefields are correctly separated.

An effective IABC is much more difficult to construct compared to the more conventional absorbing boundary condition (ABC) for external boundaries used for attenuating outgoing waves in numerical modeling. Among the external ABCs, using passive absorbing layers (e.g., C-PMLs) (see Komatitsch and Martin, 2007) is more popular than using active wavefield cancellation (e.g., exact nonreflecting boundaries; see Givoli and Cohen, 1995). However, for an internal ABC, incorporating the passive layers (i.e., C-PMLs) only inside the wavefield injection surface $Ssep$ does not work well for attenuating inward-propagating waves. This is because the number of FD grid points around the corners and edges of the interior domain is naturally limited compared to those in the exterior domain, and in that case, waves cannot be effectively attenuated using passive techniques. This problem can be almost overcome using an active absorbing boundary such as described in this paper. The only (minor) issue remaining is related to the geometry of the IABC, involving injection, recording, and emitting surfaces that are each separated by a single FD grid point (for second-order accurate implementations), as illustrated in Fig. 3. In this case, not all ingoing waves can be suppressed by the active boundary condition as some waves can still travel through the corners of this two-grid space to the other side of the injection surface $Ssep$⁠.

The application of the current methods to elastic immersive wave experimentation forms the main motivation for this work. In immersive experimentation, a physical experiment is immersed into a numerical, virtual environment by spanning densely spaced sources at the boundary of the experimental domain. These sources can emit waves such that the outward-propagating waves in the physical experiment are canceled while the physical-to-virtual interaction can be synthesized at the free surface enclosing the experimental domain. To cancel unwanted boundary reflections at the free surface of a solid object (rock cube), the source signatures or voltages needed for these boundary sources controlling the physical boundary condition of the experimental domain correspond to the normal tractions of the primary outgoing waves only, which can be obtained using wavefield injection of the free-surface data (Thomsen et al., 2021). However, although the primary wavefield can be isolated using the wavefield injection method on its own, this isolation relies entirely on knowledge of the interior property map of the target object, which is generally unavailable. In contrast, the proposed IABC methodology allows wavefield separation without knowing any information about the interior of the object, while all of the outgoing parts of the wavefields in the free-surface data, including those that have interacted with the free surface multiple times (higher-order interactions), will be separated. In this case, using the IABC for elastic immersive wave experimentation would further require the separation of first-order component from the separated outgoing wavefields through, e.g., a windowing approach. Hence, the IABC provides a practical way to build elastic immersive experiments using the free-surface recorded data only. The development of immersive experiments is ongoing work in our laboratory.

Finally, many of the problems associated with wavefield separation of the real experimental free-surface data arose from the desire to place a physical source inside the solid target object. Drilling a hole to install such a source inside will unavoidably destroy or cause some damage to the object under study. A more effective solution could involve creating a virtual source inside the rock volume by projecting elastic waves from the free surface of the solid experimental domain. These waves can be generated using arrays of piezoelectric sources placed on the free-surface boundary such that the emitted wavefields can focus at a single point inside the rock to act as a virtual source. Such a method has been explored in acoustic time-reversal experiments (Derode et al., 1995).

We propose an IABC for wavefield separation of experimental data recorded on the closed free surface surrounding an object. In a FD simulation, injecting data from such an experiment along a closed, transparent surface results in ingoing and outgoing wavefield constituents of the data radiating in their respective correct directions. However, for a closed surface, the separated ingoing waves propagating in the numerical domain cause a problem as they travel to the other sides of the closed wavefield injection surface and destructively interfere with the data injection itself. This interference prevents the separation of higher-order outgoing wavefield constituents from the recorded data. To solve this problem, internal absorbing boundaries are placed inside the wavefield injection surface to prevent the propagation of the problematic ingoing waves in the FD simulation such that all orders of outgoing waves can be reconstructed and separated at the injection surface. The separated outgoing wavefields were validated on synthetic data by kinematically tracing the wave fronts of P and S components that represent the first- and second-order interactions with the free surface.

For the LDV experiments carried out for a cube of granite rock, the IABC is applied to decompose the experimental data in a homogeneous FD model whose P and S velocities are found using the wavefield injection method. Using the IABCs allows all orders of outgoing wavefields to be retrieved from the data. In particular, these wavefield constituents contain the imprint of the refilled hole inside the rock, which generates strong-amplitude waves with a ringing tail and preferentially more S-wave energy than P-wave energy due to the soft epoxy used to fill the hole.

This development of internal absorbing boundary conditions demonstrates the significance of leveraging the power of numerical modeling for physical modeling, especially for processing data acquired in a laboratory.

We thank two anonymous reviewers for their insightful comments. We would like to thank Dr. Claudio Madonna for measuring the wavespeeds and density of the drilled granite core in the Rock Physics and Mechanics Laboratory at ETH Zürich. We thank Dr. Pascal Edme for the discussion about determining the P and S velocities of the rock. We thank Christoph Bärlocher for setting up the lift system in the laboratory. X.L. thanks Thomas Haag for the education on using the Scanning LDV system and filling the hole of the rock. X.L. thanks Dr. Henrik Thomsen for discussing closed-aperture wavefield separation. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 694407).

A significant challenge for implementing an elastic IABC in time-domain FD simulations is the computational and memory storage cost of updating the active boundary condition, which is non-local in both time and space (Broggini et al., 2017; van Manen et al., 2020). The computation of the wavefield extrapolation integrals [Eqs. (5) and (6)] is carried out at each time step of the FD simulation, and the GFs needed are pre-computed prior to performing the simulation. In the FD simulation, the source and recording positions of these GFs correspond to each grid point on the recording surface $Srec$ and each grid point on the emitting surface $Semt$⁠, and these elastodynamic GFs involve multi-component source and receiver fields. Hence, memory cost of storing these (3D) GFs is estimated as $Nemt·Nrec·Nt·6·6·4$⁠, where $Nemt$ represents the number of sources (for canceling unwanted ingoing waves) on $Semt, Nrec$ represents the number of recording points on $Srec$⁠, N_t represents total time steps of the FD simulation, the first number 6 is assigned for the six types of source fields (including body force source $fx,fy,fz$ and three-component normally oriented deformation rate vector source) of the GFs, the second number 6 is assigned for the six types of receiver fields (including particle velocity $vx,vy,vz$⁠, and three-component normal traction vector) of the GFs, and the number 4 is assigned for the 4-byte (4B) memory occupied by a single-point floating number.

A simple estimation of memory storage cost required in the FD simulation can be made from the parameters given in Table I. We obtain the number of FD grid points across the recording surface $Srec$⁠, $Nrec=77 976$ (recording points); the number of FD grid points across the emitting surface $Semt$⁠, $Nemt=75 264$ (sources); and the time step, N_t = 720 (corresponding to $tmax=0.5$ ms). This results in the (dynamic) memory required for storing the GFs, $∼600$ TB, which goes far beyond the memory capacity of a personal laptop (typically 8–64 GB) or even on a single compute node of a supercomputer (64 GB to 4 TB). Even if the message passing interface (MPI) technology can be used to run the FD simulation across multiple compute nodes, the 3D numerical simulation performed for wavefield separation would require $∼9500$ nodes, each of which, for example, has 64 GB of dynamic memory on the ETH Euler cluster.

However, in the theory for the active boundary condition in the IABC, the following two facts can be relied on to reduce the memory storage cost of GFs:

1. The GFs used in Eqs. (5) and (6) are associated with a free-space homogeneous model. The computed GF is invariant to the source and receiver locations of the GF for constant source-receiver offsets (along x, y, and z directions).

2. For a 3D cube, around 5/6 of the GFs do not need to be computed and stored due to the scaling factor α in Eq. (4). These GFs are set to zero (by multiplying α) so that transmitted waves across the emitting surface $Semt$ are not extrapolated.

Point (1) is associated with the so-called translation invariance of the GFs (van Manen et al., 2020), which can be used to compute and store the GFs for each source-receiver offset rather than for each source and receiver location. Point (2) means that memory storage cost can be further reduced by a factor of 5/6. These homogeneous GFs can be compressed without a loss in accuracy when they are used to compute the wavefield extrapolation integrals. After applying this compression, the memory cost for storing the GFs is reduced to $∼84$ GB (from $∼600$ TB).

Figure 2 shows the damping layers placed inside the active boundary $Semt$ in the IABC. These absorbing layers involve the damping profiles (e.g., d_x) made from the C-PMLs, which are commonly used for external passive absorbing boundary conditions such that a numerically simulated wavefield over a finite-size domain can resemble waves propagating in unbounded media. We refer to Komatitsch and Martin (2007) for detailed definition and derivation of C-PMLs, and here, we illustrate how to modify the original C-PMLs for the IABC.

We first briefly describe the C-PML technique that has been used for an external absorbing boundary condition. Wave propagation in a 3D space can be decomposed into plane wave components that are independently propagating along x, y, and z directions in a Cartesian coordinate. In this case, absorbing layers are placed normal to the propagation direction of each plane wave, for instance, along the x axis to the positive direction, as shown in Fig. 13, and a damping profile $dx(x)$ is set up to attenuate the incoming waves impinging on the C-PML region from the regular simulation domain (Collino and Tsogka, 2001). In this schematic, we set $dx(x)=0$ for $x<0$ in the regular domain and $dx(x)>0$ for $x> 0$ in the C-PML region. Other damping profiles $dy(y)$ and $dz(z)$ for the plane waves propagating in the positive direction follow the same principle. For plane waves propagating in a negative axis direction, the C-PML region is still placed outside the regular domain but is aligned in a reversed profile to the negative direction. The C-PML, when used in practical FD modeling, will introduce profile parameters a and b, both of which were described by a generic parameter d as in Figs. 2 and 13 (Komatitsch and Martin, 2007).

To absorb ingoing waves toward the interior of the emitting surface $Semt$⁠, we reverse all (conventional) damping profiles including a, b along both positive and negative directions of each axis (x, y, z) and incorporate the profiles into the domain enclosed by the surface $Semt$⁠, as shown in Fig. 2. Also, to maximize the efficiency of wavefield attenuation using the C-PML layers, each single profile aligned in a positive or negative direction occupies about half the extent of the cubic domain enclosed by $Semt$⁠. Figure 14 shows the interior C-PML profiles, together with the exterior C-PML profiles padded around the regular numerical domain to absorb the separated outgoing waves in the FD simulation.

To illustrate the problem of using an inexact FD model in the injection-based wavefield separation, we place a cuboid scatterer in the synthetic wave propagation experiment, with the corresponding model shown in Fig. 15(a). We first inject the free-surface data corresponding to the heterogeneous interior into a homogeneous FD model that only involves the background medium (the same as in Fig. 6). Figure 16 (top rows) shows this wavefield separation, with artefacts appearing to “leak” out of the injection surface $Ssep$ at later times (⁠$t=0.40$ and 0.50 ms) of the FD simulation. These artefacts are erroneous and contaminate the primary outgoing P and S wavefields retrieved outside the injection surface $Ssep$⁠.

The “leaked” artefacts observed in the first and second rows of Fig. 16 are caused by the missing cuboid scatterer that should be included in the model used for the FD simulation. The fourth and fifth rows of Fig. 16 show the wavefield separation when the interior of the FD model does include the cuboid scatterer and exactly matches the model used to generate the free-surface data; in this scenario, the reconstructed ingoing wavefields inside the injection surface $Ssep$ propagate exactly as the waves reflected by the free surface in the corresponding synthetic experiment. Hence, these interior-crossing waves destructively interfere with the outgoing part of the second- and higher-order free-surface interactions on the other sides, removing the source of the wave energy “leaking” out of the injection surface $Ssep$⁠. Therefore, only primary outgoing waves are obtained outside the surface $Ssep$⁠. Note that some small-amplitude wave energy still comes out of the injection surface $Ssep$ at later times of the FD simulation, as shown in the fourth row of Fig. 16. These waves are not artefacts but correspond to the primary outgoing energy that is multiply scattered by the interior cuboid inhomogeneity before arriving at the free surface. However, this multiply scattered later-arrived primary wave energy cannot be correctly retrieved in a laboratory when using the wavefield injection method only, as knowing a true elastic model is neither practical nor desirable.

The wavefield injection method can be used to optimize the numerical model parameters used in the FD simulations—the P and S wave speeds corresponding to the background rock medium. When injecting the free-surface data in the FD simulation, only the primary incident outgoing P and S waves should radiate away from the wavefield injection surface $Ssep$⁠, while all of the reconstructed ingoing waves should be perfectly confined inside the surface $Ssep$ without any energy leakage. However, this ideal case requires that the true elastic model is used inside the wavefield injection surface $Ssep$ of the FD simulation. If the numerical model used in the FD simulation is erroneous compared to the true model, the wavefield separation would be incorrect. In particular, two observations can be exploited: (1) the primary P- and S-wave fronts observed on orthogonal slices will deviate from a purely spherical wave front, or its circular projections (for an assumed homogeneous, isotropic medium) if the incorrect velocities are used. Therefore, the primary propagation velocities can be estimated by maximizing the circularity of the separated wave fronts. (2) Based on the fact that wavefield injection should only reconstruct primary outgoing waves outside of the surface $Ssep$ without significant “leaked” energy, we can estimate the velocities by minimizing the amount of such leaked energy at later times.

Figure 17 illustrates the effects of varying the shear wave velocity (V_S) in the model of the FD simulation when separating the synthetic homogeneous free-surface data (Fig. 5), and these effects motivate the two approaches described above. When varying V_S, the wave front of the separated primary S wave (first S) will be distorted compared to the (perfect) spherical wave front. Varying the P-wave velocity has a similar effect (not shown here). To illustrate the second approach, varying the S-wave velocity has a significant impact on the number of leaked artefacts outside the wavefield injection surface $Ssep$⁠, as shown in Fig. 17.

However, when we applied the first approach (i.e., optimizing the circularity of wave fronts) before decomposing the experimental data, a significant feature makes a difference: the wave front of the primary separated (S) wavefield is elliptical, as observed in Fig. 10 (first column). Interestingly, we found that this ellipticity cannot be reduced by varying the primary S-wave propagation velocities. The reason for this ellipticity is due to the fact that the piezoelectric source, embedded inside the rock volume with the surrounding epoxy, has a finite dimension, i.e., a cuboid with 8 mm in height and 1 cm in width and length (Li et al., 2019). Hence in our practices of optimizing the numerical model parameters, the second approach gave the best results, whose corresponding parameters are used in the wavefield separation of the laboratory free-surface data (Table I).

In the laboratory, the experimental rock volume is placed on the four posts of a height-adjustable metal frame, whose junctions with the rock's surface are unavoidably masked; LDV measurements in these regions are not possible, as shown in Fig. 4(a), lower inset. One straightforward method involves extrapolating the particle velocity wavefields to complement the missing data, but the extrapolated wavefields are not sufficient for wavefield separation at these missing-data regions (Koene and Robertsson, 2018; Robertsson et al., 2015). This is because the free-surface boundary condition involved in the theory of the wavefield injection technique (Sec. II B) does not hold at these contact areas, which tightly contact with the rubber at the ends of the posts. Hence, we chose not to carry out any numerical wavefield injection at those parts of the wavefield injection surface, which correspond to the (four) missing-data areas in the physical experiments. However, these areas are in the shape of 90° sectors ∼7–8 cm in radius; this size is of the same order of magnitude as the wavelengths of S waves (∼4–12 cm). Hence, our choice will unavoidably result in an undesired effect in the wavefield decomposition.

To study the effect of not having a full-aperture injection surface $Ssep$ in the FD simulation, we again carried out the wavefield decomposition using the synthetic (homogeneous) free-surface data with no wavefield injection in these rubber-masked regions (located as in the laboratory counterparts). Figure 18 shows this injection result, which can be compared to the case when using the complete synthetic data as in Fig. 6. The error caused by missing the injection around the lower four corners of the cube $Ssep$ is small compared to the separated primary outgoing P and S waves. Also, the influence of this error is local in space in the separated outgoing wavefields observed outside the injection surface $Ssep$⁠, with a relatively larger effect close to the no-data-injection regions and a negligible effect above the top face of $Ssep$⁠. Note that using the IABC in the FD simulation can further absorb the ingoing part of the undesired erroneous waves caused by not having injection in the rubber-masked regions.