Marchenko-Based Target Replacement, Accounting for All Orders of Multiple Reflections.

In seismic monitoring, one is usually interested in the response of a changing target zone, embedded in a static inhomogeneous medium. We introduce an efficient method that predicts reflection responses at the Earth's surface for different target-zone scenarios, from a single reflection response at the surface and a model of the changing target zone. The proposed process consists of two main steps. In the first step, the response of the original target zone is removed from the reflection response, using the Marchenko method. In the second step, the modelled response of a new target zone is inserted between the overburden and underburden responses. The method fully accounts for all orders of multiple scattering and, in the elastodynamic case, for wave conversion. For monitoring purposes, only the second step needs to be repeated for each target-zone model. Since the target zone covers only a small part of the entire medium, the proposed method is much more efficient than repeated modelling of the entire reflection response.

Introduction

In seismic modelling, inversion, and monitoring one is often interested in the response of a relatively small target zone, embedded in a larger inhomogeneous medium. Yet, to obtain the seismic response of a target zone at the Earth's surface, the entire medium enclosing the target should be involved in the modelling process. This may become very inefficient when different scenarios for the target zone need to be evaluated or when a target that changes over time needs to be monitored, for example, to follow fluid flow in an aquifer, subsurface storage of waste products, or production of a hydrocarbon reservoir. Through the years, several efficient methods have been developed for modelling successive responses of a medium in which the parameters change only in a target zone. Robertsson and Chapman (2000) address this problem with the following approach. First, they model the wave field in the full medium, define a boundary around the target zone in which the changes take place, and evaluate the field at this boundary. Next, they numerically inject this field from the same boundary into different models of the target zone. Because the target zone usually covers only a small part of the full medium, this injection process takes only a fraction of the time that would be needed to model the field in the full medium. This method is very well suited to model different time‐lapse scenarios of a specific subsurface process in an efficient way. A limitation of the method is that multiple scattering between the changed target and the embedding medium is not taken into account. The method was adapted by van Manen et al. (2007) to account for this type of interaction, by modifying the field at the boundary around the changed target at every time‐step of the simulation. Wave field injection methods are not only useful for efficient numerical modelling of wave fields in a changing target zone but they can also be used to physically inject a field from a large numerical environment into a finite‐size physical model (Vasmel et al., 2013).

Instead of numerically modelling the field at the boundary enclosing the target, Elison et al. (2016) propose to use the Marchenko method to derive this field from reflection data at the surface. Hence, to obtain the wave field in a changing target zone, they need a measured reflection response at the surface of the original medium and a model of the target. Their method exploits an attractive property of the Marchenko method, namely, that “redatumed” reflection responses of a target zone from above ( ) and from below ( ) can both be obtained from single‐sided reflection data at the surface and an estimate of the direct arrivals between the surface and the target zone (Wapenaar et al., 2014a).

In most of the methods discussed above, the wave fields are derived inside the changing target. Here we discuss a method that predicts reflection responses (including all multiples) at the Earth's surface for different target‐zone scenarios, from a single reflection response at the surface and a model of the changing target zone. The proposed method, which we call target replacement, consists of two main steps. In the first step, which is analogous to the method proposed by Elison et al. (2016), we use the Marchenko method to remove the response of the target zone from the original reflection response. In the second step we insert the response of a new target zone, yielding the desired reflection response at the surface for the particular target‐zone scenario. Both steps fully account for multiple scattering between the target and the embedding medium. Note that, to model different reflection responses for different target models, only the second step needs to be repeated. Hence, this process is particularly efficient when reflection responses at the surface are needed for many target‐zone scenarios. Also note that, unlike the model‐driven methods of Robertsson and Chapman (2000) and van Manen et al. (2007), our method as well as that of Elison et al. (2016) only needs a smooth model of the overburden and no model of the underburden. The required detailed information of the over‐ and underburden comes from the measured reflection response.

Similar as the other methods discussed in this introduction, we assume that the target zone is the only region in which changes occur; the overburden and underburden are assumed to remain unchanged. However, changes in a reservoir may lead to changes in the embedding medium (Hatchell & Bourne, 2005; Herwanger & Horne, 2009). When this is the case, the target zone should not be restricted to the reservoir, but it should also include the part of the embedding medium in which the changes have a noticeable effect on the waves propagating through it. Of course, the larger the target zone, the smaller the efficiency gain.

The setup of this paper is as follows. In section 2, we derive a representation of the seismic reflection response at the Earth's surface (including all orders of multiple scattering), which explicitly distinguishes between the response of the target zone and that of the embedding medium. Next, based on this representation, in section 3, we discuss how to remove the response of the target zone from the reflection response at the surface. In section 4, we discuss how the response of a changed target zone can be inserted into the reflection response at the surface. The proposed method is illustrated with numerical examples in section 5. We end the paper with a discussion (section 6) and conclusions (section 7).

Representation of the Reflection Response

We derive a representation for the reflection response at the Earth's surface, which distinguishes between the response of the target zone and that of the embedding medium. We start by dividing the subsurface into three units. The first unit, indicated as unit a in Figure 1, covers the region between the Earth's surface and boundary , the latter defining the upper boundary of the target zone. The Earth's surface (indicated by the solid line) may be considered either as a free or as a transparent surface (the latter after surface‐related multiple elimination). The Earth's surface is included in unit a. A transparent boundary (indicated by the upper dashed line) is defined at an infinitesimal distance below the Earth's surface (in the following we abbreviate “an infinitesimal distance above/below” as “just above/below”). Unit a, that is, the region above the target zone, is called the overburden. The second unit, indicated as unit b in Figure 1, represents the target zone and is enclosed by boundaries and . The third unit, indicated as unit c in Figure 1, represents the region below the lower boundary of the target zone, . Unit c, that is, the region below the target zone, is called the underburden.

Figure 1

Subdivision of the inhomogeneous subsurface into three units: an overburden (unit a), a target zone (unit b), and an underburden (unit c). Note that unit a includes the Earth's surface just above . This surface may be considered either as a free or as a transparent surface.

We assume that the media inside the units are arbitrary inhomogeneous, lossless media. Furthermore, we assume that the boundaries and do not coincide with interfaces, or in other words, we consider these boundaries to be transparent for downgoing and upgoing waves incident to these boundaries. The representation derived below could be extended to account for scattering at these boundaries, but that would go at the cost of clarity. By allowing some flexibility in the definition of the target zone, it will often be possible to choose boundaries and that are (close to) transparent.

The starting point for the derivation of the representation and the target replacement method is formed by the following one‐way reciprocity theorems in the space‐frequency domain and

(Wapenaar & Grimbergen, 1996). Here and can stand for any of the boundaries , , and . Subscripts A and B refer to two independent states. Superscripts + and − stand for downward and upward propagation, respectively. Superscript t in equation (1) denotes the transpose and superscript † in equation (2) the adjoint (i.e., the complex conjugate transpose). The vectors and represent flux‐normalized one‐way wave fields in states A and B. For the elastodynamic situation they are defined as where , , and represent P, S1, and S2 waves, respectively. For the acoustic situation, and reduce to scalar functions. The Cartesian coordinate vector x is defined as x = (x 1,x 2,x 3) (the x 3‐axis pointing downward), and ω denotes angular frequency. An underlying assumption for both reciprocity theorems is that the medium parameters in states A and B are identical in the domain enclosed by boundaries and . Outside this domain the medium parameters in state A may be different from those in state B, a property that we will make frequently use of throughout this paper. Another assumption is that there are no sources between and . Finally, an assumption that holds specifically for equation (2) is that evanescent waves are neglected at boundaries and . For a more detailed discussion of these one‐way reciprocity theorems, including their extensions for the situation that the domain between and contains sources and the medium parameters in the two states are different in this domain, see Wapenaar and Grimbergen (1996).

In the following derivations, equations (1) and (2) will frequently be applied, each time to a combination of independent wave states in two media that are identical in the domain between and . Figure 2 shows six media that will be used in different combinations. Media a, b, and c in the left column contain the units a (the overburden), b (the target zone), and c (the underburden) of the actual medium, each embedded in a homogeneous background. The gray areas indicate the inhomogeneous units (as depicted in Figure 1), whereas the white areas represent the homogeneous embedding. Reflection and transmission responses are also indicated in Figure 2. Reflection responses from above and from below are denoted by and , respectively, and the transmission responses by T + and T −. The subscripts a, b, and c refer to the units to which these responses belong. The rays are simplifications of the actual responses, which contain all orders of multiple scattering and, in the elastodynamic case, mode conversion. When the Earth's surface just above is a free surface, then the responses in unit a also include multiple scattering related to the free surface. Media A, B, and C in the right column in Figure 2 consist of one to three units, as indicated (note that medium A is identical to medium a, whereas medium C represents the entire medium). The reflection and transmission responses of these media are indicated by capital subscripts A, B, and C. In addition, the Green's functions G +,+ and G −,+ in these media between and the top boundary of the deepest unit are shown (the superscripts will be explained later). Again, all responses contain all orders of multiple scattering (and mode conversion), including surface‐related multiples when there is a free surface just above .

Figure 2

Six media with their responses. Gray areas represent the inhomogeneous units (and combinations thereof) of Figure 1. Media A (=a), B, and C include the Earth's surface just above , which may be considered either as a free or as a transparent surface. The rays stand for the full responses, including all orders of multiple scattering and, in the elastodynamic case, mode conversion.

Our aim is to derive a representation for the reflection response of the entire medium, , in terms of the reflection responses of media A (=a), b, and c. We start by deriving a representation for in terms of the reflection responses of media A and b. To this end, we substitute the quantities of Table 1 into equation (1). Let us first discuss these quantities one by one. In state B, the downgoing and upgoing fields in medium B for x at are given by

Table 1

Quantities to Derive a Representation for

	State A:	State B:
	Medium A	Medium B
	Source at x R just above S_0	Source at x S just above S_0
S_0	pA+(x,ω)→Iδ(x_H−x_H,R)	pB+(x,ω)→Iδ(x_H−x_H,S)
	+r^∩RA∪(x,x_R,ω)	+r^∩RB∪(x,x_S,ω)
	pA−(x,ω)→RA∪(x,x_R,ω)	pB−(x,ω)→RB∪(x,x_S,ω)
S_1	pA+(x,ω)→TA+(x,x_R,ω)	pB+(x,ω)→GB+,+(x,x_S,ω)
	pA−(x,ω)→O	pB−(x,ω)→GB−,+(x,x_S,ω)

Here is the Green's one‐way wave field matrix in medium B in the space‐frequency domain (Wapenaar, 1996). The source is at x , which is chosen just above . The second superscript + indicates that this source is downward radiating. The receiver is at x at . The first superscript ± indicates the propagation direction at the receiver (+ for downgoing and − for upgoing). Analogous to equation (3), the general Green's one‐way wave field matrix can, for the elastodynamic situation, be written as

Each column corresponds to a specific type of source at and each row to a specific type of receiver at x (where subscripts ϕ, ψ, and υ refer to flux‐normalized P, S1, and S2 waves, respectively). For the acoustic situation, reduces to a scalar function. The following reciprocity relations hold for the general Green's matrix

(Haines, 1988; Kennett et al., 1990; Wapenaar, 1996). In state B, the upgoing field for x at in Table 1 is given by

Note that represents a reflection response from above, denoted by , whenever the source and receiver are situated at (or just above) the same depth level. From equations (6) and (9), we find

Similarly, represents a reflection response from below, denoted by , whenever the source and receiver are situated at (or just below) the same depth level. From equations (7) and (9) we find

In state B, the downgoing field for x at in Table 1 is given by

Since x was chosen just above , the direct contribution of the flux‐normalized Green's matrix consists of a spatial delta function δ(x H−x H,), with x H=(x 1,x 2) and x H,=(x 1,,x 2,); hence, the singularity occurs at the lateral position of the source. This delta function is multiplied by I, which is a 3 × 3 identity matrix for the elastodynamic situation, to acknowledge the matrix character of , as defined in equation (5). For the acoustic situation I = 1. The second term in equation (12), , accounts for the Earth's surface just above . Here is the reflection operator of the Earth's surface from below. It turns the reflection response into a downgoing field, which, according to equation (12), is added to the direct downgoing field. When the Earth's surface is transparent, we may simply set , where O is a 3 × 3 zero matrix for the elastodynamic situation and O = 0 for the acoustic situation. When the Earth's surface is a free surface, is a pseudo‐differential operator for the elastodynamic situation. We introduce its transpose, , and adjoint, , via the following integral relations and respectively. The following properties hold (Kennett et al., 1990; Wapenaar et al., 2004)

For the acoustic situation we simply have .

In state A, the downgoing field in medium A for x at in Table 1 is given by

This time the source is at x , again just above . The receiver is at x at . Note that represents a downgoing transmission response, denoted by , whenever the source and receiver are situated above and below an inhomogeneous slab. Similarly, represents an upgoing transmission response, denoted by (note the minus sign), whenever the source and receiver are situated below and above an inhomogeneous slab. From equation (8), we find

In state A, the upgoing field for x at in Table 1 is zero because medium A is homogeneous below . The downgoing and upgoing fields in state A for x at are defined in a similar way as in state B.

Now that we have discussed all quantities in Table 1, we substitute them into equation (1). Despite the different media (medium A in state A and medium B in state B), this is justified, because between and these media are the same in both states (see Figure 2). Here and in the remainder of this paper, the operator is the same in both states (zero and thus obeying equation (15) when the Earth's surface is considered transparent, or nonzero and obeying equations (15) and (16) when the Earth's surface is considered a free surface). Using equations (10), (13), (15), and (18), setting m = 0 and n = 1 in equation (1), we obtain for x and x just above , see Figure 3.

Figure 3

Visualization of the first and second term in the representation of equation (19).

Next, we derive a representation for in equation (19). Substituting the quantities of Table 2 into equation (1), using equation (10) and setting m = 1 and n = 2, gives for x just above and just above . Because is transparent (i.e., it does not coincide with an interface), equation (20) does not alter if we take at instead of just above it. Thus, taking at , substituting equation (20) into equation (19) (with x in equation (19) replaced by ), we obtain for x and x just above . This is the sought representation for . In a similar way we find the following representation for or, upon substitution of equation (21), for x and x just above . The first term on the right‐hand side is the reflection response of the overburden (Figure 2, medium A [=a]). The second and third terms on the right‐hand side contain the reflection responses of the target zone and the underburden, respectively (media b and c in Figure 2). These terms are visualized in Figure 4.

Table 2

Quantities to Derive a Representation for

	State A:	State B:
	Medium b	Medium B
	Source at x^′ just above S_1	Source at x S just above S_0
S_1	pA+(x,ω)→Iδ(x_H−xH′)	pB+(x,ω)→GB+,+(x,x_S,ω)
	pA−(x,ω)→Rb∪(x,x^′,ω)	pB−(x,ω)→GB−,+(x,x_S,ω)
S_2	pA+(x,ω)→Tb+(x,x^′,ω)	pB+(x,ω)→TB+(x,x_S,ω)
	pA−(x,ω)→O	pB−(x,ω)→O

Figure 4

Visualization of the second and third term in the representation of equation (23).

Note that, if the subsurface would be divided into more and thinner units, the recursive derivation process could be continued, leading to additional terms on the right‐hand side of equation (23). In the limiting case (for infinitesimally thin units), the reflection responses under the integrals could be replaced by local reflection operators, the Green's functions G +,+ by transmission responses T +, and the sum in the right‐hand side would become an integral along the depth coordinate. The resulting expression would be the so‐called generalized primary representation (Fishman et al., 1987; Haines & de Hoop, 1996; Hubral et al., 1980; Kennett, 1974; Resnick et al., 1986; Wapenaar, 1996).

The representation of equation (23) is not meant as a recipe for numerical modelling. However, it is a suited starting point for the derivation of a scheme for target replacement. In equation (23), represents the reflection response from above of the target zone (unit b in Figure 1). Let denote the reflection response of a changed target zone (which we denote as unit ). The reflection response of the entire medium, with the changed target zone, is given by the following representation:

Note that, although it is assumed that the overburden and underburden are unchanged, all quantities on the right‐hand side that contain a propagation path through the target zone are influenced by the changes, which is indicated by the bars. In the following two sections, we discuss the target replacement in detail. First, in section 3 we discuss the removal of the target zone response from the original reflection response . Next, in section 4 we discuss how to insert the response of the changed target into the new reflection response .

Removing the Target Zone From the Original Reflection Response

Given the reflection response of the entire medium, , our aim is to resolve the responses of the media A (=a) and c (i.e., the overburden and underburden, Figure 5). If contained only primary P wave reflections, we could apply simple time‐windowing in the time domain to separate the reflection responses of the different units. However, because of multiple scattering (possibly including surface‐related multiples) and wave conversion, the responses of the different units overlap and cannot be straightforwardly separated by time‐windowing. Here we show that so‐called focusing functions, recently introduced for Marchenko imaging (Slob et al., 2014; Wapenaar et al., 2013), can be used to obtain the responses of media A (=a) and c.

Figure 5

Left: overburden and underburden responses, obtained from the reflection response , using the Marchenko method. Right: modelled responses of the new target zone, to be inserted between the overburden and underburden responses.

We start by defining the focusing function in medium A, with or without free surface just above (Figure 6). Here defines a focal point at boundary , that is, the lower boundary of unit a. Hence, , with x 3,1 denoting the depth of . The coordinate x is a variable in medium A. The superscript + refers to the propagation direction at x (which is downgoing in this case). The focusing function is emitted from all x at into medium A. Due to scattering in the inhomogeneous medium it gives rise to an upgoing function . The focusing conditions for x at can be formulated as with . Equation (25) defines the convergence of to the focal point at , whereas equation (26) states that the focusing function contains no upward scattered components at , because for medium A the half‐space below this boundary is homogeneous. In practical situations evanescent waves are neglected to avoid instability of the focusing function; hence, the delta function in equation (25) should be interpreted as a band‐limited spatial impulse.

Figure 6

Focusing functions and in medium A. The rays stand for the full focusing functions, including all orders of multiple scattering and, in the elastodynamic case, mode conversion.

The focusing functions and for x at and at can be obtained from the reflection response for x just above , using the Marchenko method. We only outline the main features. In Appendix A1, the following relations between , , and are derived and (with x just above and at ) for the situation that the Earth's surface is transparent. For the acoustic case, these equations can be solved for and using the multidimensional Marchenko method (van der Neut et al., 2015; Ravasi et al., 2016; Slob et al., 2014; Wapenaar et al., 2014a). The main assumption is that, in addition to , an estimate of the direct arrival of is available. This can be defined in a smooth model of the overburden. The Marchenko method uses causality arguments to separate the Green's functions from the focusing functions in the left‐hand sides of the time‐domain versions of equations (27) and (28). The multidimensional Marchenko method also holds for the elastodynamic case, except that in this case, an estimate of the direct arrival plus the forward scattered events of needs to be available (Wapenaar & Slob, 2014).

For the situation that the Earth's surface is a free surface, equations (27) and (28) have been modified by Singh et al. (2017), Slob and Wapenaar (2017), and Ravasi (2017), to account for the surface‐related multiple reflections. In these approaches, the surface‐related multiples are present in the reflection response, but not in the focusing functions. For the target replacement procedure discussed in this paper it is more convenient to use focusing functions that include surface‐related multiples. From the derivation in Appendix A1 it follows that for this situation, equation (27) remains valid (but with all quantities now including the surface‐related multiples) and that equation (28) needs to be replaced by (with x just above and at ). The set of equations (27) and (29) for the situation with free surface can be solved in a similar way as the set of equations (27) and (28) for the situation without free surface. A further discussion of the multidimensional Marchenko method to resolve from the reflection response is beyond the scope of this paper.

Assuming the focusing functions and have been found, we use these to resolve the responses of medium A. In Appendix A2, we show that the response to focusing function , when emitted from into medium A, can be quantified as follows for and at , and for x just above and at . Equation (30) describes the transmission response of medium A to the focusing function. The response at is a (band‐limited) spatial impulse (consistent with the focusing condition of equation (25)). Equation (31) describes the reflection response of medium A to the focusing function. The response at is the upgoing part of the focusing function. Both equations (30) and (31) hold for the situation with or without free surface just above . Inverting these equations yields the transmission response (which, according to equation (30), is the inverse of the focusing function ) and the reflection response of medium A, the overburden (Figure 5).

To derive the response of medium A from below, we introduce a second focusing function in medium A, with or without free surface just above (Figure 6). This time defines a focal point at boundary , that is, the upper boundary of unit a. Hence, , with x 3,0 denoting the depth of . The coordinate x is a variable in medium A. The superscript − refers to the propagation direction at x (which is upgoing in this case). The focusing function is emitted from all x at into medium A. Due to scattering in the inhomogeneous medium, it gives rise to a downgoing function . The focusing conditions for x at can be formulated as

Equation (32) defines the convergence of to the focal point at , whereas equation (33) accounts for the downward reflection of the upgoing focusing function at . This term vanishes when the Earth's surface is transparent. In Appendix A3, we show that the response to focusing function , when emitted from into medium A, can be quantified as follows: for and at , and for just below and at . Inverting these equations yields the transmission response (which, according to equation (34), is the inverse of the focusing function ) and the reflection response of medium A from below (Figure 5). In Appendix A4 we show that the focusing functions and are related to the focusing functions and , according to and (with at and at ) for the situation that the Earth's surface is transparent. For the situation that the Earth's surface is a free surface, equation (36) remains valid, and equation (37) needs to be replaced by

(with at and at ).

Next we discuss how to obtain the response of unit c, the underburden, from . We consider again equations (27) and (28) (or (29)), this time with at and replaced by . The focusing functions in medium B can be obtained from the reflection response , using the multidimensional Marchenko method, under the same assumptions as outlined above. Once these focusing functions have been found, they can be substituted into the modified equations (27) and (28) (or (29)), yielding the Green's functions , with x just above and at . Analogous to equation (20), these Green's functions are mutually related via

Inversion of equation (39) yields the reflection response for x and at (Figure 5).

We summarize the steps discussed in this section. Starting with the reflection response of the entire medium, , use the Marchenko method to derive the focusing functions and for medium A. Resolve the responses of the overburden, , , and , by inverting equations (30), (31), (34), and (35). Next, use the Marchenko method to derive the Green's functions , for at . Resolve the reflection response of the underburden, , by inverting equation (39). The resolved responses are free of an imprint of unit b, the target zone.

Inserting a New Target Zone Into the Reflection Response

Given the retrieved responses of the overburden (medium A) and underburden (unit c) and a model of the changed target zone (unit ), our aim is to obtain the reflection response of the entire medium with the new target zone (medium ). The procedure starts by numerically modelling the reflection and transmission responses of the new target zone, and (Figure 5). Next, the response is built up step‐by‐step, using equation (24) as the underlying representation. Analogous to equations (21) and (22), we rewrite equation (24) as a cascade of two representations, as follows: followed by for x and x just above . Quantities in these representations that still need to be determined are , , and .

In Appendix B1, we derive the following equation for the unknown with for x just above , and x and at . Since , , and are known, can be resolved by inverting equation (42). Substituting this into equation (40), together with the other quantities that are already known, yields .

Similarly can be resolved by inverting with for x just above , and x and at . This requires expressions for and .

In Appendix B2 we derive the following representation for for x just above and at . Note that , needed in equation (41), follows by applying equation (18).

In Appendix B3, we derive the following equation for the unknown (with x just above and at ) for the situation that the Earth's surface is transparent. For the situation that the Earth's surface is a free surface, this equation needs to be replaced by (with x just above and at ). Since and are known, can be resolved by inverting either equation (47) or (48).

We summarize the steps discussed in this section. Starting with a model of the new target zone, determine its responses and by numerical modelling. Next, resolve the Green's function of medium , , by inverting equation (42). Substitute this, together with , into equation (40), which yields the reflection response of medium , . Resolve by inverting equation (47) or (48). Substitute this into equation (45) and, subsequently, substitute the result into equation (44). Resolve by inverting equation (44). Substitute this, together with the other quantities that are already known, into equation (41), which yields the sought reflection response .

Numerical Examples

We illustrate the proposed method with two numerical examples. Although the method holds for vertically and laterally inhomogeneous media, for simplicity we consider laterally invariant media in the following examples.

In the first example, we consider the acoustic plane‐wave response of a horizontally layered medium, without free surface (which is the situation after surface‐related multiple elimination). Figure 7 shows the horizontally layered medium. The velocities are given in m/s, the mass densities in kg/m3, and the depth of the interfaces (denoted by the solid lines) in m. To emphasize internal multiples, the mass densities have the same numerical values as the propagation velocities. The layer between 1,200 and 1,400 m represents a reservoir (hence, this is the layer in which changes will take place). The target zone (unit b) includes this reservoir layer (the remainder of the target zone will, however, not undergo any changes). Figure 8a shows the numerically modelled plane‐wave reflection response at in the time domain, convolved with a Ricker wavelet with a central frequency of 50 Hz (note that we replaced the boldface symbol R by a plain R, because the acoustic response is a scalar function; moreover, we replaced ω by t because the response is shown in the time domain). The reflections from the top and bottom of the reservoir are indicated by arrows. We consider a time‐lapse scenario, in which the velocity in the reservoir is changed from 3,000 to 2,500 m/s (and a similar change is applied to the mass density). Figure 8b shows the numerically modelled time‐lapse reflection response , and Figure 8c shows the difference . Note the significant multiple coda, following the difference response of the reservoir. Our aim is to show that the time‐lapse response (Figure 8b) can be predicted from the original response (Figure 8a) by target replacement.

Figure 7

Horizontally layered medium for the plane‐wave experiment, with the three units indicated. The Earth's surface is considered transparent.

Figure 8

(a) Numerically modelled reflection response of the model of Figure 7. (b) Numerically modelled time‐lapse response. (c) The difference of the responses in (a) and (b).

Following the procedure discussed in section 3 (simplified for the 1‐D situation), we remove the response of the target zone from the reflection response . The overburden response , resolved from equation (31), is shown in the time domain in Figure 9a. Note that it contains the first two events of and a coda due to the internal multiples in the low‐velocity layer in the overburden. The underburden response , resolved from equation (39), is shown in Figure 9b. For display purposes it has been shifted in time, so that the travel times correspond with those in Figure 8a.

Figure 9

(a) The response of medium A (the overburden), retrieved from . (b) The response of unit c (the underburden), retrieved from . (c) Numerically modelled response of unit (the new target zone).

Following the procedure discussed in section 4 (simplified for the 1‐D situation), we predict the time‐lapse response. To this end, we first model the response of the new target zone, . This is shown in Figure 9c. For display purposes, it has been shifted in time so that the travel time to the top of the reservoir corresponds with that in Figure 8a. The predicted time‐lapse reflection response at the surface, , obtained with the representations of equations (40) and (41), is shown in the time domain in Figure 10a. The numerically modelled response of Figure 8b, is once more shown (as a reference) in Figure 10b. The difference of the predicted and modelled responses is shown in Figure 10c and appears to be practically zero. This confirms that the new reflection response has been very accurately predicted by the proposed method.

Figure 10

(a) The predicted time‐lapse response , constructed from the responses in Figure 9. (b) For comparison, the numerically modelled time‐lapse response. (c) The difference of the responses in (a) and (b).

For the next example, we consider a 2‐D acoustic point‐source response of a horizontally layered medium. The medium is shown in Figure 11. Note that the overburden and underburden contain more layers than in the previous example. Figure 12a shows the numerically modelled response at the surface in the time domain, for a fixed source at x =(0,0) and variable receivers at x =(x 1,,0). Because the medium is horizontally layered, the responses to sources at other positions at are simply laterally shifted versions of the response in Figure 12a. In the time‐lapse scenario, the velocity in the reservoir layer is changed from 3,000 to 2,500 m/s (and a similar change is applied to the mass density). Figure 12b shows the difference of the numerically modelled responses and . The responses in this and the following figures are displayed with a small time‐dependent gain of to emphasize the internal multiples.

Figure 11

Horizontally layered medium for the 2‐D experiment.

Figure 12

(a) Numerically modelled 2‐D reflection response. (b) Numerically modelled difference response.

We use our standard implementation of the Marchenko method (Thorbecke et al., 2017) for the estimation of the focusing functions. Next, because the medium is horizontally layered, we efficiently carry out the layer replacement method in the wave number‐frequency domain (hence, all integrals from equation (30) onward reduce to straightforward products of the transformed quantities). Figure 13a shows the overburden response , resolved from equation (31) in the wave number‐frequency domain and transformed back to the space‐time domain. Note that the internal multiples of the overburden, indicated by the arrows, have been recovered from behind the reflection response of the reservoir layer. The modelled response of the new target zone, at , is shown in Figure 13b, for a fixed source at x = (0,1400) m and variable receivers at  m. The predicted time‐lapse reflection response at the surface, of the overburden and target zone, , obtained with the representation of equation (40) in the wave number‐frequency domain, is shown in Figure 14a. The numerically modelled time‐lapse response is shown (as a reference) in Figure 14b. Next, the response of the underburden is included, using the representation of equation (41) in the wave number‐frequency domain. This yields the predicted time‐lapse reflection response at the surface of the entire medium, , see Figure 15a. The numerically modelled time‐lapse response of the entire medium is shown in Figure 15b. Although the match is not as perfect as in the 1‐D example (Figure 10c), Figure 15 shows that the 2‐D time‐lapse response has been accurately predicted. We used dip filtering to suppress artifacts related to the finite aperture and the negligence of evanescent waves. This explains the diminishing amplitudes of the early reflections at large offsets.

Figure 13

(a) The response of medium A (the overburden), retrieved from . (b) Numerically modelled response of the new target zone.

Figure 14

(a) The predicted time‐lapse response , constructed from the responses in Figure 13. (b) For comparison, the numerically modelled time‐lapse response.

Figure 15

(a) The predicted time‐lapse response , constructed from and the response of the underburden. (b) For comparison, the numerically modelled time‐lapse response.

Discussion

The numerical examples in the previous section show that under ideal circumstances the proposed method accurately predicts the time‐lapse responses. Hence, these examples validate the theory. In practice, there will be several factors that limit the accuracy. First, the direct arrivals of the focusing function , needed to initiate the Marchenko scheme, are in practice defined in estimated models of the medium. Hence, the amplitudes and travel times of these direct arrivals will not be exact. The Marchenko method is robust to small‐to‐moderate errors in the direct arrival, in the sense that it predicts the multiples in the focusing functions and Green's functions, but these predicted multiples will exhibit similar amplitude and travel time errors as the direct arrival (Broggini et al., 2014; Wapenaar et al., 2014b). The errors in and largely compensate each other in the inversion of equation (31), to obtain the overburden response . Hence, will be retrieved very accurately, despite the errors in the direct arrival (it has been previously observed that the Marchenko method for obtaining data at the surface is very robust; Meles et al., 2016; van der Neut & Wapenaar, 2016). This implies that multiples generated in the overburden are accurately separated from the response of deeper layers. The response of the overburden from below, , is obtained by inverting equation (35). Here the amplitude errors in and largely compensate each other, but travel time errors will result in an overall time shift of . A similar remark holds for the underburden response . These errors will propagate into the predicted time‐lapse response. We expect that the errors in the predicted primaries and low‐order multiples will be of the same order as the errors in the direct arrivals and that these errors will grow for higher‐order multiples.

The accuracy of the predicted time‐lapse response will further be limited by losses in the medium, inaccuracies in the deconvolution for the source wavelet, the finite length of the acquisition aperture, and incomplete sampling (particularly for 3‐D applications). Currently much research is going on to improve the Marchenko method to address these issues (Ravasi et al., 2016; Slob, 2016; Staring et al., 2017; van der Neut & Wapenaar, 2016). The proposed target replacement scheme will benefit from these developments.

The computational costs of the proposed method depend on the implementation. For the numerical examples in the previous section we took advantage of the fact that the medium is horizontally layered. We implemented the 2‐D layer replacement in the wave number‐frequency domain. This implies that the inversion of the various integral equations is replaced by a straightforward scalar inversion per wave number‐frequency combination. For laterally varying media, the integral equations should be solved in the space‐frequency domain. After discretization, this comes to a matrix inversion for each frequency component. In several cases (equations (42) and (44)) the matrix inversion can efficiently be replaced by a series expansion, which can be terminated after a few terms, depending on the number of multiples that need to be taken into account. All at all, removing the target zone (section 3) requires applying the Marchenko method at two depth levels ( and ) and five matrix inversions (per frequency component) to solve integral equations (30), (31), (34), (35), and (39). Inserting the new target zone (section 4) requires numerical modelling of the target zone response and three matrix inversions (per frequency component) to solve integral equations (42), (44), and (47). The costs for substituting the results into equations (40) and (41) are negligible in comparison with the matrix inversions. Despite the significant number of steps for the entire process, the total costs should be seen in perspective with other methods. In comparison with numerically modelling the entire time‐lapse reflection response, our method requires numerical modelling of the target zone response only. The additional costs for the Marchenko method and the matrix inversions are significant but not excessive. For example, applying the Marchenko method at two depth levels is feasible, considering the fact that some Marchenko imaging methods apply this method for a large range of depth levels in an image volume (Broggini et al., 2014; Behura et al., 2014). The trade‐off between the cost reduction for the numerical modelling and the cost increase related to the Marchenko method and the matrix inversions depends on the implementation details and needs further investigation.

Conclusions

We have proposed an efficient two‐step process to replace the response of a target zone in a reflection response at the Earth's surface. In the first step, the response of the original target zone is removed from the reflection response, using the Marchenko method. In the second step, the modelled response of a new target zone is inserted between the overburden and underburden responses. The method holds for vertically and laterally inhomogeneous lossless media. It fully accounts for all orders of multiple scattering and, in the elastodynamic case, for wave conversion. It can be employed to predict the time‐lapse reflection response for a range of target‐zone scenarios. For this purpose, the first step needs to be carried out only once. Only the second step needs to be repeated for each target‐zone model. Since the target zone covers only a small part of the entire medium, repeated modelling of the target‐zone response (and inserting it each time between the same overburden and underburden responses) is a much more efficient process than repeated modelling of the entire reflection response, but there are also additional costs related to the Marchenko method and several matrix inversions. This method may find applications in time‐lapse full wave form inversion, for example, to monitor fluid flow in an aquifer, subsurface storage of waste products, or production of a hydrocarbon reservoir. Since all multiples are taken into account, the coda following the response of the target zone may be employed in the inversion. Because of the high sensitivity of the coda for changes in the medium (Snieder et al., 2002), this may ultimately improve the resolution of the inverted time‐lapse changes. Finally, when medium changes are not restricted to a reservoir, the target zone should be taken sufficiently large to include those parts of the embedding medium in which changes take place. This will of course have a limiting effect on the efficiency gain.

Supporting Information S1

Click here for additional data file.

Data Set S1

Click here for additional data file.

Data Set S2

Click here for additional data file.

Data Set S3

Click here for additional data file.

Data Set S4

Click here for additional data file.

Table A1

Quantities to Derive Marchenko Representations

	State A:	State B:
	Medium C	Medium A
	Source at x R just above S_0	Focus at x^′ at S_1
S_0	pA+(x,ω)→Iδ(x_H−x_H,R)	pB+(x,ω)→F1,A+(x,x^′,ω)
	+r^∩RC∪(x,x_R,ω)	+r^∩F1,A−(x,x^′,ω)
	pA−(x,ω)→RC∪(x,x_R,ω)	pB−(x,ω)→F1,A−(x,x^′,ω)
S_1	pA+(x,ω)→GC+,+(x,x_R,ω)	pB+(x,ω)→Iδ(x_H−xH′)
	pA−(x,ω)→GC−,+(x,x_R,ω)	pB−(x,ω)→O

Table A2

Quantities to Derive the Response to

	State A:	State B:
	Medium A	Medium A
	Source at x^′′ just below S_1	Focus at x^′ at S_1
S_0	pA+(x,ω)→r^∩TA−(x,x^′′,ω)	pB+(x,ω)→F1,A+(x,x^′,ω)
		+r^∩F1,A−(x,x^′,ω)
	pA−(x,ω)→TA−(x,x^′′,ω)	pB−(x,ω)→F1,A−(x,x^′,ω)
S_1	pA+(x,ω)→RA∩(x,x^′′,ω)	pB+(x,ω)→Iδ(x_H−xH′)
	pA−(x,ω)→Iδ(x_H−xH′′)	pB−(x,ω)→O

Table A3

Quantities to Derive the Response to

	State A:	State B:
	Medium A	Medium A
	Source at x^′′ just above S_0	Focus at x^′ at S_0
S_0	pA+(x,ω)→Iδ(x_H−xH′′)	pB+(x,ω)→r^∩Iδ(x_H−xH′)
	+r^∩RA∪(x,x^′′,ω)
	pA−(x,ω)→RA∪(x,x^′′,ω)	pB−(x,ω)→Iδ(x_H−xH′)
S_1	pA+(x,ω)→TA+(x,x^′′,ω)	pB+(x,ω)→F2,A+(x,x^′,ω)
	pA−(x,ω)→O	pB−(x,ω)→F2,A−(x,x^′,ω)

Table B1

Quantities to Derive Representation for

	State A:	State B:
	Medium b¯	Medium B¯
	Source at x^′′ just below S_2	Source at x S just above S_0
S_1	pA+(x,ω)→O	pB+(x,ω)→G¯B+,+(x,x_S,ω)
	pA−(x,ω)→T¯b−(x,x^′′,ω)	pB−(x,ω)→G¯B−,+(x,x_S,ω)
S_2	pA+(x,ω)→R¯b∩(x,x^′′,ω)	pB+(x,ω)→T¯B+(x,x_S,ω)
	pA−(x,ω)→Iδ(x_H−xH′′)	pB−(x,ω)→O

Table B2

Quantities to Derive Equation for

	State A:	State B:
	Medium B¯	Medium B¯
	Source at x S just above S_0	Source at x^′ just below S_2
S_0	pA+(x,ω)→Iδ(x_H−x_H,S)	pB+(x,ω)→r^∩T¯B−(x,x^′,ω)
	+r^∩R¯B∪(x,x_S,ω)
	pA−(x,ω)→R¯B∪(x,x_S,ω)	pB−(x,ω)→T¯B−(x,x^′,ω)
S_2	pA+(x,ω)→T¯B+(x,x_S,ω)	pB+(x,ω)→R¯B∩(x,x^′,ω)
	pA−(x,ω)→O	pB−(x,ω)→Iδ(x_H−xH′)