Well tie calculus

Well tie calculus

Synthetic seismograms can be created by doing basic calculus on travel-time functions. Integrating slowness (the reciprocal of velocity) yields a time–depth relationship. Differentiating acoustic impedance (velocity times density) yields a reflectivity function along the borehole. In effect, the integral tells us where a rock interface is positioned in the time domain, while the derivative tells us how the seismic wavelet will be scaled.

The code and data required to reproduce the results and figures in this tutorial are in an IPython Notebook at [1]. For this tutorial, we will make use of the open-source Python libraries NumPy (for computation niceties), and Matplotlib (for graphics).

As an example, I will use the sonic and density log from the Penobscot L-30 well offshore Nova Scotia. It has a standard suite of petrophysical logs, but we’ll just be using P-wave sonic and bulk density. This vertical well is in the same volume of 3D seismic data that Matt Hall used in his horizon smoothing tutorial.

Dealing with incomplete, noisy logs

To deal with the shallow section, we first need to adjust the depths relative to sea level by subtracting the KB elevation (30.2 m) from the measured depths.

If we now integrate the sonic log, we see that time = 0 s corresponds to a depth of 347 m TVDSS. To position the top of the log at the correct travel time on the seismic section, we need to use the thickness and replacement velocities for both the water column (137.5 m at 1480 m/s), and the section above the log (179.8 m at 1600 m/s) respectively.

Both the sonic and density curves have some spikes and null values that need to be dealt with. To cope with these, I applied a rolling median filter to identify spikes, and then clipped them.

Computing the time to depth relationship

The time–depth relationship is obtained by scaling the sonic log by the sample interval (6 inches or 0.1524 m), integrating it with np.cumsum(), and shifting it to the start time of the sonic log.

scaled_sonic = 0.1524 * DT / 1e6
t_int = 2 * np.cumsum(scaled_sonic)
td = t_int + log_start_time

As you can see in Figure 1, the time–depth relationship is good down to the Logan Canyon Formation. Beneath that, however, the Lower Missisauga, Abenaki, and Mid Baccaro Formations are mis-positioned by a half-cycle or so in the fast direction. We could further adjust the time–depth relationship by stretching and squeezing between anchor points and so on, but this would break the direct equivalency to the measured slowness log. As Rachel Newrick (2012) wrote, “with each modification to the synthetic we should think about why we are applying a certain process, what the effect is, and whether it makes sense.” For now, we will live with the mistie.

Computing the reflection coefficients in time

The P-wave normal incidence reflection coefficient at a sample i is given by

$R_{i}={\frac {Z_{i+1}-Z_{i}}{Z_{i+1}+Z_{i}}}$

This doesn’t have to be computed over and over again, iterating over the impedance array in a for loop, say. Since we have the acoustic impedance function stored as a NumPy array, we can slice the array with a one-sample shift to carry out the calculation in a single operation:

RC = (Z[1:] - Z[:-1]) / (Z[1:] + Z[:-1])

To make a synthetic seismogram we’d like to convolve this reflection coefficient series with a seismic wavelet. But seismic is recorded in the time domain, and the convolution operator requires inputs to have the same sample rate, so we must resample the reflection coefficients using our time–depth relationship. First, we define a regularly sampled time vector t that goes from 0 to 3 seconds with a sample rate of 4 ms:

t = np.arange(0, 3, 0.004)

Next, to obtain the reflection coefficients as a function of time RC_t, we pass three inputs to NumPy’s interpolation function: the regularly sampled time axis t, the time to depth relationship td, and the depth-domain reflectivity series RC:

RC_t = np.interp(t, td[:-1], RC)

The synthetic seismic experiment

We need a wavelet. In the IPython Notebook, you can see how to create a Ricker wavelet w with a dominant frequency of 25 Hz and convolve that with the reflection coefficient series:

synthetic = np.convolve(w, RC_t, mode='same')

We chose this Ricker wavelet because 25 Hz sits more or less in the middle of the seismic band (not shown) within the interval encountered by the well log. Furthermore, the water bottom reflection is a peak (black) with symmetric side lobes, which indicates that the data is normal polarity and close to zero phase.

To view the tie we plot the synthetic as a wiggle trace upon the real seismic amplitudes.

When P-wave travel-time from a sonic log is appropriately scaled, the cumulative sum (integral) of slowness gives the two-way-time to depth depth relationship. The stepwise differencing (derivative) of impedance, per the reflectivity equation, gives the reflectivity function used to create the synthetic seismogram for the well-tie. Positive amplitudes (peak) are shown in black on the seismic filled on the synthetic, negative amplitudes (trough) are white.

Final thoughts

If you find yourself stretching or squeezing a time–depth relationship in order to make synthetic events better align with seismic events, take the time to compute the implied corrections to the well logs. Differentiate the new time–depth curve. How much have the interval velocities changed? Are the rock properties still reasonable? Synthetic seismograms should adhere to the simple laws of calculus — and not imply unphysical versions of the earth.

External links

• IPython Notebook — includes all the code
• Madagascar workflow - reproducible using Madagascar open-source software
• Corresponding author: Evan Bianco, evanagilegeoscience.com

find literature about Well tie calculus