A wave propagation laboratory is proposed which enables the study of the interaction of broadband signals with complex materials. A physical experiment is dynamically linked to a numerical simulation in real time through transmitting and recording transducer surfaces surrounding the target. The numerical simulation represents an arbitrarily larger domain, allowing experiments to be performed in a total environment much greater than the laboratory experiment itself. Specific applications include the study of non-linear effects or wave propagation in media where the physics of wave propagation is not well understood such as the effect of fine scale heterogeneity on broadband propagating waves.

Laboratory-based physical wavefield modeling is useful for studying wave propagation in complex realistic three-dimensional materials under controlled conditions. Physical modeling was particularly popular before the 1990s when computer power was not sufficient to model realistic 3D experiments numerically (McDonald et al., 1983). With increasing computer power, synthetic modeling of complex wave propagation in earth models ranging from acoustic to anisotropic viscoelastic media has become feasible, and largely replaced the use of physical modeling.

Current laboratory experimentation is mainly aimed at noise interferometry studies, scattered surface wave analyses, fracture characterization, acoustic spectroscopy, and full waveform inversion (Hadziioannou et al., 2009; Blum et al., 2011; Larose et al., 2010; Mikesell et al., 2012). The limited size of these experiments generally requires the use of wavelengths that are orders of magnitude smaller than the laboratory target, in order to avoid interference of boundary reflections from the edges. Therefore, the use of most current wave propagation laboratories is limited to the study of frequencies that are orders of magnitude higher than the bandwidth of for example the seismological or acoustic real-world scenarios that they represent.

We propose a new concept for wave propagation laboratories in which the physical experiment is linked with a numerical simulation in real time, exploiting novel theory of exact boundary conditions (van Manen et al., 2007). The physical experiment and the numerical simulation interact at discrete time steps, so that the physical experiment is fully immersed in the larger numerical simulation domain, and the two domains interact with and drive each other. The interaction between the two domains takes place at the edges of the physical experiment where specific boundary conditions cause undesired boundary reflections to be canceled. This allows the propagation of low frequency signals to be simulated in targets of any size, including sub-wavelength scale. Furthermore, wave interactions with medium structures outside of the physical laboratory but represented within the numerical domain are automatically incorporated. The methodology we use has analogies to wave-equation based seismic imaging methods that eliminate specific parts of the wavefield (Amundsen, 2001) and to work on acoustic wave cancellation and synthesis (Berkhout et al., 1993; Thomson, 2012).

The proposed laboratory will be particularly suited to study the propagation of broadband waves containing considerably lower frequencies than currently possible, and the interaction of laboratory waves with a variety of modeled structures and boundary conditions that surround the target. Applications include the verification of homogenization theory and procedures to upscale non-periodic media (Capdeville et al., 2010), the study of frequency dependent phenomena such as velocity dispersion and attenuation, and the study of non-linear effects such as source and receiver interactions with the target medium.

Our goal is to link a physical experiment with a numerical simulation in real time. The numerical simulation represents the wavefield propagating in a virtual external environment surrounding the physical laboratory target. The physical experiment is surrounded by discrete transducers acting as a boundary of secondary sources that cancels outgoing signals and simultaneously creates signals that enter the physical domain from outside (from the numerical domain). Transducers simultaneously record outgoing waves which enter the numerical simulation, and which may scatter back into the physical simulation at later times.

Figure 1 shows a cartoon of the proposed physical laboratory and the complementary numerical simulation domain. The key in connecting these two experiments, is to include an overlapping region in which the media in both domains are identical. Outside of this linking region, the media and their Green's functions can be arbitrarily different and may even be governed by completely different physics of wave propagation. Two surfaces are used to link the two experiments: the recording surface S rec and the emitting surface S emt . The emitting surface coincides with the bounding surface of the physical laboratory, which has known boundary conditions.

FIG. 1.

Cartoon of physical laboratory and numerical simulation domains. The source distribution on S emt cancels physical boundary reflections for outgoing waves (1), and introduces ingoing waves from interactions with the numerical model (2). In the overlapping region the medium of the physical model must be equal to that of the numerical model.

FIG. 1.

Cartoon of physical laboratory and numerical simulation domains. The source distribution on S emt cancels physical boundary reflections for outgoing waves (1), and introduces ingoing waves from interactions with the numerical model (2). In the overlapping region the medium of the physical model must be equal to that of the numerical model.

Close modal

To derive the quantities to inject on the emitting boundary S emt we use a representation for the acoustic pressure in a volume V enclosed by a surface S. The pressure p ( x , t ) in V can be written as (Fokkema and van den Berg, 1993; Aki and Richards, 2002)

(1)

where v ( x , t ) is the particle velocity, s ( x , t ) is the source-time function, G ( x , t , x , 0 ) and Γ ( x , t , x , 0 ) are the medium pressure responses at x due to monopole and dipole sources at x , respectively, and n denotes the inward pointing normal to the surface. We employ this representation with the physical laboratory being the volume V enclosed by the surface S emt . Two different states are compared: the desired state in which the physical and numerical domains are linked (state A) and the actual state in the laboratory, in which only the physical domain exists (state B). The medium properties and source function are equal in both states, but different boundary conditions are specified on S emt . Note that we are free to choose the boundary conditions for the Green's functions in the representation given by Eq. (1). In the following, we show the derivation for the case in which the physical experiment in the laboratory has a rigid boundary. However, similar reasoning can be followed for other types of boundary conditions.

In state A the two linked domains should act as one larger domain, so that S emt is only an imagined surface on which the wavefield is continuous and no other boundary conditions are specified or imposed. State B corresponds to the physical experiment in which S emt is a rigid boundary. Table I summarizes the two states and their boundary conditions on S emt .

TABLE I.

Description of the states employed for the derivation of the source distribution on S emt .

  State A (linked domains) State B (physical domain)
Field  p A ( x , t ) , v A ( x , t )   p B ( x , t ) , v B ( x , t )  
Source  s A ( x , t ) = s ( x , t )   s B ( x , t ) = s ( x , t )  
Boundary conditions on S emt   None  v B | S emt n = 0  
  State A (linked domains) State B (physical domain)
Field  p A ( x , t ) , v A ( x , t )   p B ( x , t ) , v B ( x , t )  
Source  s A ( x , t ) = s ( x , t )   s B ( x , t ) = s ( x , t )  
Boundary conditions on S emt   None  v B | S emt n = 0  

Following Eq. (1) and choosing the Green's functions of state B, these boundary conditions allow us to write the pressure in state A inside V as

(2)

because the term containing Γ B ( x , t τ , x , 0 ) goes to zero due to the boundary condition in state B. Similarly, the pressure in state B inside V can be written as

(3)

without any surface integral contributions, because the normal component of the particle velocity v B on S emt goes to zero there as well. Comparing Eq. (2) to Eq. (3), the difference between the desired (linked) state A and the laboratory physical state B is the term

(4)

which corresponds to a distribution of monopole sources on S emt with source strength v A n . A critical requirement in this procedure is therefore to inject a term weighted by the particle velocity v A [Eq. (4)] simultaneously with the energy response for this velocity actually arriving at the emitting surface. We therefore need a separate recording surface ( S rec ) from which we can extrapolate the recorded wavefield to S emt in order to predict the arriving energy ahead of time. This procedure is only possible if the media between S rec and S emt are identical in the numerical and physical domain. We use a time discrete version of Green's second identity similar to Eq. (4) of van Manen et al. (2007),

(5)

where ̂ denotes discrete time sampling. G ̂ A and Γ ̂ A are the particle velocity responses due to a monopole and a dipole impulse function at S rec , respectively. These functions are computed in advance through numerical wave propagation simulations on a model corresponding to the numerical domain outside S rec . Inside S rec the structure of the model is irrelevant since the total effect of the contributions in Eq. (5) is zero (Fokkema and van den Berg, 1993; Thomson, 2012). Note that this means that the notation in Eq. (5) is actually too strict, since the medium for which the Green's functions are calculated only needs to match state A outside S rec . At each discrete time step n these Green's functions are used to extrapolate the wavefield from S rec to S emt for all future time steps l > n . At each time step a discrete time series is added to the buffer v ̂ n ( x emt , l , n ) that then already contains future values to be used for the emission on S emt at later time steps (van Manen et al., 2007). In principle the Green's functions are not restricted to numerically simulated responses and physically measured Green's functions can also be used.

Equation (5) links the physical measurements on S rec to the numerical simulation domain and provides the velocity term to be used in Eq. (4), so that the combination of Eq. (4) and Eq. (5) results in the final source distribution to be emitted at each time step. Thus, a physical laboratory with a rigid boundary can be linked to a numerical domain exploiting monopole sources on its surface. Similarly, a laboratory with a free surface and a surface source distribution consisting of dipole sources can be used.

The wave propagation laboratory outlined in this paper requires a purely acoustic medium between S rec and S emt . This has direct applications when modeling for instance (visco)acoustic wave propagation around (visco)acoustic or (an)elastic targets immersed in a purely acoustic medium such as air or water, for example, for characterizing the response of acoustic sources in water. However, the same methodology can also be generalized and extended beyond the acoustic wave equation to situations where wave propagation in the overlapping region is governed by the elastic or even electromagnetic wave equations.

The wave propagation laboratory is demonstrated through a numerical simulation representing the physical laboratory as a two-dimensional viscoacoustic model. The physical laboratory has dimensions 2 × 2 m and is surrounded by a rigid boundary. It is embedded in a larger numerical domain of dimensions 4 × 4 m, bounded by a free surface on all sides (to generate strong reflections as a demonstration). A pressure source is located at x = 1.5 m, y = 1.8 m, with the source function of a Ricker wavelet with a center frequency of 5 kHz. Figure 2 shows the two domains and their material properties. For demonstration purposes a simple homogeneous model with a P-wave velocity of 2000 m/s is used throughout. The physical model contains a strongly attenuating target, whereas the numerical domain is governed by purely acoustic wave propagation.

FIG. 2.

Velocity and Q model for numerical example. Wave propagation in a homogeneous medium with velocity 2000 m/s is modeled, with the physical laboratory containing a strongly attenuating anomaly (Q = 20). The emitting surface is a rigid boundary, and the numerical domain is bounded by free surfaces.

FIG. 2.

Velocity and Q model for numerical example. Wave propagation in a homogeneous medium with velocity 2000 m/s is modeled, with the physical laboratory containing a strongly attenuating anomaly (Q = 20). The emitting surface is a rigid boundary, and the numerical domain is bounded by free surfaces.

Close modal

The wave equation is solved on a staggered grid with second-order temporal and spatial accuracy. Attenuation is incorporated using a single memory variable (Robertsson et al., 1994) resulting in a Q value of 20 at the center frequency of the source wavelet. For the calculation of the Green's functions employed in Eq. (5), the medium properties inside S rec are irrelevant (as described above), so a purely acoustic homogeneous internal volume was used.

For this specific example, the pressure and the particle velocity are recorded on 324 locations along S rec , and 404 monopole sources on S emt are employed. The experiment is run for 500 time steps. This requires a look-up table for the Green's functions of size 2 × N rec × N emt × n t = 2 × 324 × 404 × 500 = 1.3 × 10 8 entries.

Figure 3 shows snapshots of the pressure field at t = 1 and 2 ms. A reference solution is calculated on the linked domain for comparison (left), where no boundary conditions exist on the surface S emt and finite difference modeling is carried out on the entire model (4 × 4 m). The result of the wave propagation in the physical laboratory modeled using Eqs. (4) and (5) (middle) is equal to that of the reference model (within numerical precision). For this laboratory simulation, modeling is carried out on the smaller (2 × 2 m) domain only, and all interactions with the exterior domain are accounted for by using the methodology described above with the pre-computed Green's functions from S rec to S emt . The right column shows the simulation of a physical wave propagation laboratory of equal dimensions but without active sources on its boundary. The snapshot at t = 1 ms demonstrates that the source arrangement on S emt cancels all boundary reflections that would otherwise exist from the edges of the laboratory. The snapshot at t = 2 ms shows that the extrapolation integral in Eq. (5) fully accounts for the interaction of the outgoing waves with the numerical domain and the excitation of resulting ingoing waves.

FIG. 3.

Snapshots of modeled pressure on a reference (complete) model (left), in the simulated proposed physical laboratory (middle), and in a physical laboratory without active boundary sources (right). The corresponding velocity and Q model are displayed in Fig. 2. Top: snapshots at t = 1 ms. Bottom: snapshots at t = 2 ms.

FIG. 3.

Snapshots of modeled pressure on a reference (complete) model (left), in the simulated proposed physical laboratory (middle), and in a physical laboratory without active boundary sources (right). The corresponding velocity and Q model are displayed in Fig. 2. Top: snapshots at t = 1 ms. Bottom: snapshots at t = 2 ms.

Close modal

Evidently, this methodology can also be used in pure numerical modeling to simulate wave propagation on a sub-domain of a larger model, while accounting for all interactions between this sub-domain and the larger model. The Green's functions required for the extrapolation through the larger model can be computed efficiently using interferometric modeling techniques (van Manen et al., 2005). Related concepts were recently proposed by Grote and Kirsch (2007) and Grote and Sim (2011).

We have demonstrated the concept for a new wave propagation laboratory in which a physical experiment is immersed in a larger numerical simulation domain. The physical experiment and the numerical simulation are linked through the use of exact dynamic boundary conditions. Sources distributed on the surface of the laboratory (the so-called emitting surface) cancel the boundary reflections from outgoing waves, and excite the ingoing waves that are a result of the interaction of outgoing waves with the larger numerical simulation domain.

The proposed laboratory will enable experiments with broadband signals in complex materials that cannot be modeled synthetically or in conventional physical laboratories. Since it eliminates undesired boundary reflections from the edges of the laboratory, the propagation of low frequency signals in targets of any size can be studied. This will facilitate the study of non-linear effects, and poorly understood physics of wave propagation such as the upscaling of non-periodic media.

We would like to thank Chris Chapman, Stewart Greenhalgh, and Colin Thomson for insightful discussions during the work.

1.
Aki
,
K.
, and
Richards
,
P. G.
(
2002
).
Quantitative Seismology
, 2nd ed. (
University Science Books
,
Sausalito, CA
).
2.
Amundsen
,
L.
(
2001
). “
Elimination of free-surface related multiples without need of the source wavelet
,”
Geophysics
66
,
327
341
.
3.
Berkhout
,
A. J.
,
de Vries
,
D.
, and
Vogel
,
P.
(
1993
). “
Acoustic control by wave field synthesis
,”
J. Acoust. Soc. Am.
93
,
2764
2778
.
4.
Blum
,
T. E.
,
van Wijk
,
K.
,
Snieder
,
R.
, and
Willis
,
M. E.
(
2011
). “
Laser excitation of a fracture source for elastic waves
,”
Phys. Rev. Lett.
107
,
275501
.
5.
Capdeville
,
Y.
,
Guillot
,
L.
, and
Marigo
,
J.-J.
(
2010
). “
2-D non-periodic homogenization to upscale elastic media for P–SV waves
,”
Geophys. J. Int.
182
,
903
922
.
6.
Fokkema
,
J. T.
, and
van den Berg
,
P. M.
(
1993
).
Seismic Applications of Acoustic Reciprocity
(
Elsevier
,
Amsterdam
).
7.
Grote
,
M. J.
, and
Kirsch
,
C.
(
2007
). “
Nonreflecting boundary condition for time-dependent multiple scattering
,”
J. Comput. Phys.
221
,
41
62
.
8.
Grote
,
M. J.
, and
Sim
,
I.
(
2011
). “
Local nonreflecting boundary condition for time-dependent multiple scattering
,”
J. Comput. Phys.
230
,
3135
3154
.
9.
Hadziioannou
,
C.
,
Larose
,
E.
,
Coutant
,
O.
,
Roux
,
P.
, and
Campillo
,
M.
(
2009
). “
Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: Laboratory experiments
,”
J. Acoust. Soc. Am.
125
,
3688
3695
.
10.
Larose
,
E.
,
Planes
,
T.
,
Rossetto
,
V.
, and
Margerin
,
L.
(
2010
). “
Locating a small change in a multiple scattering environment
,”
Appl. Phys. Lett.
96
,
204101
.
11.
McDonald
,
J. A.
,
Gardner
,
G. H. F.
, and
Hilterman
,
F. J.
(
1983
).
Seismic Studies in Physical Modeling
(
IHRDC Press
,
Boston
).
12.
Mikesell
,
T. D.
,
van Wijk
,
K.
,
Blum
,
T. E.
,
Snieder
,
R.
, and
Sato
,
H.
(
2012
). “
Analyzing the coda from correlating scattered surface waves
,”
J. Acoust. Soc. Am.
131
,
EL275
EL281
.
13.
Robertsson
,
J. O. A.
,
Blanch
,
J.
, and
Symes
,
W. W.
(
1994
). “
Viscoelastic finite-difference modeling
,”
Geophysics
59
,
1444
1456
.
14.
Thomson
,
C. J.
(
2012
). “
Research note: Internal/external seismic source wavefield separation and cancellation
,”
Geophys. Prospect.
60
,
581
587
.
15.
van Manen
,
D.-J.
,
Robertsson
,
J. O. A.
, and
Curtis
,
A.
(
2005
). “
Modeling of wave propagation in inhomogeneous media
,”
Phys. Rev. Lett.
94
,
1643011
.
16.
van Manen
,
D.-J.
,
Robertsson
,
J. O. A.
, and
Curtis
,
A.
(
2007
). “
Exact wave field simulation for finite-volume scattering problems
,”
J. Acoust. Soc. Am.
122
,
EL115
EL121
.