Search

One of the main research questions regarding archaeological molluscan purple pigments and dyes is whether it is possible to determine which malacological species produced these colorants. For this determination of the zoological provenance of the pigment, a multicomponent analysis must be performed, which can only be obtained from the HPLC technique—the optimal method for identifying all the detectable colorants in a sample. In order to find any trends in the compositions of the dye components from various species of purple-producing sea snails, a statistical formulation is needed. Though principal component analysis (PCA) is a powerful statistical tool that has been used in the analysis of these components, it is based on an algorithm that combines all the componential values and produces new two-dimensional parameters whereby the individualities of the original dye component values are lost. To maintain the integrity of the dye compositions in the purple pigments, a very simple formulation was first published in 2008 and applied to a limited number of samples. This property is known as DMI (short for Di-Mono Index), and for each sample, it is simply the ratio of the peak area of DBI relative to that of MBI, evaluated at the standard wavelength of 288 nm, which has been used for such peak calculations. Currently, considerably more modern and archaeological pigments have been analyzed via HPLC; thus, in the current study, the DMI has been expanded to characterize these purple pigments. Furthermore, a ternary diagram comprising the blue, violet, and red components that can be found in purple colorants is presented for both modern and archaeological purple pigments from the three Muricidae species known in antiquity to produce purple pigments. This triangular diagram is intuitive, retains the integrity of the original dyes, and is presented here for the first time. Both the DMI and the ternary diagram can discern whether a Hexaplex trunculus species or perhaps the Bolinus brandaris or Stramonita haemastoma species were used to produce the pigment. Further, these two representations can also determine whether the IND-rich or the DBI-rich varieties, or both, of H. trunculus were used to produce the purple pigment, either as a paint pigment or as a textile dye.

International audience; Research into the sustainability of natural, potentially renewable, resources is one of the major issues of our time. It naturally includes the quest for sustainable sources of colorants for textiles, cosmetics, and food. In industrialized countries, natural dyeing with plants and a few species of coccid insects was practiced on a large scale for centuries before synthetic colorants were developed. Therefore, historical documents on the growing of dye plants and dyeing processes offer a relevant basis from which to start reconsidering the potential of natural colorants in our time. However, written sources need to be completed by experimental archaeologists to allow a scientific understanding of the biochemical reactions at work in the historical processes described. The results of such interdisciplinary research can then inspire contemporary programs to revive the production of natural dyes. The long history of dyeing blue with woad, Isatis tinctoria L., is revisited here as an illustration of the fruitful complementarity of sources and approaches. This article presents a step-by-step re-assessment of the production chain of woad as described in historical texts, from the growing of the plant to its use as a source of indigo in the woad and indigo vats. The experimental reconstitution of the processing of woad leaves into couched woad allowed us to follow the evolution of the composition and proportions of indigoid colorants in the leaves by HPLC analyses. Additionally, HPLC analyses allowed a comparison of the respective indigoid contents of couched woad and sukumo, the form of indigo dye resulting from another couching process, traditionally used in Japan for dyers’ knotweed, Persicaria tinctoria (Ait.) H. Gross. The reconstitution of the 18th century woad and indigo vat process allowed investigations into the bacterial flora associated with the use of couched woad in vat liquors, which were found to contain different indigo-reducing bacteria, including two distinct strains of a new indigo-reducing species.

SUMMARYReflector-normal angles and reflector-curvature parameters are the principal geometric attributes used in seismic interpretation for characterizing the orientations and shapes, respectively, of geological reflecting surfaces. Commonly, the input data set for their computation consists of fine 3-D grids of scalar fields representing either the seismic-driven reflectivities (e.g. amplitudes of 3-D seismic migrated volumes) or model-driven reflectivities, computed, for example, from the derived elastic impedance parameters. Conventionally, the computation of curvature parameters at each gridpoint is based on analysing the local change in the inline/crossline dips, considering the potential existence of a local quadratic reflecting surface in the vicinity of that point. This assumption breaks down for subsurface points in the vicinity of either complex reflecting surfaces (e.g. brittle/rough/tilted synclines/anticlines, ridges/troughs and saddles) and/or sharp, discontinuous geological features (e.g. fault edges/tips, pinch-outs, fracture systems, channels and small geobodies), where the values of the computed curvature become extremely high. However, while these high values can indicate the existence of non-reflecting objects, they do not deliver their specific geometric characteristics. In this study, we present a novel method that better characterizes the shapes of these complex geological features by extending the assumption of local surfaces (2-D surfaces in 3-D space) into local hypersurfaces (3-D hypersurfaces in 4-D space), with their corresponding (three rather than two) principal (and effective) curvature parameters. We demonstrate the advantages of our method by comparing the conventional dip-based surface curvature parameters with the hypersurface curvature parameters, using a synthetic model/image with different types and shapes of geological features and a seismic image of real data containing a complex fault and hidden buried channels.

SUMMARY We present an original, generic and efficient approach for computing the first and second partial derivatives of ray (group) velocities along seismic ray paths in general anisotropic (triclinic) elastic media. As the ray velocities deliver the ray element traveltimes, this set of partial derivatives constructs the so-called kinematic and dynamic sensitivity kernels which are used in different key seismic modelling and inversion methods, such as two-point ray bending methods and seismic tomography. The second derivatives are useful in the solution of the above-mentioned kinematic problems, and they are essential for evaluating the dynamic properties along the rays (amplitudes and phases). The traveltime is delivered through an integral over a given Lagrangian defined at each point along the ray. In our approach, we use an arclength-related Lagrangian representing a reciprocal of the ray velocity magnitude. Although this magnitude cannot be explicitly expressed in terms of the medium properties and the ray direction components, its derivatives can still be formulated analytically using the corresponding arclength-related Hamiltonian that can be explicitly expressed in terms of the medium properties and the slowness vector components; this requires first to obtain (invert for) the slowness vector components, given the ray direction components. Computation of the slowness vector and the ray velocity derivatives is considerably simplified by using an auxiliary scaled-time-related Hamiltonian obtained directly from the Christoffel equation and connected to the arclength-related Hamiltonian by a simple scale factor. This study consists of two parts. In Part I, we consider general anisotropic (triclinic) media, and provide the derivatives (gradients and Hessians) of the ray velocity, with respect to (1) the spatial location and direction vectors and (2) the elastic model parameters. The derivatives are obtained for both quasi-compressional and quasi-shear waves, where other types of media, characterized with higher symmetries, can be considered particular cases. In Part II, we apply the theory of Part I explicitly to polar anisotropic media (transverse isotropy with tilted axis of symmetry, TTI), and obtain the explicit ray velocity derivatives for the coupled qP and qSV waves and for SH waves. The derivatives for polar anisotropy are simplified (as compared to general anisotropy), obviously yielding more effective computations. The ray velocity derivatives are tested by checking consistency between the proposed analytical formulae and the corresponding numerical ones.

Highlights • Multilevel Composition is an innovative method involving color composition and co-rendering of multilevel attribute maps. • It is useful for characterizing multi-depth geological features based on their spatiotemporal distribution within three-dimensional seismic data. • The technique produces a single image map, in which inter-window/layer depth information is coded in colors for reliable representation of the actual geology. • In the eastern Nile fan, it was applied to visualize and resolve the complexities of buried clastic deep-water depositional elements. • On the Omakere Ridge, it successfully illuminated seafloor seeps and reveals their link to deeper fluid-bearing intervals. Abstract Advanced seismic data and multi-attribute visualization techniques, such as color blending of attributes, have considerably enhanced the capability of interpreters to characterize geological features in three-dimensional (3D) seismic reflection datasets. However, high resolution investigation of complex, vertically linked geological features such as channel systems and fluid conduits, remains challenging. These features may appear in the dataset as pronounced attribute anomalies, such as high-amplitude or spectrally or structurally enhanced seismic reflectivity bands, at several depth levels. Vertical linkages between these features, however, may not be readily established. We have developed an innovative method, Multilevel Composition, for an intuitive display of vertically connected features. Our method involves the composition of attribute maps from three different depth/time windows or slices onto a single map, in which inter-window/layer depth information is coded in colors. Multilevel Composition starts with the identification of suitable seismic attributes, such as high amplitudes in the examples displayed here, to map features of geological interest. At least one reference horizon is then identified and mapped in the vicinity of the target window of interest. Three sub-windows are then defined with respect to the reference horizon(s) based on the vertical and spatial distribution of the geological features. Relevant seismic attributes are computed for each of the sub-windows, and the resulting maps, one from each sub-window, are assigned basic color channels and are co-rendered to reveal multilevel linkages between these features. We demonstrate the efficacy of this method by applying it to two 3D seismic datasets, one illuminating deep-water depositional elements in the eastern Nile fan, eastern Mediterranean and the other targeting seafloor seeps and underlying gas migration systems beneath the Omakere Ridge, offshore New Zealand. The new method is simple and should be easy to implement to enhance seismic interpretation workflows.

Summary This work presents a novel technology for azimuth-dependent facies analysis (Facies Analysis versus Azimuth — FACIVAZ) to improve the prediction of hydrocarbon-saturated permeable fractures in terrigenous carbonate reservoirs. The analysis is performed in the depth domain along high-resolution, full-azimuth, angle domain common image gathers created by the EarthStudy 360™ Local Angle Domain (LAD) imaging system. The amplitude and phase preservation of the seismic reflectivities obtained by this imaging system is crucial to the proposed analysis. Prior to the facies analysis, the general orientation and intensity of the target fracture systems are analysed and characterized by azimuth-dependent velocity and amplitude analyses (VVAZ and AVAZ) performed along these LAD gathers. The remaining effects of the azimuth-dependent (and frequency-dependent) absorption and dispersion on the LAD gather events are then detected and further connected to the rate of the existing oil-saturated fractures within the reservoirs. The examples presented in this article show the effectiveness of the proposed FACIVAZ technology in accurately predicting the distribution of seismic facies in target production areas associated with oil-saturated fractured reservoirs in Western Siberia and Middle Volga. The results strongly agree with the corresponding facies characteristics measured in the boreholes along the reservoir area, and therefore serve as valuable information for the drilling decisions of new wells.

SUMMARY The inverse problem of finding the slowness vector from a known ray direction in general anisotropic elastic media is a challenging task, needed in many wave/ray-based methods, in particular, solving two-point ray bending problems. The conventional resolving equation set for general (triclinic) anisotropy consists of two fifth-degree polynomials and a sixth-degree polynomial, resulting in a single physical solution for quasi-compressional (qP) waves and up to 18 physical solutions for quasi-shear waves (qS). For polar anisotropy (transverse isotropy with a tilted symmetry axis), the resolving equations are formulated for the slowness vectors of the coupled qP and qSV waves (quasi-shear waves polarized in the axial symmetry plane), and independently for the decoupled pure shear waves polarized in the normal (to the axis) isotropic plane (SH). The novelty of our approach is the introduction of the geometric constraint that holds for any wave mode in polar anisotropic media: The three vectors—the slowness, ray velocity and medium symmetry axis—are coplanar. Thus, the slowness vector (to be found) can be presented as a linear combination of two unit-length vectors: the polar axis and the ray velocity directions, with two unknown scalar coefficients. The axial energy propagation is considered as a limit case. The problem is formulated as a set of two polynomial equations describing: (i) the collinearity of the slowness-related Hamiltonian gradient and the ray velocity direction (third-order polynomial equation) and (ii) the vanishing Hamiltonian (fourth-order polynomial equation). Such a system has up to twelve real and complex-conjugate solutions, which appear in pairs of the opposite slowness directions. The common additional constraint, that the angle between the slowness and ray directions does not exceed ${90^{\rm{o}}}$, cuts off one half of the solutions. We rearrange the two bivariate polynomial equations and the above-mentioned constraint as a single univariate polynomial equation of degree six for qP and qSV waves, where the unknown parameter is the phase angle between the slowness vector and the medium symmetry axis. The slowness magnitude is then computed from the quadratic Christoffel equation, with a clear separation of compressional and shear roots. The final set of slowness solutions consists of a unique real solution for qP wave and one or three real solutions for qSV (due to possible triplications). The indication for a qSV triplication is a negative discriminant of the sixth-order polynomial equation, and this discriminant is computed and analysed directly in the ray-angle domain. The roots of the governing univariate sixth-order polynomial are computed as eigenvalues of its companion matrix. The slowness of the SH wave is obtained from a separate equation with a unique analytic solution. We first present the resolving equation using the stiffness components, and then show its equivalent forms with the well-known parametrizations: Thomsen, Alkhalifah and ‘weak-anisotropy’. For the Thomsen and Alkhalifah forms, we also consider the (essentially simplified) acoustic approximation for qP waves governed by the quartic polynomials. The proposed method is coordinate-free and can be applied directly in the global Cartesian frame. Numerical examples demonstrate the advantages of the method.

This paper is the second in a sequel of two papers and dedicated to the computation of paraxial rays and dynamic characteristics along the stationary rays obtained in the first paper. We start by formulating the linear, second‐order, Jacobi dynamic ray tracing equation. We then apply a similar finite‐element solver, as used for the kinematic ray tracing, to compute the dynamic characteristics between the source and any point along the ray. The dynamic characteristics in our study include the relative geometric spreading and the phase correction due to caustics (i.e. the amplitude and the phase of the asymptotic form of the Green's function for waves propagating in 3D heterogeneous general anisotropic elastic media). The basic solution of the Jacobi equation is a shift vector of a paraxial ray in the plane normal to the ray direction at each point along the central ray. A general paraxial ray is defined by a linear combination of up to four basic vector solutions, each corresponds to specific initial conditions related to the ray coordinates at the source. We define the four basic solutions with two pairs of initial condition sets: point–source and plane‐wave. For the proposed point–source ray coordinates and initial conditions, we derive the ray Jacobian and relate it to the relative geometric spreading for general anisotropy. Finally, we introduce a new dynamic parameter, similar to the endpoint complexity factor, presented in the first paper, used to define the measure of complexity of the propagated wave/ray phenomena. The new weighted propagation complexity accounts for the normalized relative geometric spreading not only at the receiver point, but along the whole stationary ray path. We propose a criterion based on this parameter as a qualifying factor associated with the given ray solution. To demonstrate the implementation of the proposed method, we use several isotropic and anisotropic benchmark models. For all the examples, we first compute the stationary ray paths, and then compute the geometric spreading and analyse these trajectories for possible caustics. Our primary aim is to emphasize the advantages, transparency and simplicity of the proposed approach.

Advanced seismic data and multi-attribute visualization techniques such as the red-green-blue (RGB) have considerably augmented the capability of interpreters to characterize geological features on three-dimensional (3D) seismic reflection datasets. However, high resolution investigation of complex features remains challenging, requiring additional approaches to improve all round interpreters’ performance. Intervals of interest are commonly associated with intricate geological features, including channel systems and fluid migration pathways, which may be concealed in the dataset as multilevel high-amplitude seismic reflectivity bands. Delineating such features onto a single, intuitive map is often an arduous task. This may prove even more demanding where such features are spatially connected and occur near discontinuous and difficult-to-interpret horizons. To aid this task, we have developed an innovative technique involving RGB-blending and composition of multilevel amplitude maps, multilevel composition. Multilevel composition involves identification of high-amplitude features of geological interest within the dataset and defining their window of occurrence . This is followed by the interpretation of at least one reflecting horizon within/around the target window and the division of the window into three sub-windows based on the spatiotemporal distribution of the geological features. Amplitude-accentuating seismic attributes are computed for the sub-windows, the resulting maps are assigned to red (shallowest level), green (intermediate level) and blue (deepest level) colors and are co-rendered in the RGB color space. This results in a single map in which inter-window/layer depth information is coded in colors for reliable representation of the actual geology . We demonstrate the efficacy of this technique by applying it to characterize classic deep-water depositional elements in the eastern Nile fan, eastern Mediterranean, and to investigate seafloor seeps and underlying fluid plumbing systems beneath the Omakere Ridge, offshore New Zealand, using high-resolution 3D seismic data. The new technique is simple and easy to execute and enhances existing seismic interpretation workflows.

A seismic survey was conducted in a production oilfield located in Serbia, confined to the Pre-Neogene (Paleozoic) base of the Pannonian basin. It was assumed that additional significant unrecovered residual reserves still exist in this oil field, as well as additional similar undiscovered reservoirs. The further characterization of the existing reservoirs, and the identification and characterization of the new ones, required the implementation of advanced seismic imaging technology. Hence, a new project was designed composed of the following steps: Obtaining the highest possible seismic resolution in the area, and creating an updated, high-definition subsurface model that includes the structural complexities of the geological layers and the azimuthal anisotropic effects (e.g., fracture systems), especially within the target layers. This enables the identification and characterization of the target productive zones, making it possible to accurately design and plan the well placement. A modern, full-azimuth seismic survey was performed with a fairly regular distribution of the source-receiver offsets and azimuths. The dominant fold (number of traces per shot) was about 120. The seismic sources consisted of groups of vibrators with a linear sweep signal, where the frequency range was 6 to 96 Hz and the time duration 15 seconds. Emerson’s EarthStudy 360 full-azimuth angle domain imaging system (Koren and Ravve, 2011) was chosen to facilitate the above-mentioned tasks. This is an advanced subsurface imaging system operating directly in the Local Angle Domain (LAD). The high-resolution images with the unique full-azimuth directional and reflection angle common image gathers obtained by this seismic migration technology make it possible to better define the structural subsurface model and furthermore, to detect fine interlayer fracture systems at the target areas. Both regional faults and low-amplitude sub-seismic faults (fracture indicators at these regions) were mapped. The main imaging characteristics were correlated with existing production wells in the area

Kinematical characteristics of reflected waves in anisotropic elastic media play an important role in the seismic imaging workflow. Considering compressional and converted waves, we derive new, azimuthally dependent, slowness‐domain approximations for the kinematical characteristics of reflected waves (radial and transverse offsets, intercept time and traveltime) for layered orthorhombic media with varying azimuth of the vertical symmetry planes. The proposed method can be considered an extension of the well‐known ‘generalized moveout approximation’ in the slowness domain, from azimuthally isotropic to azimuthally anisotropic models. For each slowness azimuth, the approximations hold for a wide angle range, combining power series coefficients in the vicinity of both the normal‐incidence ray and an additional wide‐angle ray. We consider two cases for the wide‐angle ray: a ‘critical slowness match’ and a ‘pre‐critical slowness match’ studied in Parts I and II of this work, respectively. For the critical slowness match, the approximations are valid within the entire slowness range, up to the critical slowness. For the ‘pre‐critical slowness match’, the approximations are valid only within the bounded slowness range; however, the accuracy within the defined range is higher. The critical slowness match is particularly effective when the subsurface model includes a dominant high‐velocity layer where, for nearly critical slowness values, the propagation in this layer is almost horizontal. Comparing the approximated kinematical characteristics with those computed by numerical ray tracing, we demonstrate high accuracy.

Predicting the permeability of fractured reservoirs is valuable for both reservoir assessment and drilling planning. Characterization of such systems requires advanced amplitude analysis, mainly based on seismic imaging results of the recorded wavefield. A significant amount of work has been done on the subject over the past few decades. Thomsen (1995) showed the effect of seismic amplitude variations for different fractured media. Pisetsky and Fedorov (1998) showed the influence of the size of cracks and the seismic wave length on the change of reflection coefficient. Ruger (1998) used the Zoeppritz equation to formalize the reflection coefficient variation as a function of the wavefront azimuthal direction with respect to the fracture azimuth. An approximation of the latter by Tsvankin and Grechka (1998) became the basis for AVAz inversion analysis. Studies have also been performed on frequency absorption properties of fluid-saturated and dry fractured layers. Goloshubin (2002) observed that reflections are stronger for fluid-saturated fractured layers than for dry fractured layers at low-frequency ranges. Geek (2008) studied high-frequency absorption vs. wavefront azimuth for dry fractures in the acoustic frequency range for P-waves, and showed that the absorption is insensitive to the azimuth for S-waves. Kozlov (2006) made a comprehensive generalization of the physics of frequency absorption effects in the case of dry and fluid-saturated fractures. In his study he showed that in fractured reservoirs, high fluid permeability has an observable effect of absorption of the low-frequency spectral range for reflected P-waves. This observation was also noted in several exploration areas of Western Siberia (Upper and Middle Jurassic deposits, Davydova, 2004). Despite the research detailed above, predicting and delineating oil-saturated fracture systems using low-frequency absorption analysis has not yet become a standard practice in the industry. This is in large part due to the approximation made by standard processing and imaging methods in creating azimuthal data. To make this process more effective, a true amplitude full-azimuth imaging procedure that decomposes the fully recorded seismic wavefield in the local angle domain is required. This procedure is carried out in-situ, in the depth migrated domain, preserving both reflection azimuth and structural azimuth. Such a technology is the ideal solution for extracting the main faults, building the Structural Tectonic Skeleton (STS), exposing the main fracture orientation, and eventually measuring anisotropic frequency absorption effects. This paper presents a study of carbonate reservoirs in an oilfield in the Middle Volga region of Russia, and suggests a workflow based on imaging and processing in the Local Angle Domain to predict prospective areas of oil-saturated permeable fractured reservoirs.

SUMMARY Considering general anisotropic (triclinic) media and both, quasi-compressional (qP) and quasi-shear (qS) waves, in Part I of this study, we obtained the ray (group) velocity gradients and Hessians with respect to the ray locations, directions and the elastic model parameters along ray trajectories. Ray velocity derivatives for anisotropic elastic media with higher symmetries were considered particular cases of general anisotropy. In this part, Part II, we follow the computational workflow presented in Part I, formulating the ray velocity derivatives directly for polar anisotropic media (transverse isotropy with tilted axis of symmetry, TTI) for the coupled qP waves (quasi-compressional waves) and qSV waves (quasi-shear waves polarized in the ‘axial’ plane) and for SH waves (shear waves polarized in the ‘normal’ plane). The acoustic approximation for qP waves is considered a special case. In seismology, the medium properties, normally specified at regular 3-D fine gridpoints, are the five material parameters: the axial compressional and shear wave velocities, the three (unitless) Thomsen parameters and two geometric parameters: the polar angles defining the local direction (the tilt) of the medium symmetry axis. All the parameters are assumed spatially (smoothly) varying, so that their spatial gradients and Hessians can be reliably numerically computed. Two case examples are considered; the first represents compacted shale/sand rocks (with positive anellipticity) and the second, unconsolidated sand rocks with strong negative anellipticity (manifesting a qSV triplication). The ray velocity derivatives obtained in this part are first tested by comparing them with the corresponding numerical (finite difference) derivatives. Additionally, only for validation purpose, we show that exactly the same results (ray velocity derivatives) can be obtained if we transform the given polar anisotropic model parameters (five material and two geometric) into the 21 stiffness tensor components of a general anisotropic (triclinic) medium, and apply the theory derived in Part I. Since in many practical wave/ray-based applications in polar anisotropic media only the spatial derivatives of the axial compressional wave velocity are taken into account, we analyse the effect (sensitivity) of the spatial derivatives of the other parameters on the ray velocity and its derivatives (which, in turn, define the corresponding traveltime derivatives along the ray).

Summary Accuracy of initial velocity models has an impact on subsequent processing and use of seismic data, including curved-ray time migration, depth migration and tomography. Generating instantaneous velocities from stacking or rms velocity functions in vertical time involves an inversion step in which physical phenomena and geological characteristics of the subsurface can be taken into account by specific equations in a constrained Dix inversion method. In the neighbouring domain of geomodelling, interpolating sparse values with geological constraints on any type of discrete model, such as triangulated surfaces or volumetric grids, is a well-researched field. Building on that experience, the constrained Dix inversion method can be extended to integrate new types of constraints, thus improving the accuracy and plausibility of initial velocity models. In this paper, we first review the constrained Dix inversion method then show how constraints can be extended using a mathematical framework initially designed for geomodelling. Finally, we review the effects of adding these new types of constraints on inversion results.

Summary Considering 3D heterogeneous media and general anisotropy, and a parameterized stationary ray path obtained by the kinematic Eigenray, we analyse the second traveltime variation and derive the linear, second-order, vector-form Jacobi dynamic ray tracing (DRT) ODE. It delivers paraxial rays defined by shift vectors normal to the ray direction in its vicinity. The solution is obtained by applying the same finite element scheme (with the Hermite polynomial interpolation) used in the kinematic Eigenray. The resolving matrix of the linear algebraic DRT equation set coincides with the global traveltime Hessian already computed in the kinematic Eigenray stage. Two different solutions of the Jacobi DRT ODE with their corresponding point-source initial conditions, related to the chosen ray coordinates (RC), are needed to compute the ray Jacobian representing to the signed cross-area of the corresponding ray tube. For the chosen RC, we derive an original relationship between the ray Jacobian and the relative geometric spreading. The proposed Eigenray method has been tested using a benchmark numerical example with orthorhombic elliptic factorized inhomogeneous anisotropy. This model has an analytic solution for the ray path configuration, traveltime, arclength, and parameter sigma. The results demonstrate the high accuracy obtained by the proposed method.

Considering horizontally layered transversely isotropic media with vertical symmetry axis and all types of pure-mode and converted waves we present a new wide-angle series approximation for the kinematical characteristics of reflected waves: horizontal offset, intercept time, and total reflection traveltime as functions of horizontal slowness. The method is based on combining (gluing) both zero-offset and (large) finiteoffset series coefficients. The horizontal slowness is bounded by the critical value, characterised by nearly horizontal propagation within the layer with the highest horizontal velocity. The suggested approximation uses five parameters to approximate the offset, six parameters to approximate the intercept time or the traveltime, and seven parameters to approximate any two or all three kinematical characteristics. Overall, the method is very accurate for pure-mode compressional waves and shear waves polarised in the horizontal plane and for converted waves. The application of the method to pure-mode shear waves polarised in the vertical plane is limited due to cusps and triplications. To demonstrate the high accuracy of the method, we consider a synthetic, multi-layer model, and we plot the normalised errors with respect to numerical ray tracing.

Sedimentary layers affected by vertical compaction and strong lateral tectonic stresses are often characterized by low anisotropic symmetry (e.g., tilted orthorhombic [TOR]/monoclinic or even triclinic). Considering all types of pure-mode and converted waves, we derive the normal moveout (NMO) series coefficients of near normal-incidence reflected waves in arbitrarily anisotropic horizontally layered media, for a leading error term of order six. The NMO series can be either a function of the invariant horizontal slowness (slowness domain) or the surface offset (offset domain). The NMO series coefficients, referred to also as effective parameters, are associated with the corresponding azimuthally varying NMO velocity functions. We distinguish between local (single-layer) and global (overburden multilayer) effective parameters, which are related by forward and inverse Dix-type transforms. We derive the local effective parameters for an arbitrary anisotropic (triclinic) layer, which is the main contribution of this paper. With some additional geologic constraints, the local effective parameters can then be converted into the interval elastic properties. To demonstrate the applicability of our method, we consider a synthetic layered model in which each layer is characterized with TOR symmetry. The corresponding global effective model loses the symmetries of the individual layers and is characterized by triclinic symmetry.

Considering all types of pure-mode and converted waves, we derive the azimuthally dependent, fourth-order normal moveout (NMO) velocity functions, and hence the corresponding effective anellipticity functions, for horizontally layered orthorhombic media. We emphasize that this paper does not suggest a new nonhyperbolic traveltime approximation; rather, it provides exact expressions of the NMO series coefficients, computed for normal-incidence rays, which can then be further used within known azimuthally dependent traveltime approximations for short to moderate offsets. We do not assume weak anisotropy or acoustic approximation for P-waves. At each layer, the elastic parameters, thickness, and azimuth of the orthorhombic vertical symmetry planes are considered to be different. We distinguish between two different azimuths: slowness azimuth (part 1 of this paper) and offset azimuth (part 2 of this paper). In part 1, the slowness-azimuth domain NMO is approximated as a series of either infinitesimal horizontal slowness (slowness-azimuth/slowness domain) or infinitesimal offsets (slowness-azimuth/offset domain). Similarly, in part 2, we distinguish between two offset-azimuth domains: offset-azimuth/slowness and offset-azimuth/offset. Note that the azimuthally dependent NMO velocity functions of each of the four cases are different. The validity of the method is tested by introducing our derived azimuthally dependent, fourth-order effective anellipticity, into the well-known azimuthally dependent, asymptotic nonhyperbolic traveltime approximation, in which we compare the traveltime approximation versus exact numerical ray tracing for short to moderate offsets. It is clearly shown that for these types of azimuthally anisotropic layered models, the fourth-order terms are essential even for relatively small horizontal-slowness values or short offsets.

Based on the theory derived in part 1, in which we obtained the azimuthally dependent fourth-order normal-moveout (NMO) velocity functions for layered orthorhombic media in the slowness-azimuth/slowness and the slowness-azimuth/offset domains, in part 2, we extend the theory to the offset-azimuth/slowness and offset-azimuth/offset domains. We reemphasize that this paper does not suggest a new nonhyperbolic traveltime approximation; rather, it provides exact expressions of the NMO series coefficients, computed for normal-incidence rays, which can then be further used within known azimuthally dependent traveltime approximations for short to moderate offsets. The same type of models as in part 1 are considered, in which the layers share a common horizontal plane of symmetry, but the azimuths of their vertical symmetry planes are different. The same eight local (single-layer) and global (overburden multilayer) effective parameters are used. In addition, we have developed an alternative set of global effective parameters in which the “anisotropic” effective parameters are normalized, classified into two groups: two “azimuthally isotropic” parameters and six “azimuthally anisotropic” parameters. These parameters have a clearer physical interpretation and they are suitable for inversion purposes because they can be controlled and constrained. Next, we propose a special case, referred to as “weak azimuthal anisotropy,” in which only the azimuthally anisotropic effective parameters are assumed to be weak. The resulting NMO velocity functions are considerably simplified, reduced to the form of the slowness-azimuth/slowness formula. We verify the correctness of our method by applying it to a multilayer orthorhombic medium with strong anisotropy. We introduce our derived, fourth-order slowness-azimuth/offset domain NMO velocity function into the well-known nonhyperbolic asymptotic traveltime approximation, and we compare the approximate traveltimes with exact traveltimes obtained by two-point ray tracing. The comparison shows an accurate match up to moderate offsets. Although the accuracy with the weak azimuthal anisotropic formula is inferior, it can still be considered reasonable for practical use.

The well-known asymptotic fractional four-parameter traveltime approximation and the five-parameter generalised traveltime approximation in stratified multi-layer transversely isotropic elastic media with a vertical axis of symmetry have been widely used for pure-mode and converted waves. The first three parameters of these traveltime expansions are zero-offset traveltime, normal moveout velocity, and quartic coefficient, ensuring high accuracy of traveltimes at short offsets. The additional parameter within the four-parameter approximation is an effective horizontal velocity accounting for large offsets, which is important to avoid traveltime divergence at large offsets. The two additional parameters in the above-mentioned five-parameter approximation ensure higher accuracy up to a given large finite offset with an exact match at this offset. In this paper, we propose two alternative five-parameter traveltime approximations, which can be considered extensions of the four-parameter approximation and an alternative to the five-parameter approximation previously mentioned. The first three short-offset parameters are the same as before, but the two additional long-offset parameters are different and have specific physical meaning. One of them describes the propagation in the high-velocity layer of the overburden (nearly horizontal propagation in the case of very large offsets), and the other characterises the intercept time corresponding to the critical slowness that includes contributions of the lower velocity layers only. Unlike the above-mentioned approximations, both of the proposed traveltime approximations converge to the theoretical (asymptotic) linear traveltime at the limit case of very large (“infinite”) offsets. Their accuracy for moderate to very large offsets, for quasi-compressional waves, converted waves, and shear waves polarised in the horizontal plane, is extremely high in cases where the overburden model contains at least one layer with a dominant higher velocity compared with the other layers. We consider the implementation of the proposed traveltime approximations in all classes of problems in which the above-mentioned approximations are used, such as reflection and diffraction analysis and imaging.

Summary We present an efficient and stable procedure for estimating second- and fourth-order azimuthally-dependent effective parameters from full-azimuth residual moveouts. The residual moveouts are automatically picked at depth image points along full-azimuth angle domain reflection angle gathers. It is assumed that the azimuthally varying residual moveouts are due to fracture systems within compacted sand/shale sediment layers which were not accounted for in the seismic migration. The extracted (up to eight) effective parameters can then be used to obtain local (layer) effective parameters, characterizing the intensity and orientation of the fracture systems at each layer. Finally, the local effective parameters can be inverted to obtain interval anisotropic (e.g., orthorhombic) model parameters to be used in orthorhombic seismic migration.

We study the azimuthally dependent hyperbolic moveout approximation for small angles (or offsets) for quasi-compressional, quasi-shear, and converted waves in onedimensional multi-layer orthorhombic media. The vertical orthorhombic axis is the same for all layers, but the azimuthal orientation of the horizontal orthorhombic axes at each layer may be different. By starting with the known equation for normal moveout velocity with respect to the surface-offset azimuth and applying our derived relationship between the surface-offset azimuth and phase-velocity azimuth, we obtain the normal moveout velocity versus the phase-velocity azimuth. As the surface offset/azimuth moveout dependence is required for analysing azimuthally dependent moveout parameters directly from time-domain rich azimuth gathers, our phase angle/azimuth formulas are required for analysing azimuthally dependent residual moveout along the migrated local-angle-domain common image gathers. The angle and azimuth parameters of the local-angle-domain gathers represent the opening angle between the incidence and reflection slowness vectors and the azimuth of the phase velocity ψphs at the image points in the specular direction. Our derivation of the effective velocity parameters for a multi-layer structure is based on the fact that, for a one-dimensional model assumption, the horizontal slowness ph and the azimuth of the phase velocity ψphs remain constant along the entire ray (wave) path. We introduce a special set of auxiliary parameters that allow us to establish equivalent effective model parameters in a simple summation manner. We then transform this set of parameters into three widely used effective parameters: fast and slow normal moveout velocities and azimuth of the slow one. For completeness, we show that these three effective normal moveout velocity parameters can be equivalently obtained in both surface-offset azimuth and phase-velocity azimuth domains.

Hard-to-recover hydrocarbon reservoirs lie at great depths, in complex geological conditions. They are characterized by complex structures, low fluid permeability and a low oil and gas recovery ratio. In Russia today, about 60% of the potential oil and gas fields are located in this type of reservoir. These include hydrocarbon deposits in the Paleozoic basement (pre-Jurassic basement) of Western Siberia, subsalt carbonate sediments under salt dome tectonics, and carbonate and terrigenous deposits in the Volga region and Eastern Siberia. The exploration of these reservoirs benefits from a new, full-azimuth angle domain approach to seismic processing and imaging. This new technology can provide a more detailed depth image of the entire structural-tectonic reservoir skeleton, and a more accurate forecast of the main rock properties of the reservoirs. Conventional seismic depth imaging tools, such as ray-based or beam-based Kirchhoff migrations, applied to rich azimuth seismic data, normally generate multi-azimuth offset-domain common image gathers (CIGs). These are further used for anisotropic velocity model determination and for the characterization of reservoir properties, such as fracture systems. In these types of migrations, the input data is first binned into specific surface offset/ azimuth geometrical groups, such as offset vector tiles (OVT), azimuthal sectors or planner spirals, depending on the acquisition pattern. Each set of binned data is then independently migrated, with the final CIGs being simply a collection of individual images. However, in many cases, particularly when studying hydrocarbon reservoirs below complex geological areas and along steep inclined layers, the offset/azimuth CIGs do not provide the required information (in terms of accuracy and resolution, for example) to achieve the above mentioned goals. Unlike subsurface imaged events along the angle domain CIGs, which indicate ‘true’ local reflectivities, the reflection image events along the offset domain CIGs can be only considered a rough approximation of the ‘true’ reflectivities. Obviously, the accuracy and reliability of the offset domain CIGs are strongly compromised when imaging below complex geological areas with complex wave phenomena. One of the main drawbacks of offset domain imaging, especially in complex geological areas, is its inability to deal with the actual multi-pathing waves which are naturally handled within angle domain imaging.

Summary Considering reciprocity where the traveltime is an even function of the offset or horizontal-slowness, the fourth-order normal moveout (NMO) series are governed by the normal-incidence time and eight effective parameters: three second-order and five fourth-order. Local effective parameters are related to the individual layers, while the global effective parameters are related to the overburden multi-layer model. Local and global parameters are related by forward and inverse Dix-type transforms. The NMO formulae are different in the slowness and offset domains, but the eight parameters are the same in both cases. We suggest a new set of intuitive normalized effective parameters, classified into two “azimuthally isotropic” and six “azimuthally anisotropic” parameters. We provide feasible ranges for the normalized parameters, thus allowing their used for controlled inversion.

Summary Considering all types of pure-mode and converted waves, we derive the azimuthally-dependent normal moveout (NMO) series coefficients of near normal-incidence reflection waves in general anisotropic (triclinic) horizontally layered media, for a leading error term of order six. The NMO series can be either a function of the invariant horizontal-slowness (slowness domain) or the surface-offset (offset domain). The NMO series coefficients of different orders, also referred to as effective parameters, are associated with the corresponding azimuthally-dependent NMO velocity functions. We distinguish between local (single-layer) and global (overburden multilayer) effective parameters, where the local and global effective parameters are related by forward and inverse Dix-type transforms. We first consider the case in which the reciprocity assertion for incidence and reflected waves holds, i.e. pure-mode waves for general anisotropic horizontally-layered media, and converted waves for anisotropic horizontally-layered models sharing a common horizontal symmetry plane. Considering reciprocity, the odd-power coefficients of the NMO series cancel, and the remaining coefficients are zero-offset time, three second-order and five fourth-order effective parameters. Next we consider converted waves in general anisotropic media, where reciprocity no longer holds. Twelve additional parameters are required: two firstorder, four third-order and six fifth-order effective parameters.

The main goal of the seismic surveys conducted in the target area in Kazakhstan was to image and detect a major Devonian carbonate barrier reef and to characterize the density and orientation of its main fracture zones. This area is unique in that the carbonate layers occur at great depths, from 5700 to 8500 m, with a complex structure of overburden layers of interbedded clay, sandstone, salt and others (Table 1). These different lithology plays and morphology rocks create strong vertical and lateral velocity variations, resulting in a complex seismic wave phenomenon. In addition, the target carbonate structures contain heterogeneous objects and require high-quality processing of the recorded data. Under such conditions, it is crucial to use full-azimuth, long-offset and dense (high fold) acquisition patterns. Advanced processing sequence tools and high-end depth migrations are required to handle both strong heterogeneity and azimuthal anisotropy effects. Traditional Kirchhoff migrations, even the most accurate ones (wavefront reconstruction, beam), have not been able to provide the required image quality and level of detail required at the target zones (Figure 1). Conventional Kirchhoff migrations generate surface offset-azimuth/offset domain common image gathers (CIG). In this particularly complex area, the correlation between the surface offset-azimuths and the actual, in situ, subsurface slowness-azimuths (azimuth of the incidence/reflected ray pairs at the image points) is relatively poor, leading to significant errors in the estimation of fracture orientation. Additionally, Kirchhoff migrations do not account for multi-pathing (multi-wave path solutions between image points and surface source/receiver locations), which is essential when imaging in areas involving complex wave phenomena. Kirchhoff beam migrations map surface beams with given directions backward to the subsurface. Since they account for the multi-pathing, they provide better results. However, they normally generate the same type of surface offset-azimuth/offset CIGs and therefore cannot be accurately used for azimuthal studies. They also cannot ensure sufficient subsurface illumination, especially at the complex target subsurface regions.

Small-scale subsurface features, such as natural fractures, act as scattering sources for seismic waves propagating through the subsurface. The wavefield generated by those source points is identified as diffraction energy. The amplitude of this type of energy is much smaller than the recorded events reflected from actual interfaces between different geological layers. Moreover, diffraction energy is normally suppressed by conventional processing and standard imaging algorithms, where summations and averaging processes are applied. The common objective in such processing workflows is to focus on the high specular amplitudes in order to enhance the continuity of seismic reflection events for improving the structural mapping of the subsurface. Our goal is to complement the traditional seismic interpretation workflow by integrating information relative to diffraction energy as another seismic attribute to be interpreted. The technique applied in this paper is based on a depth imaging algorithm that maps and bins the recorded surface information into multi-dimensional, local angle domain (LAD) common image gathers. The advantage of this system is its unique ability to decompose the wavefield into reflection and diffraction energy directly at the image locations. This paper provides a brief overview of the technology and illustrates its benefits when applied to the Eagle Ford and Barnett shale reservoirs, where seismic data can be of moderate quality, leading to accurate, high-resolution, and highcertainty seismic interpretation for risk-managed field development.

We consider a case where a 3D depth migration has been performed in the local angle domain (LAD) using rich-azimuth seismic data (e.g., conventional land surveys). The subsurface geologic model is characterized by considerable azimuthally anisotropic velocity variations. The background velocity field used for the migration can consist of azimuthally independent, e.g., vertical transverse isotropy, and/or azimuthally dependent (e.g., orthorhombic), velocity layers. The resulting 3D full-azimuth reflection angle gathers generated by the LAD migration represent in situ high-resolution amplitude preserved reflectivities associated with opening angles between incident and reflected slowness vectors in the specular directions. Residual moveouts (RMOs) automatically picked on these 3D image gathers along major horizons can indicate considerable residual periodic azimuthal variations. This situation is typical in depth imaging applied to unconventional shale plays, where the background velocity model doesn’t yet account for the aligned stress/fracture systems that exist in some of the target layers. We use the azimuthally dependent, phase-angle RMOs to update the interval parameters of the background model, accounting for the azimuthal anisotropy effect. Until now, this problem was mainly treated in the unmigrated time-offset domain, which is limited in describing the actual in situ changes of the velocity field with azimuths. The subsurface full-azimuth phase-angle domain RMOs provide better physical parameters to analyze the in situ azimuthal variations of the anisotropic media. Our method is grounded in a newly derived generalized Dix-based theory, where locally the background and updated models are assumed to be 1D anisotropic velocity models. At each lateral location, the orthorhombic axis [Formula: see text] points in the vertical direction across all layers, but the azimuthal orientations of the orthorhombic layers change from layer to layer. An effective model for such a layered structure (background or updated) is represented by a single layer with a vertical time identical to that of the whole package, effective fast and slow normal moveout (NMO) velocities, and an effective azimuthal orientation of the slow NMO velocity. Our approach begins with computation of these effective parameters for the background model and conversion of the high-resolution RMOs into a dense set of updated, effective, azimuthally dependent NMO velocities, which are then converted into three effective parameters of the updated model. Next, we apply a generalized Dix-based inversion approach to estimate the local NMO parameters for each updated layer. Finally, we convert the local parameters into interval azimuthally varying anisotropic model parameters (e.g., TTI, orthorhombic, or tilted orthorhombic) within each layer. The 1D Dix-based approach presented in this work should not be considered an alternative to more accurate 3D global inversion approaches, such as global anisotropic tomography. However, the proposed method can be effectively used for moderately laterally varying models, and some of the principal physical rules derived for the 1D model can be further used to improve the formulation and geophysical constraints applied to 3D global inversion methods.

SUMMARY We present a method to improve imaging in the local angle domain (LAD) decomposition and imaging system. This system uses the entire recorded data to generate trueamplitude, angle-dependent or angle and azimuth dependent imaging gathers (Koren and Ravve 2011). These gathers have the ability to distinguish the wavefields by their directional components: Specular (continuous structural surfaces) and diffraction (discontinuous objects such as small-scale fractures and faults). The high-energy values associated with the specular directions can be used to enhance the continuous objects to obtain a diffractionfree, sharpened image of highly complex areas. We propose that the specular enhancement in the LAD system be used to re-evaluate existing land and marine (including narrow-azimuth legacy) seismic data to obtain more detailed high-resolution images without the need to acquire additional 3D data about existing assets.