Michael Jervis

Determining bed boundaries from inversion of EM logging data using general measures of model structure and data misfit

In this paper we propose a new method to locate bed boundaries by carrying out a 1-D nonlinear inversion of electromagnetic (EM) logging data. We first solve for 1-D resistivity structure in which the earth is modeled using layers of constant thickness. This thickness is determined based upon the tool resolution and the desired resolution from the user. We use general measures of data misfit and model structure in the inversion to construct piecewise‐constant models through the iteratively reweighted least‐squares (IRLS) procedure, and we minimize a generic global model objective function subject to data constraints. The general measure includes the traditional lp norm as a special case. The recovered piecewise‐constant models simulate traditional models comprised of a few uniform layers and hence permit easy determination of bed boundaries. The generic model objective function allows the incorporation of prior geological information and provides measures of the closeness to the reference model and the am...

Reconstruction of resistivity structure from 2-D inversions of resistivity logging data using equality and inequality constraints

The inversion technique is a powerful tool in estimating residual and movable hydrocarbon and has been widely used in the oil industry with great su ccess. However, geophysical inverse problems are ill-posed b ecause of incomplete measurements in the survey, noise in the data, various approximations in the physical models, and numerical errors in the computations. The solutions from these ill-posed problems are nonunique and unstable. If the problem is nonlinear, then the inversion is further complicated by local minima. To overcome these difficulties, it is necessary to regularize the inverse problems and apply appropriate constraints on the solutions. In this paper, we develop a model-norm-based inverse algorithm to estimate 2-D resistivity distribution around the borehole using resistivity logging data. This algorithm features a model objective function that contains a priori information about the geological structure and an outside penalty function that augments the unconstrained inversion with equality and inequality constraints on the solutions. This penalty function is introduced through a weighting parameter that is increased from a starting value to a theoretically infinite value. The linear system of equations to be solved at each iteration in the inversion is preconditioned to avoid unstable and slow convergence caused by the large weighting on the penalty term. In this paper, we show the usefulness of two particular types of inequality constraints, the range and the relative constraints. The range constraint defines the feasible region for the model parameters, and the relative constraint dictates the relative strength among the model parameters. We apply our algorithm to both synthetic and field data.

Analysis of microseismic events during a multistage hydraulic stimulation experiment at a shale gas reservoir

Despite the current easing in demand for increased oil production linked to the global downturn in crude prices, energy demand continuously increases and the long-term demand will require maximizing the productivity of reservoirs and a search into the exploitation of new resources in increasingly challenging environments. In this study, we present the results from the monitoring of the very first multistage stimulation experiment at a shale gas reservoir in Saudi Arabia, presenting an analysis of the microseismicity induced during the treatment. Our aim was to analyse microseismic events to better understand fracture growth and the role of pre-existing fractures in these reservoirs. Microseismic (MS) event monitoring is used to track the creation of fractures during and after the stimulation, and therefore to evaluate the effect of the reservoir stimulation. The monitoring includes a downhole array of 12 3C-sensors that were deployed in a vertical well with a 30.5 m level spacing. A total of 415 MS events were located and analysed, with the results outlining induced fractures extending consistently with an average azimuth of N335° E, normal to the horizontal section of the treatment well. This implies that there are no changes in the local stress direction along the treatment well either in situ or induced along the treatment. There are significant changes in total length and aspect ratio (length/width) of the fractures induced in the different stages. These variations could be attributed to in situ fracturing, local rock heterogeneity or the influence of the treatment parameters. In general, early and late stages of stimulation show the longest fracture networks, with events induced further away from the initiation point. We found no immediate relationship between treatment parameters (peak pressure and pumping rates) and fracture extension. Sensitivity analysis using Monte Carlo simulation methods shows a higher location uncertainty for events located at the early stages, thus limiting the interpretation from monitored seismicity in the early stages. An analysis of magnitude distribution with distance shows a decrease in sensitivity of one degree of magnitude for every 375 m, and a maximum viewing distance of approximately 700 m for the current set-up. The low number of located events does not provide a complete enough dataset for a robust analysis of changes in b-value (slope in linear part of magnitude distribution) during the treatment: however, magnitude distributions, corrected for array sensitivity, provide a useful variable for the validation of geomechanical models currently being developed for the reservoir.

Robust and efficient 3-D traveltime calculations in arbitrary media

Traveltimes are calculated for 3D acoustic media with arbitrary velocity variations. The scheme used here is developed from the 2D dynamic programming approach of Schneider et al. (1992) as modified by Faria and Stoffa (1994). This approach, based on Fermat`s principle, uses simple calculus techniques and ordering of the calculations to determine first arrival times on a rectangular grid. The initial condition requires that traveltimes are computed to the nearest neighbor points surrounding the source location. The computational wavefronts then proceed outward from the source in the 4 principal horizontal directions before the computations are repeated in the vertical directions and again in the horizontal directions to account for turning rays. Though the amount of computation is large for 3D media, this algorithm is easily vectorized. In addition, the model does not require any smoothing, which is an advantage over upwind finite-difference schemes, when the velocity model contains salt or other high velocity media.

On the use of microtremors for hydrocarbon detection

The link of spectral anomalies of microtremors to underlying hydrocarbon reservoirs is very controversial, as field experiments support both positive and negative opinions, and there is not a solid theory supporting this work hypothesis. We conducted field tests at different sites, with and without oil and gas presence, to add new experimental data to the ongoing studies. Microtremor information may become repeatable (and so physically meaningful) only when the observation duration exceeds a few days, but even in this case, factors such as topography and active faults may severely bias the signal. Ocean waves impinging the coasts provide natural background noise, which stands out clearly when the observation time exceeds a dozen days or so, in such a way that human noise is stacked out statistically over time. Microtremors recorded in (relatively) deep wells may provide useful information about ongoing production in a reservoir, and may link well data and seismic surveys, as their interferometric analysis can provide information comparable to Vertical Seismic Profiles.

High resolution crosswell reflection imaging in the presence of anisotropy: The Santa Rosa gas field, Eastern Venezuela

Seismic anisotropy in sand/shale sequences causes structural imaging problems when an isotropic earth is assumed. Three high resolution crosswell seismic profiles were acquired in the Santa Rosa gas field in eastern Venezuela. The main objective is to gain detailed structural images in a seismically difficult data zone, where the deep target horizons are obscured by a near-surface gas column. Isotropic imaging techniques were applied to the crosswell data with relatively poor results. The imaging target comprises sands from 5 to 50 ft thick in a 1800 ft thick deltaic, mainly shaley sequence unconformably overlying a predominantly fluvial Oligocene sand-rich sequence. The abundance of shales and fine layering has resulted in significant anisotropy. We have performed perhaps the first anisotropic reflection imaging of crosswell seismic data. Due to the structure and well deviations we use a fully 3-D modeling and imaging framework to produce both velocity and crosswell reflection images. Initial tomographic inversion revealed strong evidence of anisotropy, which correlated well with sonic log information. Differences between vertical and horizontal velocities range from 1623% in the upper shale-rich sequence, dropping to less than 5% in the lower sands.

3-D Traveltime Calculations

This equation is the basis for both the 2-D and 3-D traveltime code. Figure 1 shows that the update time t is the time of the minimum raypath, t0 for example, added to the time from z0 to the update point at (x,z2). In other words t = t0 + s[(z2 – z0) + Δx2]1/2 To find the minimum raypath location z0, this equation is differentiated w.r.t. z0 giving dt/dz0 = z0w/t0 – s(z2 – z0)[(z2 – z0) + Δx2]-1/2

Simultaneous recovery of the source wavelet and seismic impedances by nonlinear elastic inversion

Previous work in nonlinear seismic waveform inversion has assumed a known source (Crase, 1987). I am including the source time function as a parameter in a preconditioned gradient algorithm along with the p-wave and s-wave impedances. Since local methods of optimization are used, the reference or background model of seismic velocities is assumed known. Tarantola (1988) shows that the gradients of the misfit function with respect to the elastic parameters, density and the source time function can be expressed in terms of a backpropagating residual field. For either a moment density (marine data) or a force density source (land data), the source gradient is simply the field recorded at the source location during backpropagation of the data residuals i.e propagation of the differences between the observed data and the data calculated using the current model of the subsurface. Either pressure (marine case) or displacement (land case) is recorded in the appropriate time window as the residuals are input into the receiver locations as secondary sources. By perturbing each parameter subspace, or parameter type, with it’s respective gradient and solving the forward problem again for each type of parameter type, a new inverse problem is formulated for solution of the step lengths. The optimum step lengths can be determined using the subspace methods of Kenneth and Williamson (1987). We have a subspace defined by the individual descent vectors (gradients), for the seismic impedances and the source, and an optimal combination of the these vectors can be found which gives the greatest reduction in misfit (Sambridge et al. 1990).