Shohei Minato∗and Ranajit Ghose, Delft University of Technology SUMMARY
Hydraulic properties of a fractured reservoir are often con-trolled by large fractures. In order to seismically detect and characterize them, a high-resolution imaging method is nec-essary. We apply a non-linear imaging condition to image fractures, considered as non-welded interfaces. We derive the imaging condition from the general correlation-type represen-tation theorem, assuming compliances to be real-valued func-tions. We investigate the P-wave image due to P-wave sources and the effect of source illumination. We present here numer-ical modeling results for (1) a single dry fracture, (2) orthogo-nally intersecting fractures, and (3) deviated multiple fractures in a multi-layered subsurface. Our results show that the non-linear terms in the imaging condition help to cancel the arte-facts appearing in case of conventional imaging condition, and improve the resolution of the final image. For a multi-layered model, the one-sided source illumination from the earth’s sur-face can image the multiple fractures. However, an incomplete source illumination restricts the removal of all the artefacts. In this case, installing additional borehole sources can greatly improve the resolution of the fracture distribution.
INTRODUCTION
Seismic detection and characterization of large fractures are important because they dominate the hydraulic properties of a fractured reservoir. When fractures are large in size compared to the seismic wavelength, they generate scattered waves due to an incident seismic wave and the effective medium theory cannot be used. Willis et al. (2006) used scattered waves to obtain the spatial orientation and distribution of multiple frac-tures. Recently, Minato and Ghose (2013, 2014) have shown that it is possible to obtain further details of the individual frac-tures, and proposed an approach to characterize the spatially heterogeneous elastic compliance along the fracture plane through wavefield extrapolation. The spatial heterogeneity along the fracture is a key determinant for fracture-associated hydraulic properties. The proposed approach requires precise informa-tion of the posiinforma-tion of fractures, and hence asks for a method to image the fractures with high resolution. In this vein, Mi-nato and Ghose (2014) used the reverse-time migration with a conventional imaging condition.
Recently, a non-linear imaging condition has been proposed for high-resolution imaging of the boundary of contrasting acous-tic or elasacous-tic constants (Fleury and Vasconcelos, 2012; Ravasi and Curtis, 2013). In this case, the imaging condition is a non-linear function of the perturbed wavefield, where the perturbed wavefield is defined as the difference between total wavefield and background wavefield. This new imaging condition al-lows achieving higher resolution in the migrated image than with the conventional imaging condition which is linear with respect to the perturbed wavefield (Ravasi and Curtis, 2013).
The non-linear imaging condition is derived using correlation-type reciprocity theorem applied to the perturbed wavefield. Then the migration image can be seen as offset and zero-time scattered amplitudes.
Unlike a boundary between contrasting elastic or acoustic prop-erties, a fracture is often represented by a linear-slip bound-ary across which the stress is continuous but the displacement is discontinuous (e.g., Schoenberg, 1980; Pyrak-Nolte et al., 1990). This interface is generalized as a non-welded interface (e.g., Wapenaar et al., 2004). By appropriately including the effect of the non-welded interface into the reciprocity theo-rem, a general correlation-type representation theorem can be derived (Wapenaar, 2007). This leads to the possibility of de-riving a non-linear imaging condition when the medium con-tains non-welded interfaces.
In this study, we apply the non-linear imaging condition to im-age fractures, considered as non-welded interfaces. We first derive the non-linear imaging condition from the general cor-relation type representation theorem. Next, we show numeri-cally the importance of the non-linear imaging condition and the effect of source illumination. We present 3 examples: (1) a single fracture, (2) orthogonally intersecting fractures, and (3) deviated multiple fractures in a multi-layered subsurface.
THEORY
The non-linear imaging condition is derived from the correlation-type representation theorem which relates the wavefield be-tween two different states (states A and B). We use the general representation theorem (Wapenaar, 2007) which includes the effect of a non-welded interface.
The boundary condition for a non-welded interface can be writ-ten in the general form as,
[Mu] =− jωY < Mu >, (1) where [·] and < · > indicate the jump and the average, respec-tively, across the non-welded interface. For elastic wave, we have uT= ( vT,−τττT1,−τττT2,−τττT3), where v andτττiare the
par-ticle velocity and the traction vectors. The matrices in equation 1 are defined as,
M = ( I 0 0 0 0 nb 1I nb2I nb3I ) and Y = ( 0 Sb ρρρb 0 ) , (2) where the superscript b denotes the boundary parameters. That is, nbi, Sb andρρρbare the vector normal to the interface, the boundary compliance tensor and the boundary density tensor, respectively.
We consider the general correlation-type representation theo-rem (equation 59 in Wapenaar, 2007) with the two position x′ and x′′inside a closed surface∂D (Figure 1). When we con-sider a horizontal fracture and a rotationally invariant
Fracture imaging
ance matrix (Schoenberg, 1980), the boundary parameters re-duce to (nb)T= (0, 0, 1)T, Sb= diag(ηT,ηT,ηN) andρρρb= 0,
withηT andηN being the tangential and the normal
compli-ance, respectively. In this case, the general correlation-type representation theorem yields to:
¯ v∗i, j(x′, x′′) + vi, j(x′, x′′) = − I ∂D [ ¯ v∗m,i(x, x′)τmn, j(x, x′′) + vm, j(x, x′′) ¯τmn,i∗ (x, x′) ] nnd2x + jω ∫ ∂Dint [ ∆ηTτ¯13,i∗ (x, x′)τ13, j(x, x′′) +∆ηTτ¯23,i∗ (x, x′)× τ23, j(x, x′′) +∆ηNτ¯33,i∗ (x, x′)τ33, j(x, x′′) ] d2x, (3)
where we assume that the two states have identical elastic con-stants but different fracture compliances. ¯vi, j(x′, x′′) and vi, j(x′, x′′)
indicate i−th component of particle velocity at x′due to j−th direction of point force at x′′ in state A and state B, respec-tively. ¯τpq,i and τpq, j are their corresponding stress values.
∆ηT and ∆ηN are the perturbation in compliances between
state A and state B. We assume the compliances to be real-valued functions. The second integral in equation 3 is taken along (any number of) non-welded interfaces defined as∂Dint.
Following the procedure similar to Vasconcelos et al. (2009), we assume the total wavefield to be vi, j= v0i, j+ vSi, j, where v0i, j
and vSi, jare the background wavefield and the perturbed wave-field due to the non-zero compliance perturbation, respectively. When we consider the two states in equation 3 to be both per-turbed fields,∆ηTand∆ηNvanish, and we have:
vSi, j∗(x′, x′′) + vSi, j(x′, x′′) = − I ∂D [ v0m,i∗(x, x′)τmn, jS (x, x′′) + v0m, j(x, x′′)τmn,iS∗ (x, x′) ] nnd2x − I ∂D [ vSm,i∗(x, x′)τmn, j0 (x, x′′) + vSm, j(x, x′′)τmn,i0∗ (x, x′) ] nnd2x − I ∂D [ vSm,i∗(x, x′)τmn, jS (x, x′′) + vSm, j(x, x′′)τmn,iS∗ (x, x′)]nnd2x. (4) This equation is the basis of the non-linear imaging condition that we propose here. It indicates that the perturbed wave-field at any position inside∂D can be retrieved from the values along the closed surface∂D. The equation 4 is the same as that in Ravasi and Curtis (2013); however, these authors do not consider a non-welded interface. We obtain the same equation due to our assumption of the compliances as real-valued func-tions. Although we have considered horizontal fractures to de-rive equation 3, equation 4 is valid for fractures of arbitrary orientation, because∆ηT and∆ηNalways vanish. We obtain
the imaging condition at imaging point xPas the zero-time and
zero-offset amplitude of the perturbed wavefield (Ravasi and Curtis, 2013): IΨN,ΨI(xP)≈ 2 ρcK ∫ ∞ 0 I ∂D [ G0Ψ∗ K,ΨI(xP, x)G S ΨN,ΨK(xP, x) ] d2xdω + 2 ρcK ∫ ∞ 0 I ∂D [ GSΨ∗ I,ΨK(xP, x)G 0 ΨN,ΨK(xP, x) ] d2xdω + 2 ρcK ∫ ∞ 0 I ∂D [ GSΨ∗ I,ΨK(xP, x)G S ΨN,ΨK(xP, x) ] d2xdω, (5)
where we have replaced the particle velocity v0,Si, j by the po-tential field G0,SΨ
I,ΨN, with source and receiver represented by potentials (I, N, K take the value of 0,1,2 or 3 with 0 being the P-wave potential). Note that, equation 5 is derived already af-ter the far-field approximation.
When a P-wave source is activated and no S-wave source is available as in case of most marine seismic experiments, ig-noring the terms including the S-wave sources in equation 5 gives the approximate imaging condition for the P-wave im-age due to only P-wave sources:
IPP(xP)≈ 2 ρcP ∫ ∞ 0 I ∂D [ G0Ψ0∗,Ψ0(xP, x)GSΨ0,Ψ0(xP, x) ] d2xdω + 2 ρcP ∫ ∞ 0 I ∂D [ GSΨ0∗,Ψ0(xP, x)G0Ψ0,Ψ0(xP, x) ] d2xdω + 2 ρcP ∫ ∞ 0 I ∂D [ GSΨ0∗,Ψ0(xP, x)GSΨ0,Ψ0(xP, x) ] d2xdω. (6) One can see that the first and the second integrals in equation 6 contain the crosscorrelation between the background wavefield and the perturbed wavefield, which indicates that they are lin-ear with respect to the perturbed wavefield. However, the third integral contains the autocorrelation of the perturbed wave-field, which suggests that it is non-linear with respect to the perturbed wavefield. Note that, in real applications, the back-ground wavefield G0Ψ
I,ΨNand the perturbed wavefield G S ΨI,ΨN on the right hand side of equation 6 may be estimated by cal-culating the wavefield using a macro-velocity model and by extrapolating the receiver responses, respectively.
3 1
n n
b
Figure 1: Configuration of general correlation-type represen-tation theorem.∂Dintindicates a non-welded interface. n and
nbare vectors normal to∂D and ∂Dint, respectively.
NUMERICAL MODELING
We test the non-linear imaging condition for imaging fractures using 2D numerical modeling. In this study, we particularly focus on the source illumination and on the contribution of the non-linear term in the imaging condition to resolution. There-fore, we assume that both background and perturbed wave-fields are already estimated and we ignore the discussion on extrapolation of receiver responses.
Single dry fracture
The first example shows that a 200-m long vertical fracture
is located inside a closed surface∂D of radius 400 m (Fig-ure 2a). 40 transient, point P-wave sources are installed along ∂D. The values of normal and tangential compliances are taken from earlier field experiments (Worthington and Hudson, 2000). We assume a dry fracture, i.e.,ηN=ηT= 1.1× 10−9
(m/Pa). The homogeneous background represents sandstone with VP= 4000 (m/s), VS= 2350 (m/s) andρ = 2500 (kg/m3).
We calculate the seismic response with 100 Hz Ricker wavelet using the staggered-grid finite difference time domain (FDTD) method (Coates and Schoenberg, 1995), which implicitly eval-uates the linear-slip boundary condition.
Figures 2(b) and (c) show a snapshot of the background P-wave potential field G0Ψ0,Ψ0(x, xS1) with the first source at xS1
and the perturbed field GSΨ0,Ψ0(x, xS1), respectively. The total
wavefield can be obtained from G0Ψ0,Ψ0+ GSΨ0,Ψ0. Figure 2(c) shows the transmission and the reflection due to the fracture in the perturbed wavefield. The conventional imaging condi-tion is the crosscorrelacondi-tion of the background wave field and the perturbed wavefield (first and second term of equation 6). The result of using the conventional imaging condition (us-ing 1 source) is shown in Figure 3(a). The negative ampli-tude region on the left side of the fracture (blue area in Fig-ure 3a) comes from the correlation between the incident wave and the perturbed transmitted wave. The positive amplitudes in Figure 3a have peaks along the right side of the fracture (red color), which is a result of correlation between the inci-dent wave and the reflected wave. The non-linear term in the imaging condition is due to the autocorrelation of the perturbed wavefield (third term of equation 6) and is illustrated in Fig-ure 3(b). The summation of these two panels is the result of us-ing the non-linear imagus-ing condition (Figure 3c). One can see that the non-linear imaging condition shows more clearly the peak positive amplitudes along the fracture than the conven-tional imaging condition. However, they are spatially shifted from the true fracture position (dashed line in Figure 3c) and artefacts appear on both right and left sides of the fracture. These artefacts are reduced by increasing the source illumina-tion. Figure 4(a) is the result of conventional imaging condi-tion using all sources in Figure 2(a). The peak negative ampli-tudes are correctly located along the fracture position. How-ever, the artefact due to the correlation of the incident wave with the transmitted wave still remains on both right and left sides of the fracture (Figure 4a). Stacking the non-linear term (Figure 4b) contributes to canceling these artefacts and results in higher resolution (Figure 4c). Comparing the results be-tween Figure 3(c) and Figure 4(c), we can infer that the posi-tive peak amplitudes, which are spatially shifted from the true fracture position using only 1 source (Figure 3c), contribute to the sidelobes in the fracture image (Figure 4c).
Intersecting orthogonal fractures
The next example illustrates the case of two orthogonally in-tersecting fractures. A vertical and a horizontal fracture, with the same fracture compliances as in the previous example, are considered within the closed surface in Figure 2(a). We use the effective medium theory (Nichols et al., 1989) to represent the intersection point of two fractures. Figure 5(a), (b) and (c) show results with conventional imaging condition, using
400 m 200 m 0 100 200 0 100 200 300 400 1(m) 3 (m) 0 100 200 0 100 200 300 400 1(m) (a) (b) (c)
Figure 2: (a) Fracture and source distribution on a vertical plane. (b) Snapshot of the background P-wave potential field
G0Ψ0,Ψ0(x, xS1) at t = 0.128 s in the region enclosed by the
dashed line in (a). (c) Same to (b) but for the perturbed field
GSΨ0,Ψ0(x, xS1). 0 100 200 0 100 200 300 400 1(m) 3 (m) (a) (b) (c) 0 100 200 0 100 200 300 400 1(m) 0 100 200 0 100 200 300 400 1(m)
Figure 3: (a) Conventional imaging condition with only 1 source at xS1. (b) Non-linear term in the non-linear imaging
condition. (c) Summation of (a) and (b). Dashed line shows the horizontal position of the fracture.
0 100 200 0 100 200 300 400 1(m) 3 (m) (a) (b) (c) 0 100 200 0 100 200 300 400 1(m) 0 100 200 0 100 200 300 400 1(m)
Figure 4: (a), (b) and (c) are same as Figure 3(a), (b) and (c), respectively, but with all sources.
0 100 200 0 100 200 300 400 1(m) 3 (m) (a) (b) (c) 0 100 200 0 100 200 300 400 1(m) 0 100 200 0 100 200 300 400 1(m)
Figure 5: (a), (b) and (c) are same as Figure 4(a), (b) and (c), respectively, but for a set of orthogonal fractures.
Fracture imaging
the non-linear term, and using non-linear imaging condition, respectively. As in the case of a single fracture, considering the non-linear term significantly reduces the artefacts in the conventional imaging. Note that the intersection of the two fractures has a larger absolute amplitude than other positions of the fracture. This is because this point is more compliant than others due to intersection of two fractures.
Deviated multiple fractures in a multi-layerd subsurface The final example presents a more practical situation. We con-sider a five-layered earth, inspired by the numerical example of Willis et al. (2006). The fractured reservoir is located in the third layer (Figure 6). The velocity and density of each layer are shown in Table 1. The fractures are deviated by 50 degree from the horizontal axis, with a separation of 100 m between individual fractures. We install 27 sources with 100 Hz Ricker wavelet along the earth’s surface (yellow stars in Figure 6) to simulate a surface-seismic experiment. We consider P-wave sources; we ignore the free surface to focus on the P-wave im-age. Furthermore, we calculate the background wavefield by using the smooth velocity model.
Figure 7(a) shows the result with conventional imaging condi-tion in the region enclosed by the white dashed line in Figure 6. Due to one-sided illumination from the earth’s surface, the artefacts below the fractures (blue area in Figure 7a) distorts the shape and the extent of the fractures. The non-linear term (Figure 7b) contributes to canceling these artefacts, as in the case of the single fracture experiment. The result with the non-linear imaging condition (Figure 7c) shows more clearly the fractures as linear, local amplitude depressions (blue color). Due to one-sided illumination, the sidelobes are clearer at shal-lower parts of the fractures (red color). This result is similar to that for a single fracture with 1 source (Figure 3c). Note that the amplitude spread around x3= 400 m is mainly due to the
interaction between the layer boundary, the fractures and the smooth background model.
We finally check the improvement as a result of further source illumination. We install 15 subsurface sources along a vertical borehole (blue stars in Figure 6). Figure 7(d) shows the result of the proposed non-linear imaging condition using these ad-ditional borehole sources. Due to better illumination from the right, the fractures on the right side in Figure 7(d) are clearer than in Figure 7(c).
CONCLUSION
We propose a non-linear imaging condition for high-resolution imaging of fractures. The fractures are represented by the linear-slip boundary condition, as they are generalized as non-welded interfaces. We show that the imaging condition as in previous studies which do not consider non-welded interfaces can be retrieved from the general correlation-type representa-tion theorem by assuming the compliances to be real-valued functions. We perform 2D numerical modeling to test the re-sults for (1) single dry fracture, (2) orthogonally intersecting fracture and (3) deviated multiple fractures in a multi-layered subsurface. The results show that the non-linear term in the
imaging condition contributes to canceling the artefacts that appear in case the conventional imaging condition is used, and improves the resolution of the final image. In our five-layer model, one-sided source illumination from earth’s surface can image the multiple fractures. Due to one-sided illumination, the sidelobes become more prominent at shallower parts of the fractures. 0 200 400 800 1000 600 0 500 1000 1500 2000 1(m) 3 (m)
Figure 6: Five-layered model. The third layer includes the de-viated multiple fractures (white lines). A surface source array (yellow stars) and a subsurface source array (blue stars) are considered. VP(m/s) VS(m/s) ρ(kg/m3) Layer 1 3000 1756 2200 Layer 2 3500 2060 2250 Layer 3 4000 2353 2300 Layer 4 3500 2060 2250 Layer 5 4000 2353 2300
Table 1: Isotropic velocity and density in each layer. For Layer 3 (fractured medium), the background value is shown.
1(m) 500 1000 1500 (b) (d) 3 (m) 250 500 750 (c) 1(m) 3 (m) 250 500 750 500 1000 1500 (a)
Figure 7: The result using (a) the conventional imaging con-dition, (b) the non-linear term and (c) the non-linear imaging condition, using a surface source array (yellow stars in Fig-ure 6). (d) Same as (c) but using an additional subsurface source array (blue stars in Figure 6).
ACKNOWLEDGMENTS
This work is supported by The Netherlands Research Centre for Integrated Solid Earth Science.
EDITED REFERENCES
Note: This reference list is a copy-edited version of the reference list submitted by the author. Reference lists for the 2014 SEG Technical Program Expanded Abstracts have been copy edited so that references provided with the online metadata for each paper will achieve a high degree of linking to cited sources that appear on the Web.
REFERENCES
