Effects of variations in fluid properties and fracture geometry on dispersion,anisotropy, and reflection in media with planar fractures

Our objective was to discover the effect of variations in fluid properties and fracture geometry on the velocity of seismic wave propagation in fluid‐saturated media with parallel planar fractures. We used numerical models calculated by analytical solutions to examine the behaviour of P‐wave phase velocity dispersion in the normal direction to layering, in non‐porous and porous media with planar fractures. We also examined the anisotropy of low frequency phase and group velocities of fast and slow P‐waves and angular‐dependent reflection coefficients in media with planar fractures, under conditions of saturation by fluids with varying bulk moduli, densities, and fracture apertures. We defined several parametre , , and characterising dispersion, characterising anisotropy, characterising the difference between fast and slow modes, and R₀ and characterising reflection. Our results show that the behaviour of dispersion shows wider stopbands in the case of gas saturation. Concavity indicator of dispersion for gas saturation was greater than that for liquid saturation and is usually greater than one. Anisotropy is more sensitive to bulk modulus contrast than to density contrast between the solid and the fluid, and is more sensitive to density contrast than to bulk modulus contrast. The case of gas saturation usually had a greater negative R₀ and a greater value of compared with those of brine and heavy and light oil saturations. Our results are helpful in distinguishing fluid types saturating geophysical fractures and estimating the aperture and spacing of planar fractures. In seismic exploration, bulk modulus and fluid density can provide useful information in distinguishing among brine, oil, and gas; fracture geometry is important to estimate the permeability of reservoirs.     

Measuring changes in fracture properties from temporal variations in anisotropic attenuation of microseismic waveforms

We investigate fracture‐induced attenuation anisotropy in a cluster of events from a microseismic dataset acquired during hydraulic fracture stimulation. The dataset contains 888 events of magnitude ?3.0 to 0.0. We use a log‐spectral‐amplitude‐ratio method to estimate change in over a half‐hour time period where fluid is being injected and an increase in fracturing from S‐wave splitting analysis has been previously inferred. A Pearson's correlation analysis is used to assess whether or not changes in attenuation with time are statistically significant. P‐waves show no systematic change in during this time. In contrast, S‐waves polarised perpendicular to the fractures show a clear and statistically significant increase with time, whereas S‐waves polarised parallel to the fractures show a weak negative trend. We also compare between the two S‐waves, finding an increase in with time. A poroelastic rock physics model of fracture‐induced attenuation anisotropy is used to interpret the results. This model suggests that the observed changes in t* are related to an increase in fracture density of up to 0.04. This is much higher than previous estimates of 0.025 ± 0.002 based on S‐wave velocity anisotropy, but there is considerably more scatter in the attenuation measurements. This could be due to the added sensitivity of attenuation measurement to non‐aligned fractures, fracture shape, and fluid properties. Nevertheless, this pilot study shows that attenuation measurements are sensitive to fracture properties such as fracture density and aspect ratio.     

Experimental quantification of the seismoelectric transfer function and its dependence on conductivity and saturation in loose sand

Under certain circumstances, seismic propagation within porous media may be associated to the conversion of mechanical energy to electromagnetic energy, which is known as a seismo‐electromagnetic phenomenon. The propagation of fast compressional P‐waves is more specifically associated to the manifestations of a seismoelectric field linked to the fluid flows within the pores. The analysis of seismoelectric phenomena, which requires the combination of the theory of electrokinetics and Biot's theory of poroelasticity, provides us with transfer function that links the coseismic seismoelectric field E to the seismic acceleration . To measure the transfer function, we have developed an experimental setup enabling seismoelectric laboratory observation in unconsolidated quartz sand within the kilohertz range. The investigation focused on the impact of fluid conductivity and water saturation over the coseismic seismoelectric field. During the experiment, special attention was given to the accuracy of electric field measurements. We concluded that, to obtain a reliable estimate of the electric field amplitude, the dipole from which the potential differences are measured should be of much smaller length than the wavelength of the propagating seismic field. Time‐lapse monitoring of the seismic velocities and seismoelectric transfer functions were performed during imbibition and drainage experiments. In all cases, the quantitative analysis of the seismoelectric transfer function was in good agreement with theoretical predictions. While investigating saturation variations from full to residual water saturation, we showed that the ratio undergoes a switch in polarity at a particular saturation , which also implies a sign change of the filtration, traducing a reversal of the relative fluid displacement with respect to the frame. This sign change at critical saturation stresses a particular behaviour of the poroelastic medium: the dropping of the coseismic electric field to zero traduces the absence of relative pore/fluid displacements representative of a Biot dynamically compatible medium. We concluded from our experimental study in loose sand that the measurements of the coseismic seismoelectric coupling may provide information on fluid distribution within the pores and that the reversal of the seismoelectric field may be used as an indicator of the dynamically compatible state of the medium.     

Characterizing the degree of amplitude‐variation‐with‐offset nonlinearity in seismic physical modelling reflection data

The nonlinearity of the seismic amplitude‐variation‐with‐offset response is investigated with physical modelling data. Nonlinearity in amplitude‐variation‐with‐offset becomes important in the presence of large relative changes in acoustic and elastic medium properties. A procedure for pre‐processing physical modelling reflection data is enacted on the reflection from a water‐plexiglas boundary. The resulting picked and processed amplitudes are compared with the exact solutions of the plane‐wave Zoeppritz equations, as well as approximations that are first, second, and third order in , , and . In the low angle range of 0°–20°, the third‐order plane‐wave approximation is sufficient to capture the nonlinearity of the amplitude‐variation‐with‐offset response of a liquid‐solid boundary with , , and ρ contrasts of 1485–2745 m/s, 0–1380 m/s, and 1.00–1.19 gm/cc respectively, to an accuracy value of roughly 1%. This is in contrast to the linear Aki–Richards approximation, which is in error by as much as 25% in the same angle range. Even‐order nonlinear corrective terms are observed to be primarily involved in correcting the angle dependence of , whereas the odd‐order nonlinear terms are involved in determining the absolute amplitude‐variation‐with‐offset magnitudes.     

Orthotropic anisotropy: on contributions of elasticity parameters to a difference in quasi-P-wave-squared velocities resulted from propagation in two orthogonal symmetry planes

We investigate the dependence of quasi P-wave phase velocity propagating in orthotropic media on particular elasticity parameters. Specifically, due to mathematical facilitation, we consider the squared-velocity difference, , resulted from propagation in two mutually perpendicular symmetry planes. In the context of the effective medium theory, may be viewed as a parameter evaluating the influence of cracks – embedded in the background medium – parallel to one or both aforementioned planes. Our investigation is both theoretical and numerical. Based on Christoffel's equations, we propose two accurate approximations of . Due to them, we interpret the aforementioned squared-velocity difference as being twice more dependent on , than on . To describe the magnitude of the dependence, we consider the proportions between the partial derivatives of . Further, it occurs that is influenced by the ratio of vertically propagating quasi P-wave to vertically propagating quasi S-wave. Anomalously high might be caused by the low P/S ratio, which in turn can be an indicator of the presence of gas in natural fractures or aligned porosity. Also, we carry out numerical sensitivity study, according to which is approximately twice more dependent on than on , twice more sensitive to than to , and equally dependent on as on . The dependence on and can be neglected, especially for small phase angles. We verify the approximations and perform the sensitivity study, using eight examples of the elasticity tensors.     

Some remarks on Q-compensated sparse deconvolution without knowing the quality factor Q

The subsurface media are not perfectly elastic, thus anelastic absorption, attenuation and dispersion (aka Q filtering) effects occur during wave propagation, diminishing seismic resolution. Compensating for anelastic effects is imperative for resolution enhancement. Q values are required for most of conventional Q-compensation methods, and the source wavelet is additionally required for some of them. Based on the previous work of non-stationary sparse reflectivity inversion, we evaluate a series of methods for Q-compensation with/without knowing Q and with/without knowing wavelet. We demonstrate that if Q-compensation takes the wavelet into account, it generates better results for the severely attenuated components, benefiting from the sparsity promotion. We then evaluate a two-phase Q-compensation method in the frequency domain to eliminate Q requirement. In phase 1, the observed seismogram is disintegrated into the least number of Q-filtered wavelets chosen from a dictionary by optimizing a basis pursuit denoising problem, where the dictionary is composed of the known wavelet with different propagation times, each filtered with a range of possible values. The elements of the dictionary are weighted by the infinity norm of the corresponding column and further preconditioned to provide wavelets of different values and different propagation times equal probability to entry into the solution space. In phase 2, we derive analytic solutions for estimates of reflectivity and Q and solve an over-determined equation to obtain the final reflectivity series and Q values, where both the amplitude and phase information are utilized to estimate the Q values. The evaluated inversion-based Q estimation method handles the wave-interference effects better than conventional spectral-ratio-based methods. For Q-compensation, we investigate why sparsity promoting does matter. Numerical and field data experiments indicate the feasibility of the evaluated method of Q-compensation without knowing Q but with wavelet given.     

A fast algorithm for 3D azimuthally anisotropic velocity scan

The conventional velocity scan can be computationally expensive for large‐scale seismic data sets, particularly when the presence of anisotropy requires multiparameter scanning. We introduce a fast algorithm for 3D azimuthally anisotropic velocity scan by generalizing the previously proposed 2D butterfly algorithm for hyperbolic Radon transforms. To compute semblance in a two‐parameter residual moveout domain, the numerical complexity of our algorithm is roughly as opposed to of the straightforward velocity scan, with N being the representative of the number of points in a particular dimension of either data space or parameter space. Synthetic and field data examples demonstrate the superior efficiency of the proposed algorithm.     

Time‐lapse pre‐stack seismic data registration and inversion for CO2 sequestration study at Cranfield

Pre‐stack seismic data are indicative of subsurface elastic properties within the amplitude versus offset characteristic and can be used to detect elastic rock property changes caused by injection. We perform time‐lapse pre‐stack 3‐D seismic data analysis for monitoring sequestration at Cranfield. The time‐lapse amplitude differences of Cranfield datasets are found entangled with time‐shifts. To disentangle these two characters, we apply a local‐correlation‐based warping method to register the time‐lapse pre‐stack datasets, which can effectively separate the time‐shift from the time‐lapse seismic amplitude difference without changing the original amplitudes. We demonstrate the effectiveness of our registration method by evaluating the inverted elastic properties. These inverted time‐lapse elastic properties can be reliably used for monitoring plumes.     

Non‐convex compressed sensing with frequency mask for seismic data reconstruction and denoising

Compressed Sensing has recently proved itself as a successful tool to help address the challenges of acquisition and processing seismic data sets. Compressed sensing shows that the information contained in sparse signals can be recovered accurately from a small number of linear measurements using a sparsity‐promoting regularization. This paper investigates two aspects of compressed sensing in seismic exploration: (i) using a general non‐convex regularizer instead of the conventional one‐norm minimization for sparsity promotion and (ii) using a frequency mask to additionally subsample the acquired traces in the frequency‐space () domain. The proposed non‐convex regularizer has better sparse recovery performance compared with one‐norm minimization and the additional frequency mask allows us to incorporate a priori information about the events contained in the wavefields into the reconstruction. For example, (i) seismic data are band‐limited; therefore one can use only a partial set of frequency coefficients in the range of reflections band, where the signal‐to‐noise ratio is high and spatial aliasing is low, to reconstruct the original wavefield, and (ii) low‐frequency characteristics of the coherent ground rolls allow direct elimination of them during reconstruction by disregarding the corresponding frequency coefficients (usually bellow 10 Hz) via a frequency mask. The results of this paper show that some challenges of reconstruction and denoising in seismic exploration can be addressed under a unified formulation. It is illustrated numerically that the compressed sensing performance for seismic data interpolation is improved significantly when an additional coherent subsampling is performed in the domain compared with the domain case. Numerical experiments from both simulated and real field data are included to illustrate the effectiveness of the presented method.     

The rotated Cartesian coordinate method to remove the axial singularity of cylindrical coordinates in finite‐difference schemes for elastic and viscoelastic waves

When modelling the propagation of 3D non‐axisymmetric elastic and viscoelastic waves in cylindrical coordinates using the finite‐difference time‐domain method, a mathematical singularity occurs due to the presence of terms in the elastic and viscoelastic wave equations. For many years, this issue has been impeding the accurate numerical solution near the axis. In this work, we propose a simple but effective method for the treatment of this numerical singularity problem. By rotating the Cartesian coordinate system around the z‐axis in cylindrical coordinates, the numerical singularity problems in both 2D and 3D cylindrical coordinates can be removed. This algorithm has three advantages over the conventional treatment techniques: (i) the excitation source can be directly loaded at , (ii) the central difference scheme with second‐order accuracy is maintained, and (iii) the stability condition at the axis is consistent with the finite‐difference time‐domain in Cartesian coordinates. This method is verified by several 3D numerical examples. Results show that the rotating the Cartesian coordinate method is accurate and stable at the singularity axis. The improved finite‐difference time‐domain algorithm is also applied to sonic logging simulations in non‐axisymmetric formations and sources.     

Robust f‐x projection filtering for simultaneous random and erratic seismic noise attenuation

Linear prediction filters are an effective tool for reducing random noise from seismic records. Unfortunately, the ability of prediction filters to enhance seismic records deteriorates when the data are contaminated by erratic noise. Erratic noise in this article designates non‐Gaussian noise that consists of large isolated events with known or unknown distribution. We propose a robust f‐x projection filtering scheme for simultaneous erratic noise and Gaussian random noise attenuation. Instead of adopting the ?₂‐norm, as commonly used in the conventional design of f‐x filters, we utilize the hybrid ‐norm to penalize the energy of the additive noise. The estimation of the prediction error filter and the additive noise sequence are performed in an alternating fashion. First, the additive noise sequence is fixed, and the prediction error filter is estimated via the least‐squares solution of a system of linear equations. Then, the prediction error filter is fixed, and the additive noise sequence is estimated through a cost function containing a hybrid ‐norm that prevents erratic noise to influence the final solution. In other words, we proposed and designed a robust M‐estimate of a special autoregressive moving‐average model in the f‐x domain. Synthetic and field data examples are used to evaluate the performance of the proposed algorithm.     

Zur Bestimmung des gesamten anorganischen Kohlenstoffes in natürlichen Gewässern durch Titration mit Salzsäure

The knowledge of total inorganic carbon concentration (c) is important for characterizing natural waters. It is usually measured by the titration alkalinity (“m-value”) and pH which depend on temperature and ionic strength. This paper demonstrates that Ca (and Mg) can influence the calculation of from titration alkalinity, too. Errors result from neglecting this influence. In such cases the share of ion pair CaCO amounts to more than 50 % of . General relationships among the influencing factors are given by tables calculated with the help of a BASIC computer programme for calculation from titration alkalinity, pH and Ca concentration.     

Gewinnung biokinetischer Daten von Altwässern mit dem Respirationsfermentor. Teil 1: Methodische Grundlagen und die Bestimmung der spezifischen Wachstumsund Sauerstoffaufnahmerate

Es wird eine Fermentationsvorrichtung für die Messung des Sauerstoffverbrauchs bei wählbaren Sauerstoffpartialdrucken vorgestellt, die in einem geschlossenen Gaskreislauf einen Laboratoriumsfermentor mit Meß- und Regeleinrichtungen sowie Möglichkeiten der Probenahme aufweist. Die Druckregression durch die Kohlendioxidadsorption steuert die synchron registrierte Nachlieferung von Sauerstoffgas. Durch diskontinuierliche Kultivation eines adaptierten Belebtschlammes wurden von einem Braunkohlenpyrolyseabwasser die Sauerstoffverbrauchskurven aufgenommen. Es wurde ein polyauxischer Kurvenverlauf erhalten, wobei auf den Abbau der niederen Fettsauren, im wesentlichen Azetat, 70 … 80% des Gesamtsauerstoffverbrauchs entfielen. Hierfür werden biokinetische Daten bestimmt; ; ; μ_max = 0,56 h^–1 (ausgewählte Werte). Die bei der Erfassung und Errechnung der Daten auftretenden Probleme wurden diskutiert.     

Identifying gravity anomalies caused by granitic intrusions in Nanling mineral district,China: a multifractal perspective

Several power‐law relationships of geophysical potential fields have been discussed recently with renewed interests, including field value–distance () and power spectrum–wavenumber () models. The singularity mapping technique based on the density/concentration–area (C–A) power‐law model is applied to act as a high‐pass filter for extracting gravity and magnetic anomalies regardless of the background value and to detect the edges of gravity or magnetic sources with the advantage of scale invariance. This is demonstrated on a synthetic example and a case study from the Nanling mineral district, Southern China. Compared with the analytic signal amplitude and total horizontal gradient methods, the singularity mapping technique provides more distinct and less noisy boundaries of granites than traditional methods. Additionally, it is efficient for enhancing and outlining weak anomalies caused by concealed granitic intrusions, indicating that the singularity method based on multifractal analysis is a potential tool to process gravity and magnetic data.     

Transmission + reflection anisotropic wave-equation traveltime and waveform inversion

A transmission + reflection wave-equation traveltime and waveform inversion method is presented that inverts the seismic data for the anisotropic parameters in a vertical transverse isotropic medium. The simultaneous inversion of anisotropic parameters and ε is initially performed using transmission wave-equation traveltime inversion method. Transmission wave-equation traveltime only provides the low-intermediate wavenumbers for the shallow part of the anisotropic model; in contrast, reflection wave-equation traveltime estimates the anisotropic parameters in the deeper section of the model. By incorporating a layer-stripping method with reflection wave-equation traveltime, the ambiguity between the background-velocity model and the depths of reflectors can be greatly mitigated. In the final step, multi-scale full-waveform inversion is performed to recover the high-wavenumber component of the model.  We use a synthetic model to illustrate the local minima problem of full-waveform inversion and how transmission and reflection wave-equation traveltime can mitigate this problem. We demonstrate the efficacy of our new method using field data from the Gulf of Mexico.     

Separation of multi-mode surface waves by supervised machine learning methods

Logistic regression, neural networks and support vector machines are tested for their effectiveness in isolating surface waves in seismic shot records. To distinguish surface waves from other arrivals, we train the algorithms on three distinguishing features of surface-wave dispersion curves in the domain: spectrum coherency of the trace's magnitude spectrum, local dip and the frequency range for a fixed wavenumber k in the spectrum. Numerical tests on synthetic data show that the kernel-based support vector machines algorithm gives the highest accuracy in predicting the surface-wave window in the domain compared to neural networks and logistic regression. This window is also used to automatically pick the fundamental dispersion curve. The other two methods correctly pick the low-frequency part of the dispersion curve but fail at higher frequencies where there is interference with higher-order modes.     

Multi-parameter reflection waveform inversion for acoustic transversely isotropic media with a vertical symmetry axis

Full waveform inversion in transversely isotropic media with a vertical symmetry axis provides an opportunity to better match the data at the near and far offsets. However, multi-parameter full waveform inversion, in general, suffers from serious cycle-skipping and trade-off problems. Reflection waveform inversion can help us recover a background model by projecting the residuals of the reflected wavefield along the reflection wavepath. Thus, we extend reflection waveform inversion to acoustic transversely isotropic media with a vertical symmetry axis utilizing the proper parameterization for reduced parameter trade-off. From a radiation patterns analysis, an acoustic transversely isotropic media with a vertical symmetry axis is better described by a combination of the normal-moveout velocity and the anisotropic parameters η and δ for reflection waveform inversion applications. We design a three-stage inversion strategy to construct the optimal resulting model. In the first stage, we only invert for the background by matching the simulated reflected wavefield from the perturbations of and δ with the observed reflected wavefield. In the second stage, the background and η are optimized simultaneously and the far-offset reflected wavefield mainly contribute to their updates. We perform Born modelling to compute the reflected wavefield for the two stages of reflection waveform inversion. In the third stage, we perform full waveform inversion for the acoustic transversely isotropic media with a vertical symmetry axis to delineate the high-wavenumber structures. For this stage, the medium is described by a combination of the horizontal velocity , η and ε instead of , η and δ. The acoustic multi-parameter full waveform inversion utilizes the diving waves to improve the background as well as utilizes reflection for high-resolution information. Finally, we test our inversion algorithm on the modified Sigsbee 2A model (a salt free part) and a two-dimensional line from a three-dimensional ocean bottom cable dataset. The results demonstrate that the proposed reflection waveform inversion approach can recover the background model for acoustic transversely isotropic media with a vertical symmetry axis starting from an isotropic model. This recovered background model can mitigate the cycle skipping of full waveform inversion and help the inversion recover higher resolution structures.     

3D angle-domain common-image gathers for migration velocity analysis

Angle‐domain common‐image gathers (ADCIGs) are an essential tool for migration velocity analysis (MVA). We present a method for computing ADCIGs in 3D from the results of wavefield‐continuation migration. The proposed methodology can be applied before or after the imaging step in a migration procedure. When computed before imaging, 3D ADCIGs are functions of the offset ray parameters (p, p) ; we derive the geometric relationship that links the offset ray parameters to the aperture angle γ and the reflection azimuth φ. When computed after imaging, 3D ADCIGs are directly produced as functions of γ and φ. The mapping of the offset ray parameters (p, p) into the angles (γ, φ) depends on both the local dips and the local interval velocity; therefore, the transformation of ADCIGs computed before imaging into ADCIGs that are functions of the actual angles is difficult in complex structure. By contrast, the computation of ADCIGs after imaging is efficient and accurate even in the presence of complex structure and a heterogeneous velocity function. On the other hand, the estimation of the offset ray parameters (p, p) is less sensitive to velocity errors than the estimation of the angles (γ, φ). When ADCIGs that are functions of the offset ray parameters (p, p) are adequate for the application of interest (e.g. ray‐based tomography), the computation of ADCIGs before imaging might be preferable. Errors in the migration velocity cause the image point in the angle domain to shift along the normal to the apparent geological dip. By assuming stationary rays (i.e. small velocity errors), we derive a quantitative relationship between this normal shift and the traveltime perturbation caused by velocity errors. This relationship can be directly used in an MVA procedure to invert depth errors measured from ADCIGs into migration velocity updates. In this paper, we use it to derive an approximate 3D residual moveout (RMO) function for measuring inconsistencies between the migrated images at different γ and φ. We tested the accuracy of our kinematic analysis on a 3D synthetic data set with steeply dipping reflectors and a vertically varying propagation velocity. The tests confirm the accuracy of our analysis and illustrate the limitations of the straight‐ray approximation underlying our derivation of the 3D RMO function.     

Effect of micro‐inhomogeneity on the effective stress coefficients and undrained bulk modulus of a poroelastic medium: a double spherical shell model

Although most rocks are complex multi‐mineralic aggregates, quantitative interpretation workflows usually ignore this complexity and employ Gassmann equation and effective stress laws that assume a micro‐homogeneous (mono‐mineralic) rock. Even though the Gassmann theory and effective stress concepts have been generalized to micro‐inhomogeneous rocks, they are seldom if at all used in practice because they require a greater number of parameters, which are difficult to measure or infer from data. Furthermore, the magnitude of the effect of micro‐heterogeneity on fluid substitution and on effective stress coefficients is poorly understood. In particular, it is an open question whether deviations of the experimentally measurements of the effective stress coefficients for drained and undrained elastic moduli from theoretical predictions can be explained by the effect of micro‐heterogeneity. In an attempt to bridge this gap, we consider an idealized model of a micro‐inhomogeneous medium: a Hashin assemblage of double spherical shells. Each shell consists of a spherical pore surrounded by two concentric spherical layers of two different isotropic minerals. By analyzing the exact solution of this problem, we show that the results are exactly consistent with the equations of Brown and Korringa (which represent an extension of Gassmann's equation to micro‐inhomogeneous media). We also show that the effective stress coefficients for bulk volume α, for porosity n_? and for drained and undrained moduli are quite sensitive to the degree of heterogeneity (contrast between the moduli of the two mineral components). For instance, while for micro‐homogeneous rocks the theory gives n_? = 1, for strongly micro‐inhomogenous rocks, n_? may span a range of values from –∞ to ∞ (depending on the contrast between moduli of inner and outer shells). Furthermore, the effective stress coefficient for pore volume (Biot–Willis coefficient) α can be smaller than the porosity ?. Further studies are required to understand the applicability of the results to realistic rock geometries.     

Research Note: Experimental measurements of Q‐contrast reflections

While seismic reflection amplitudes are generally determined by real acoustical impedance contrasts, there has been recent interest in reflections due to contrasts in seismic‐Q. Herein we compare theoretical and modelled seismic reflection amplitudes for two different cases of material contrasts. In case A, we examine reflections from material interfaces that have a large contrast in real‐valued impedance () with virtually no contrast in seismic‐Q. In case B, we examine reflections from material interfaces that have virtually no contrast in but that have very large seismic‐Q contrasts. The complex‐valued reflection coefficient formula predicts non‐zero seismic reflection amplitudes for both cases. We choose physical materials that typify the physics of both case A and case B. Physical modelling experiments show significantly large reflections for both cases – with the reflections in the two cases being phase shifted with respect to each other, as predicted theoretically. While these modelling experiments show the existence of reflections that are predicted by theory, there are still intriguing questions regarding the size of the Q‐contrast reflections, the existence of large Q‐contrast reflections in reservoir rocks and the possible application of Q‐reflection analysis to viscosity estimation in heavy oilfields.