Seismic Modeling of Complex Geological Structures

Seismic forward modeling is seismic forward realization of a given geological model (Carcione et al., 2002; Fagin, 1991; Krebes, 2004; Sayers & Chopra, 2009). Two main stages of seismic modeling are geological model building, and numerical computation of seismic response for the model. It describes the forward process of propagating waves from sources to scatterers down in the subsurface and back to the receivers. The quality of the computed seismic response is partly related to the type of model that is built. Therefore the model building approaches become equally important as seismic forward realization methods. Models are considered to be representations of real objects (Ellison, 1993) and can be 1D, 2D, or 3D. 1D models are usually generated at well locations to predict the seismic response of the geological model and further to investigate the link between the geological beds at the well to the real reflection seismic data (seismic to well tie analysis).


Introduction
Seismic forward modeling is seismic forward realization of a given geological model (Carcione et al., 2002;Fagin, 1991;Krebes, 2004;Sayers & Chopra, 2009).Two main stages of seismic modeling are geological model building, and numerical computation of seismic response for the model.It describes the forward process of propagating waves from sources to scatterers down in the subsurface and back to the receivers.The quality of the computed seismic response is partly related to the type of model that is built.Therefore the model building approaches become equally important as seismic forward realization methods.Models are considered to be representations of real objects (Ellison, 1993) and can be 1D, 2D, or 3D.1D models are usually generated at well locations to predict the seismic response of the geological model and further to investigate the link between the geological beds at the well to the real reflection seismic data (seismic to well tie analysis).
The increasing amount of data which new technologies (such as advanced multi-component 3D seismic surveys) are able to provide, together with the development of more powerful and numerically efficient computing systems, have led to the rapid growth of subsurface modeling techniques (e.g.Alaei & Petersen, 2007;Mallet, 2008).Model building techniques developed significantly over the past decades.Khattri & Gir (1976) used a series of lithological elements through a cyclic succession (for example sand and shale) to create different 1D seismic models.The seismic response of such models have been predicted using ray theory approach.May & Hron (1978) carried out zero offset ray tracing for primary P waves to predict seismic response of a series of simple 2D geological models including stratigraphic wedge, unconformity, anticline, reef, normal fault, growth fault, thrust fault, salt dome flank, and overhang salt dome.The 2D models consisted of homogenous layers separated with curved interfaces.Gjøystdal et al (1985) introduced solid modeling technique to build 3D models of complex geological structures.The model consisted of a series of columns or solids and the properties such as P and S wave velocities and density varied continuously within solids and discontinuously across model interfaces or boundaries.They have used the 3D models to run dynamic ray tracing.Open model building technique (Åstebøl, 1994 as cited in Vinje et al., 1999) unlike the solid modeling technique may contain holes and cracks in interfaces.A wide range of computer aided design (CAD) methods have been developed to build complex geological models.Mallet (2008) gave a thorough review of theses methods.Patel & McMechan (2002) provided an algorithm to create 2D geological models from controlled horizons and well log data.
The main goal of this chapter is to introduce seismic forward modeling as a powerful tool to investigate the seismic wave propagation in different geological settings with a special reference to complex geological structures.The source of complexities of seismic wave transmission and reflection in subsurface will be explained.Different model building approaches will be described with examples.Three different seismic forward realizations including asymptotic (ray tracing methods), integral, and direct (e.g.finite difference algorithms) methods will be presented.

Sources of complexities of seismic wave propagation
Most of the problems in seismic wave propagation of geological settings are due to the complexities in structure (structure dependent complexities) and rock types (structure independent or stratigraphic complexities).The term 'complex' is used for those geological settings which cannot be easily imaged (Fagin, 1991) due to special characteristics of structural or stratigraphic complexities.Examples of structural complexities are: steep dipping beds, faults with steep dips, complex faulted folding, folds with complex geometry, closely spaced folds and faults.Complex faulted/folded salt basins are good examples of complex geological settings.Near surface problems add more complexity to the seismic wave propagation in particular to the land seismic data with variable topography.Some of the near surface problems are: i) seismic data distortion due to near surface velocity variations, ii) topographic variations, iii) irregular data coverage caused by rugged topography, and iv) illumination problem caused by near surface complex velocity fields.In seismology, illumination is the amount of seismic wave energy that falls on a reflector and thus available to be reflected (Sheriff, 2004).The complication is caused by propagation of body waves through the complex near surface layers and source generated noise that are trapped in the near surface (Al-Ali & Verschuur, 2006).

Near surface problems
A near surface, low velocity layer (LVL) causes delay of seismic travel times.The term low velocity layer is often used for material above water table or to geologically unconsolidated deposits on harder consolidated rocks (Cox, 1999;Marsden, 1993).This seismic weathering layer despite its terminology is different from geologic weathered layer.The variability in both thickness and velocity of the near surface layers is the main source of problem.The LVL is usually above the water table and the pore spaces of rocks are filled with air rather than water which considerably lower the seismic velocity.Corrections must be applied to seismic travel time to compensate for the delay caused by the LVL.These corrections are part of the static corrections applied to seismic data and there are several methods available such as up-hole based statics, and first break statics.The main assumption behind the conventional static corrections is that raypath through a relatively simple near surface is almost vertical and therefore a vertical time shift can be used to reference the acquired data to a flat datum (Cox, 1999).
Topographic variation is one of the most complicating factors affecting reflection seismic imaging.Vertical near surface propagation assumption explained earlier may not be valid in case of rugged acquisition topography.The angular dependence of statics should be considered otherwise the diffractions would not be handled correctly in the subsequent imaging steps.Most of the conventional imaging methods require data has to be collected on a level or datum and with regular grid.In the case of land seismic, data are acquired along irregularly-sampled surface with varying topography.Redatuming with static shift can be used to remove the topographic variations.The objective is to determine the reflection events arrival times which, would have been observed if all recording were made on a flat datum.The limitation of conventional static corrections is known from before (e.g.Shtivelman & Canning, 1988) and alternative methods such as wave-equation datuming has been used instead (Al-Ali & Verschuur, 2006;Bevc, 1997;Reshef, 1991).Fig. 1 (Yang et al., 2009) shows a Prestack Depth Migrated (PSDM) seismic image from the Chinese Foothills using conventional static and wave equation based datuming.This example illustrates that it is necessary to compensate for the effect of complex propagation.Fig. 1.PSDM seismic image from Chinese Foothills after static (left) and after wave equation datuming (right).The image qulaity improved using wave equation datuming in particluar in the deeper section (Yang et al., 2009) The complexity of near surface can also cause poor illumination of deeper targets.Seismic wave propagation in near surface beds composed of incompetent rock types such as gypsum is complicated due to the sever heterogeneity of the rocks.Internal faulting and folding of such layers will add to the complexity of the wave propagation (Alaei & Pajchel, 2006).Fig. 2 shows a Prestack Time Migrated (PSTM) section from Zagros fold belt (Alaei, 2006).Incompetent material exposed at the surface of the line cause significant illumination problem for the deeper targets.Unusually high velocities at the near surface can also cause illumination problems.An example for that could be high velocity limestone near the sea floor in the Norwegian Barents Sea that act as a strong scatterer and complicate the wave propagation.

Subsurface problems
Subsurface complexities also complicate the wave propagation and vary depending on the rock type and dominant structural patterns.Among the various geological settings, saltrelated structures and structures of fold and thrust belts cause greater challenges for the propagation of seismic waves compared to other geologic settings.However, seismic modeling has been used to improve imaging in the salt basins and fold and thrust belts (Fagin, 1991).

Salt-related complexities:
Complex structure and strong velocity contrast of salt with sediments around in salt-related geological settings is a great challenge for most of the seismic imaging algorithms (Albertin et al., 2001;Ray et al., 2004;Seitchick et al., 2009).Signal to noise ratio is usually low in the vicinity of salt bodies in particular below the salt.Examples of such settings are Gulf of Mexico, Nordkapp Basin in the Norwegian Barents Sea, and Santos Basin offshore Brazil.Seismic wave propagation through such large velocity contrast and structural complexity is associated with many wave phases including primary reflections, diffractions, and diffracted reflections.Seismic modeling has been extensively used to plan accurate seismic acquisition surveys over complex salt related structures (Gjøystdal et al., 2007) and improve seismic processing flows (Aminzadeh et al., 1997;Gjøystdal et al., 2007;Huang et al., 2010).Fig. 3 shows a seismic image from the Nordkapp Basin, Norwegian Barents Sea.
Fold and thrust belt complexities: Fold and thrust belts (such as Zagros fold belt, Canadian Rocky Mountain, and Andean fold belt) are dominated by a series of thrust faults and steeply dipping rock units.Fold geometry, internal structure complexity, highly dipping layers, and faulting associated with folding complicate the wave propagation (Alaei, 2005;Lines et al., 2000).Reflection seismic images from fold and thrust belts have frequently failed to give the correct picture of the subsurface structures when tested by drilled wells (Lingrey, 1991).Due to velocity and structural complexity, rays are bent and there are nonhyperbolic arrival times in addition to the hyperbolic arrival times.Fig. 4 (Alaei, 2006) shows an example of the ambiguity in seismic images from the structures in the Zagros fold and thrust belt.Parts of the stratal geometry is clear while the central part (indicated by yellow circle) shows a lack of reflection signal.Different seismic imaging algorithms, acquisition designs, and interpreted geologic models of fold and thrust belts can be tested using seismic forward modeling technique.(Alaei, 2006) Fault shadows: Seismic wave propagation is complicated under fault planes (usually footwall zone) which cause an unreliable seismic image of the zone.This zone of poor illumination is called fault shadow.Seismic imaging algorithms that doesn't take into account lateral velocity variations above imaging points fail to provide correct image under fault planes.
Lateral lithology variations: Lithological variations within rock units can cause strong lateral velocity variations which can be associated with relatively simple structures.Seismic wave propagation is complex in such settings despite the relatively simple structure (Alaei, 2005).

Applications of seismic modeling
Seismic modeling is useful in a wide range of applications in exploration and earthquake seismology.It plays an important role in almost all aspects of exploration seismology such as seismic data acquisition, processing, interpretation, and reservoir characterization.It increases the reliability of seismic data analysis.
Applications in seismic acquisition: In seismic acquisition, seismic forward modeling reduces the risk in seismic exploration by providing quantitative information to design better 3D surveys (e.g.Gjøystdal et al., 2007;Laurain et al., 2004;Robertsson et al., 2007).In complex geological settings seismic forward modeling can be used to test different acquisition parameters and subsurface models to achieve the optimum data collection strategy.Illumination problems of target horizons have been addressed using 2D and 3D modeling studies.The results of illumination studies can be directly applied to survey layout design.There are two categories of illumination studies used for the feasibility purposes including global approach that provides information over the whole target interface and local approach that gives information at one point in time and space (Laurain et al., 2004).Subsalt imaging has been a challenge for exploration seismology for many years and the application of seismic modeling has considerably improved the acquisition survey design for subsalt imaging (e.g.Regone, 2007).The modeling studies showed that wide angle azimuth acquisition surveys provide better illumination from subsalt structures.Seismic modeling studies have been carried out to improve seismic data acquisition over complex geologic settings of fold and thrust belts (Alaei, 2005).Fig. 5 shows raypaths from one shot record of a complex faulted anticline setting from the Zagros fold thrust belt.The acquisition geometry includes an off-end source-receiver array.In off-end source-receiver array, the seismic source is at one side of the array and receivers are deployed at the other side of the array.Complex structural settings cause poor coverage of raypaths at deeper levels.The modeling shows that the large offset may partially improve the target illumination but there are still large areas of the subsurface with poor coverage.Although the 2D seismic surveys are the dominant acquisition pattern over mountainous terrains of complex geological settings, the example shown in fig. 5 clearly indicates that 2D seismic acquisition fails to provide good quality images from the subsurface and instead alternative methods such as 3D seismic acquisition can be used.However no single technology can improve the image as much as detailed analysis of survey parameters through seismic modeling.
Applications in seismic processing: Seismic modeling has been used to test different processing algorithms and flows.An important role of seismic modeling is to calibrate migration methods (Gray et al., 2001).
It can be used to optimize the processing sequences particularly in complex situations.Because of the important role of seismic modeling in seismic processing a number of synthetic models have been generated and widely used to test processing sequences.Some examples are the SEG/ EAGE 3D salt/overthrust model (Aminzadeh et al., 1997), Marmousi 2D model (Versteeg, 1994), The Society of Exploration Geophysicists Advanced Modeling Program (SEAM) (Pangman, 2007), Husky model (Stork et al., 1995;W.J. Wu et al., 1998) and Spratt Foothills model (Lines et al., 2001).Some of these models were used to test new imaging algorithms (e.g.R.S. Wu et al., 2008).Several seismic processing techniques such as multiple removal, velocity estimation (e.g.Chen & Du, 2010), migration (Moser & Howard, 2008), and seismic inversion (Jang et al., 2009;Hu et al., 2011) have been tested using the Marmousi synthetic data.SEG/EAGE 3D salt and overthrust models and associated synthetic seismic data have been used to test different migration velocity estimation and seismic imaging methods (e.g.Operto et al., 2000;Xu et al., 2004).The SEG/EAGE salt model is similar to the salt features of Gulf of Mexico and the overthrust model is similar to structures of thrust belts from South America.
Applications in seismic interpretation: Seismic forward modeling can be used to relate the response of an interpreted geologic model to real data.One application is the development of geological models to investigate the structural and stratigraphic problems faced during the seismic interpretation (Chopra & Sayers, 2009).It can be used to check the validity of interpretation particularly in complex situations.Seismic image data quality of complex geological settings is often poor that the reliance on structural styles in complex geological settings is necessary in view of the fact that the quality of the seismic images of such settings is poor.Parts of the stratal geometry maybe clearly shown while other parts show either a lack or a confusing overabundance of reflection signals (Lingrey, 1991).Seismic modeling can be used to investigate the validity of models representing different structural styles and find the best match with the real seismic data (Alaei, 2006;Alaei & Petersen, 2007;Lingrey, 1991;Morse et al., 1991).

Model building approaches
The integration of different data types for model definition in space and time is increasing.The model building methods can be divided into two categories: Interface based methods (e.g.Alaei, 2005) and grid based methods (e.g.Mallet, 2002).The model type can considerably influence the quality of the seismic realizations from the model.Fagin (1991) suggested a range of questions to avoid errors caused by constructing improper models.These questions are: Sideswipes (structural features that lies off the 2D profile) can not be simulated using 2D models.However, the 3D modeling has the capability of simulating sideswipes.Seismic response of 3D models can be viewed in different directions including time slice and mapped on the geological surfaces of the model.The model size depends on the modeling objectives.Target interface size, and depth are some of the main factors controlling the model size.It can be very large to study regional structural settings (e.g.Alaei, 2005) or small scale to investigate numerical simulations of petrophysical properties of rocks (e.g.Saenger et al., 2007).
The process starts with building the geometry of model and followed by propagating different properties such as velocity and density within different units of the model.Geometry of model is composed of stratigraphic surfaces (horizons) and faults irrespective of modeling approach.Examples of horizons are top and base of reservoir rocks, unconformities, top and base of salt, and surfaces that correspond to significant velocity variations.Faults are structural surfaces that juxtapose rocks of different properties and cause seismic wave scattering.These components shall be selected based on geological and modeling objectives.Modeling objectives that have to be included in addition to geological objectives are those which satisfy seismic wave propagation through the model.For example if there is significant velocity variation above target horizon (overburden), additional surfaces or interfaces should be included to properly simulate the seismic wave propagation through the variable velocity overburden.Layers representing velocity inversions such as thrust faults and base of salt bodies are important for modeling as they cause defocusing of seismic waves.

Interface based modeling
Interface based model building approach starts with defining the model dimensions and is followed by selecting horizons and faults of the model.The structure is constructed by interfaces (curves or lines).The curves are composed of points in depth or time domain.A minimum number of points are required to build an interface using an interpolation algorithm such as spiline method.Some of the seismic realization methods (e.g.Ray tracing methods) needs continuous second derivative.For Ray tracing methods variations in the interface geometry should be small relative to the dominant wavelength in the seismic signal.Curvature radius of the interfaces is an attribute that can be used to define a threshold for the interface smoothness.The minimum curvature radius of model interfaces should be larger than the dominant wavelength in the seismic signal.Curves representing horizons should usually be long enough to cross the model lateral boundaries.Horizons can either cross the model lateral boundaries or other horizons above or below (for example unconformities).The intersection of interfaces with either each other or model boundaries is necessary for defining blocks between the interfaces (solid model).The cross cutting horizons add to the complexity of the model.In complex models it is useful to start building the large scale architecture of the model first and then add more details into the model.The area bounded by interfaces (horizons or faults) and model boundaries corresponds to layers or blocks which include seismic properties (P and S wave velocities and density).Fig. 6 illustrates a model with 14 interfaces and corresponding blocks.There are two possibilities to define faults in an interface-based model (Fagin, 1991).They can be modeled as separate surfaces that cut the stratigraphic surfaces, or they can be represented as offsets in the modeled horizons.Although defining the fault planes as surfaces cutting the stratigraphic units is difficult, it allows us to follow the reflections from the fault plane.It is useful to provide information about velocity before we describe the sources of velocity data for the modeling purpose.Seismology in general and exploration seismology in particular is overflowing with velocities (Margrave, 2003).To name a few, there are instantaneous velocity, average velocity, interval velocity, root mean square (rms) velocity, migration velocity, stacking velocity, phase velocity, and group velocity.The type of velocity that is used for seismic forward modeling is the interval velocity which is simply derived by dividing the thickness of a particular layer by the travel time through the layer.
Sources of interval velocity are sonic wire line logs, checkshot surveys, and Vertical Seismic Profile (VSP) data.The thickness of the time intervals varies from 1 to 3 feet in sonic logs to hundreds of meters in checkshot surveys.Checkshot data provide travel times from source that is usually located at the land or sea surface to receivers located in the borehole and can be used to estimate interval velocity.VSP data acquired in the same way as checkshot data but includes closely (and usually evenly) spaced measurement points.VSP data can be considered as high resolution checkshot data that unlike the checkshot survey that uses only the first break data uses the entire trace information.
The interval velocity data can be derived from seismic prestack gathers.When the subsurface layers are horizontal and velocity varies more in vertical direction, reflections from interfaces are described by hyperbolas (e.g.Binodi, 2006).The change in receiver to source (offset) distance causes a delay in reflection arrival time known as moveout.For a multilayer subsurface the travel time at an offset x is (Dix, 1955): Where T 0 is the reflection travel time at zero offset and T x is the reflection travel time at offset x.The Vrms formula for n layers is: where V(t i ) is the interval velocity of layer i and t i is the time thickness of layer i.The denominator of the formula corresponds to the total two way travel time to the base of the nth layer.The interval velocity is the one that can be directly used in modeling studies.The equation 3 can be solved for interval velocity (Dix, 1955), The velocity that is measured from seismic gathers is moveout or stacking velocity (V NMO ) which under certain conditions (stratified flat isotropic settings) is equivalent to rms velocity.However in complex geological settings with dipping layers and lateral velocity variations velocities measured from the seismic gathers can not be directly used to estimate interval velocity through Dix equation.Levin (1971) provided the following equation to account for the dip using V NMO : cos( ) where θ is the dip angle.When the geological model is too complex and lateral velocity variations is too strong equation 3 can no longer provide accurate estimate of interval velocity and advanced model based methods must be used to estimate interval velocity.Reflection tomography is one of these methods that estimate interval velocity by using an inversion procedure to fit modeled travel times to measured travel times.Fig. 7 (Alaei, 2005) shows a regional 2D model from Zagros fold and thrust belt southwest of Iran that is 81 by 17km.The model is built using interface-based model building approach.

Grid or cell-based modeling
Constructing cell-based geological models has received a lot of attention in the past decades.
In the model building process of the subsurface, model elements including faults and horizons are modeled as triangulated surfaces (Mallet, 2002).In Discrete Model (Mallet, 2008) geological model is represented by a set of points called 'nodes' that are linked to their neighbors.The nodes together with the linked neighbors generate a gird.Each node of the grid is associated with both coordinates (x,y,z) and values of physical properties (such as velocity or density).Patel & McMechan (2003) used well log data and control horizons to build a gridded model from physical properties such as seismic velocity.Inverse distance weighting or linear interpolation has been used to extend the well log information into the 2D model.The geometry of the control horizons is used to control the spatial extent of the properties.To obtain data for building any model with this method it is required to provide sufficient wells to sample every element in the model and enough control horizons to define the lateral extent of the structures.
Petersen (1999) proposed a modeling approach -compound modeling -to construct geological models.The compound model is composed of compound cells and each cell occupies an area.Different physical property distributions are assigned to each cell.The properties can vary within each cell.The model is consistent with geological evolution since the final distribution of properties emulates geological processes over time.In complex structural models where the sequence of events is important this feature of compound modeling will make it possible to differentiate between different stages of deformation (faulting, folding or fault-related folding, erosion etc.).If, for example, the result of a geological process such as a deformation phase, orogenic event or a sedimentation-related process is overprinted by the result of another process, the properties belonging to the latest stage (in 'time') replace the previous one for a specific position in 'space'.This ensures the time and space consistency of the geological model.The property distributions involve ranking.So, in the case of several property functions for one position in space, the one with highest priority derived from ranking will be selected.The geometry is controlled by curves of parametric description and properties by 1D functions of depth.Some characteristics of the curves are: (1) made of isolated points (x, z); (2) continuous by spline interpolation (spline means that the curve is continuous to second order); (3) x and z are functions of a common parameter, so that the curve may take any shape.The property distribution in space relates directly to space by property cells.The characteristics of property cells are the curve, the property values and curve orientation.The compound model is transformed into a grid using corresponding setting parameters that have been applied for different compound cells of the model.The grid point positions along x and z can be set according to the model requirements for the grid realization (gridding).The internal geometry of geological units gives some information about the deposition and post-deposition history of the units.(Vail, 1987) Any modeling attempt without taking into consideration such details will not represent the real geology.Sedimentary bodies are rarely equi-dimensional, so proper modeling requires knowledge of the orientation of the geological units.Fig. 8 shows an example of the importance of internal orientation of geological units in modeling.If one just models the whole unit shown in fig.8 as one block without attention to internal structure and orientation, it will not represent the real situation.Therefore, a successful modeling approach is the one that can include such geometrical details in the model so that the output will be geologically consistent.It is possible with curves and the hierarchical approach in compound modeling to build any kind of internal geometry and orientation such as truncations, onlap, downlap, and complex small-scale faulting and folding inside the geological units.Fig. 9 (Alaei & Petersen, 2007) shows the regional 2D Zagros model shown in fig.7 that is constructed using Compound modeling approach.The model includes regional as well as small scale structural and stratigraphic details.The velocity model is based on the integration of different available data, including check shot data from 10 wells.All available density logs from the wells used in the model and for deeper layers constant values has been used.Fig. 9. 2D regional model (80x17km) from Zagros fold and thrust belt.It includes regional structural elements as well as small-scale stratigraphic detail.The color represents the scaled acoustic impedance (Alaei, 2006)

Seismic forward realizations
Seismic forward realizations can be carried out following the construction of model geometry and populating proper seismic properties.The goal is to predict seismic response of subsurface model recorded on a group of receivers.Seismic modeling methods can be classified into three main categories (Carcione et al., 2002): i)asymptotic, ii) integral-equation and iii) direct methods.

Asymptotic methods
Asymptotic methods (ray tracing methods) have been frequently used in seismic modeling and imaging.They do not take into account the full wavefield (e.g.Ćervenŷ, 2001).In these methods, the wavefield is considered as a series of certain events, with characteristic travel time and associated amplitude.Raypaths are traced either by solving a certain differential equation that can be extracted from seismic wave equation (girded models) or by using analytic results within layers and explicit Snell's law calculations (interface based models).Raypaths are unbent in a constant velocity layer, bend at velocity interfaces (in accordance with Snell's law), and reflect at an angle equal to incidence angle at impedance interfaces.Snell's law is the relation that governs the transmission and reflection of raypaths at velocity interfaces and is used to calculate the raypath bending at velocity interfaces, V i and V i+1 are the velocities at medium i and i+1 and α and θ are incidence and transmission angles respectively.Rays follow the geometrical rule of transmission/reflection (Snell law) also called geometric rays.Rays that follow the law of edge diffraction at a point are called diffracted rays (e.g.Klem-Musatov et al., 2008;Kravtsov & Orlov, 1993).
For a constant density variable velocity scalar wave equation 2 2 2 2 1 () an approximate high frequency solution can be written as P and T are functions of position and are smooth scalar functions.If we take the derivate of the equation and considering the high frequency assumption we get, These equations are basic equations in the asymptotic methods and are called eikonal and transport equation (e.g.Carcione et al., 2002;Ćervenŷ, 2001).The eikonal equation is a nonlinear partial differential equation of first order for arrival time, T. The transport equation represents a linear partial differential equation of first order in P(x).The eikonal equation describes the travel time behavior of seismic waves under high frequency condition (Bleistein et al., 2000).Kinematic ray tracing includes travel time computation of the rays and only requires seismic velocity of geological model while amplitude calculation (dynamic ray tracing) requires both velocity and density of the model to be defined.
There are different ray tracing modes depending on the source and receiver arrays (acquisition mode) and seismic modeling objectives which can be categorized into two main groups of zero offset and offset methods.Offset ray tracing includes a series of seismic traces recorded with different receivers but same source.Different source-receiver arrays can be used such as split spread (source in the middle of the receivers) and off-end (source at one side of receivers) arrays depending on the position of the source with respect to the receivers.This ray tracing mode simulates the same way that real seismic data are acquired and has been widely used to test different processing stages that are carried out on prestack common shot gather data.Processing steps such as dip filtering, and some prestack migration algorithms can be tested using the offset ray tracing of geological models (Fagin, 1991).Seismic response of a single source and receivers is called shot gather record.Shot gather ray tracing is one of the most common geometries used to simulate prestack seismic response from subsurface models.Shot gather ray tracing with off-end survey geometry was Such non hyperbolic event will not be properly imaged if poststack seismic imaging methods are applied (Biondi, 2006).
In the zero offset ray tracing modes there is no offset between the source and receiver.Two main zero offset ray tracing modes are normal incidence ray tracing and image ray tracing.
Normal incidence ray tracing is one of the methods used in the modeling of complex geological structures (e.g.Alaei, 2005).It simulate a Common Mid Point(CMP) stack section.CMP stack section is generated by combining (stacking) CMP gathers.CMP gathers are traces that would be recorded by a coincident source and receiver at each location and can be generated by resorting shot gather data.Stacking process attenuates random effects and improves the signal to noise ratio.However in case of complex geological settings with non hyperbolic moveout, stack section fail to give the correct image of the subsurface.
Fig. 11 shows raypaths from normal incidence ray tracing of part of the Zagros model (fig.7 x=30 to x=70km).The source/receiver spacing was 60m which is similar to real 2D data acquired from the same area.The objective is to get the stacked unmigrated image of the complex model.The raypaths are normal to model interfaces.The raypath distribution illustrates data density available from different parts of the model interfaces and can be used to identify defocused areas.Additional detailed seismic modeling can be applied to the defocused zones.The second zero offset ray tracing mode is image ray tracing.Image ray tracing is an approximation of migrated data that can be called simulated time migration (Hubral, 1977).
The raypath that represents the minimum on the time reflection surface emerges perpendicular to the recording surface.Rays are traced from points regularly distributed along the model top, and every time they hit an interface, two-way times are calculated.
Reflections are considered as a continuum of diffractions and each diffraction hyperbola collapses on its least travel time peak.It can also be used to locate reflections more accurately by converting time migrated data to depth along the image rays (Thorn, et al., 1986, Johansen, et al., 2007).Imaging steeply dipping beds of complex structures is a challenge to different seismic imaging methods.Synthetic seismic response of image ray tracing with 6km aperture.The arrow illustrates the location of steeply dipping reflector.Increasing aperture from 3 km (real data) to 6 km (synthetic data) improved the signal continuity (Alaei, 2006) fault plane.Image ray tracing with 60m source/receiver spacing applied to the model.Imaging the steeply dipping part of the thrust plane (ramp part) is particularly important.In order to improve seismic image from the complex thrust faulted structure, a range of different migration aperture have been tested.Migration aperture is the range of spatial data considered in seismic migration.Fig. 12 (right) shows seismic response of 6km aperture.The aperture used in the processing of the real data was about 3km.The modeling results illustrate that the aperture used in the processing of the real data was not sufficient to image the steep flank of the structure properly.Increasing aperture during migration of real data will decrease the signal to noise ratio which can be improved using post migration noise cancellation filters.

Integral-equation methods
The second group of seismic modeling methods are integral-equation methods.Integralequation methods of seismic modeling are based on integral representation of the seismic wavefield spreading from point sources (Huygens principal).There are two forms of integral methods: volume integral and boundary integral.Integral representation of scalar seismic wave equation is (Carcione et al., 2002), 0 (, /) (,) 4 Where D is the region in space where the source term q is present and x is position vector.C 0 is the wave speed of sound.Boundary integral methods have been used to investigate the scattering of elastic waves by cracks and cavities (Bouchon, 1987;Rodrıguez-Castellanos et al., 2006;Bouchon & Sánchez-Sesma, 2007) and hydrofractures (Liu et al., 1997).Integral equation methods have been used to model wavefield scattering caused by small scale cracks or inclusions (Muijres et al., 1998;Herman & Mulder, 2002 ).Herman & Mulder (2002) used the integral method with two crack boundary conditions to a homogenous model containing 4000 cracks.Transmitted pressure of two different boundary conditions including compliant crack and rigid crack is illustrated in fig.13.The compliant crack is characterized by zero pressure, and the rigid crack by zero normal particle velocity.The velocity in the model is 3000 m.s -1 .

Direct methods
The last category of seismic forward modeling methods are the direct methods that involve numerical solution of wave equation.Direct methods such as Finite Difference (FD) (Alterman & Karal, 1968;Claerbout, 1985) and finite element (De Basabe & Sen, 2009) require the model to be discretized into a finite number of points and therefore sometimes are called grid methods.The methods also called full wave equation method since it implicitly provides the full wave field.They have the ability to accurately model seismic waves in arbitrary heterogeneous media.FD method is a numerical method for solving differential equations that can be applied to seismic wave equation to calculate displacement at any point in geological models.Seismic wavefield is computed at each grid point by approximating derivatives of the wave equation with finite difference formulas and solving the resulting difference equation Fig. 13.Transmitted pressure for two different choices of the crack boundary condition.The velocity of the embedding equals 3,000 m/s, the crack half width 1 m and the number of cracks 4,000 (Herman & Mulder, 2002) recursively.It includes Taylor series expansions of functions near the grid point.Explicit finite difference methods involve the estimation of wavefield at present time using the wavefield at past times.In implicit FD methods, the present values of the wavefield depend on past and future values.The mathematical formulation of finite difference seismic modeling can be found in several articles (Carcione et al., 2002;Marfurt, 1984;Margrave, 2003;Moczo et al., 2007).Let consider function U(x) with continuous first derivative.Forward, backward, and center-difference equations of the function U(x) are The forward and backward-FD equations are first order approximations to the first derivative and the difference in the value of first derivative and the right hand side in equations 11 and 12 is the truncation error.Finite difference operators can be used to predict a function.For example if we know the function U(x) and its derivative at a point x 0, the function at x 0 +Δx can be derived from equation It is easy to find an approximation to a derivative using Taylor series and the equation can be considered as a truncated Taylor series.An approximation for the second derivative can be derived by applying equations 11 and 12 and a frequently used operator is The equation is a centered approximation of the second derivative that is frequently used in grid-displacement FD methods.
As stated earlier the FD methods require that the geological(computational) domain is characterized by a set of discrete space-time grids.There are three different types of grids depending on the spatial distribution of functions including displacement/particle velocity and stress tensor components in the grids.In the conventional grids all functions are approximated at the same grid positions.In a partly-staggered-grid the displacement components are located at one grid position and the stress tensor components are located at other grid positions.In a staggered grid, each function (displacement and/or particlevelocity component and each shear stress-tensor component) has its own grid position.FD methods have been widely used for seismic modeling of different geological settings (e.g.Alaei & Petersen, 2007;Regone, 2007).
The Seismic Unix implementation of acoustic finite difference modeling (Cohen & Stockwell 2002) was used to generate two shots from different parts of the Zagros model (fig.9).This program uses an explicit second-order differencing method.The source was a 30 Hz Ricker wavelet.Split-spread source/receiver array was used as survey geometry with a maximum offset of 6 km.
The first shot is located on the flank of the anticline located at x=20 km of the model (fig.9) with a thick overburden section.One of the objectives of the modeling is to investigate the effect of thick overburden on the deeper target section.Fig. 14 (left) shows the real seismic data around the first synthetic shot which is a 2D poststack time migrated image.It shows packages of reflectors at some locations and disturbed or overabundance of reflectors in the remaining part of the section, typical in complex fold belts with irregular topography.The shallow picked horizon on the seismic section (fig.14) is located in the thick overburden.
The marked area below the second picked horizon is a complex area that is not imaged as good as the upper part.Fig. 14 (right) shows the common shot gather for the first shot.It is clear that the seismic wavefield in the model is complex.The same picked events in the seismic section (fig.14) have been picked in the shot gather.This single shot that is located on the flank of  (Alaei & Petersen, 2007) the anticline shows how the complexity of the overburden geology affects the seismic response.The types of events that can be seen in the shot record (particularly the diffractions and nonhyperbolic events inside the marked area in fig.14 (right) are significant challenges for pre-or poststack time-processing methods.Advanced seismic velocity analysis and migration methods are required to image these types of complexities.
The second shot is located at X=63 km (fig.9).The main purpose is to investigate the seismic wavefield at the top of a faulted anticline with a complex overburden.The last FD modeling example is from another faulted anticline from Zagros fold and thrust belt.The model (fig.16a) is built using Compound modeling and a single shot has been generated using explicit second-order differencing method (Cohen & Stockwell, 2002).To show the real complexity of the seismic wavefield, a snapshot of the single shot is illustrated in fig.16b.Different wave modes including first arrival reflections and diffractions are illustrated.

Summary
The main goal of this chapter is to introduce seismic forward modeling as a powerful tool to investigate the seismic wave propagation in different geological settings with a special reference to complex geological structures.Seismic forward modeling is the computation of seismic response of a geologic model and has been widely used in both earthquake and exploration seismology.The source of complexities of seismic wave transmission and reflection in subsurface have been explained.I provided some applications of seismic modeling in exploration seismology including its applications in seismic data acquisition, processing, and interpretation of complex structures.Interface and grid based model building approaches presented with examples from complex structures of Zagros fold and thrust belt.Three different seismic forward realizations including asymptotic (ray tracing methods), integral, and direct (e.g.finite difference algorithms) methods have been presented.

Fig. 2 .
Fig. 2. PSTM seismic image from Zagros fold and thrust belt.The incompetent beds exposed at the surface (central part of the line) cause illumination problem for the deeper targets.The picked line (green) illustrates the boundary between competent and incompetent rocks.(Alaei, 2006)

Fig. 3 .
Fig. 3. Seismic image from Nordkapp Basin, Norwegian Barents Sea.The image is complex around and under the salt bodies

Fig. 5 .
Fig. 5. Raypaths from a single shot gather ray tracing with 7km offset from a source located at x=70km.Structural complexity caused poor subsurface coverage between x=70 and x=72km www.intechopen.comShould the model be 2D o 3D?How large should the model be?How many and which surfaces should the model contain?Where should the model properties (interval velocity and density) be obtained?How should the model properties vary between model interfaces?How much complexity (structural or stratigraphic) should be portrayed in the model?

Fig. 6 .
Fig. 6.Interface-based model building approach applied to a faulted anticline.14 interfaces shown in the figure.The shallowest interface represents the topographic surface When the geometry of the model is constructed, seismic forward realizations require properties to be assigned to each of the model layers.These properties include P and S wave velocities and density and can be constant or vary within model layers.The variation can be horizontal or vertical.The representation of properties within each layer reflects geological settings.For example in a siliciclastic sequence properties vary with depth representing compaction trends.The sources of velocities and densities are well data and reflection seismic data.

Fig. 7 .
Fig.7.Geometry of a regional complex 2D model (81x17km) from Zagros fold and thrust belt together with P wave velocity distribution.Well velocity data used to define velocities of model blocks(Alaei, 2005) A strategy for modeling clastic reservoirs was explained byBryant & Flint (1993).It includes five major steps: (1) definition of the space occupied by the modeled interval; (2) recognition of geological units within the model space; (3) assignment of geometries to the units; (4) arrangement of the units within the model space (architecture); (5) assignment of properties to the units.Two common approaches for the third step, assignment of geometry and orientation, are proposed.1. Modeling of discrete objects such as shale in sand or sand in shale.2. Modeling based on continuous variation.This is based on a Boolean method.Both methods use cell-based systems.

Fig. 8 .
Fig. 8. Seismic patterns of a stratigraphic sequence.The reflection terminations at different locations of the unit indicate different geological processes(Vail, 1987) carried out at point x=49.5km of the Zagros model shown in fig. 7. Seismic response of the single shot (fig.10) is generated by convolving travel time data derived from ray tracing with a zero phase synthetic Ricker wavelet of 35 Hz dominant frequency.The nonhyperbolic event shown in fig. 10 corresponds to the repeated layer of the anticline.Due to velocity and structural complexity, rays are bent and there are non-hyperbolic arrival times.

Fig. 10 .
Fig. 10.Synthetic shot gather record at x=49.5km of Zagros model shown in fig. 7. Maximum offset is 4500m.The event indicated by the arrow represent a nonhyperbolic event in the shot gather record(Alaei, 2005)

Fig. 11 .
Fig. 11.Raypath image of Normal incidence ray tracing from part of the Zagros 2D model (fig.7) Fig. 12. 2D seismic image from faulted structure of Zagros fold and thrust belt (left).Synthetic seismic response of image ray tracing with 6km aperture.The arrow illustrates the location of steeply dipping reflector.Increasing aperture from 3 km (real data) to 6 km (synthetic data) improved the signal continuity(Alaei, 2006)

Fig. 14 .
Fig. 14. 2D poststack time migrated seismic image from Zagros fold and thrust belt (left).The marked area in the deeper part is complex and has poor quality image.The right figure illustrates the shot gather.The source position is shown by arrow(Alaei & Petersen, 2007) Fig. 15 (left)  shows the real seismic data around the shot and the arrow shows the location of the shot.The first shallow picked horizon is in the overburden.The other picked horizons are all top reservoir rock units but faulted and repeated at the top of the structure.The overburden deepens and thickens significantly towards the southwest flank of the structure.Fig.15(right) shows the common shot gather.This shot is from the same location as the shot gather record illustrated in fig.10using ray based method.There is much more detail imaged with the FD modeling (fig.15).The shallowest picked horizon is strong in the far offsets and weaker in the near offsets.It represents the high impedance contrast in the overburden.The picked events between 1 s and 2 s are faulted and folded top reservoir rock unit.The arrow shows the faulted leg of top reservoir.It is not easy to explain the complexity of the seismic image in the real data in the overburden section of the southwest flank of the structure using this synthetic shot.

Fig. 15 .
Fig. 15.Real seismic data around the shot at x=63 km of fig. 9 (left).The first shallow picked horizon is in the overburden.The other picked horizons are all faulted (repeated) top reservoir.Common shot gather of the second shot (right).The two picked events between 1s and 2s are faulted and folded top reservoir.The arrow shows the faulted leg that has been confirmed by drilling result(Alaei & Petersen, 2007)

Fig. 16
Fig. 16.a) Geological model of a complex faulted fold from Zagros fold belt.The slowness values are scaled.b) Snapshot from a shot located at the centre of the model.Complexity of the wavefield illustrated by several diffracted events