Seismic Refraction Modeling Using Finite Difference Method and its Implications in the Understanding of the First Arrivals

Y.P. Goyes, S. Khurama, G. Reina, G. Jimenez Perm State University. 15 Bukireva Str., Perm 614990, Russia. E-mail: goyes.yesid@gmail.com Faculty of Geology, Santander Industrial University, Cra 27 Calle 9 Ciudad Universitaria, Bucaramanga, Santander680002 Colombia E-mail: skhurama@uis.edu.co Research Group in Basic and Applied Geology (GIGBA), Santander Industrial University, Cra 27 Calle 9 Ciudad Universitaria, Bucaramanga, Santander 680002, Colombia. E-mail: gjimenezd@gmail.com (Статья поступила в редакцию 2 сентября 2016 г.)


Introduction
As described in Dobrin (1976), the seismic modeling has been an important tool for analyzing various problems mainly in the exploration of hydrocarbons. Their algorithms have been improved in terms of the numerical approximation and propagation in Seismic Refraction Modeling Using Finite Difference Method and its Implications … 257 the elastic mediums (Sayers, 2009). According to Bridle (2016), and due to the inclusion of geophysical methods in the determination of exploration geotechnical data, the seismic modeling was taken to study the behavior of waves on different geometries and acoustic contrast between layers (Budhu, 2010). For example, such as the contact between soil and residual colluvial deposits, where the acoustic impedance is low, and the contact between the colluvial deposits and bedrock, where the impedances high (Redpath, 1973).
In this work, we study the propagation of head waves in an elastic medium over a contact between unconsolidated sediments and bedrock, which enables to see the changes in the first arrivals of the P wave on the surface, showing significant changes in the traveltime curves in the areas with tilt refractors (Lankston, 1988). Traditionally, the modeling of the seismic refraction focuses in calculating the first arrival time (FAT), which is performed using mainly the time-term method (TTM). An example of this method is shown in (Iwasaki, 2002) and (Berry et al, 1966). In this work, to study the first arrivals picking, the method of finite differences (FDM) was used, because it allows one to analyze the times (FAT) and the record of the wave propagation for each geophone. Finally, a comparison using linear fitting, coefficient of determination, and residual mean square between the results obtained by FDM and TTM modeling was made.

Methodology
The modeling by time-term method (TTM) was conducted using the module of the Plotrefra (SEISIMAGER), which includes the capability of creating a custom velocity model for forward modeling purposes. The first arrival time was calculated by time-term method, according to (Bath, 1978). This technique is a linear least-squares approach to determine the best discrete-layer solution to the data. The math behind this technique is comparatively simple. A two layers single receiver-source model is showed in Fig. 1 with the next configuration; V1 and V2 are the velocities for the layer 1 and 2 respectively; ic is the critical angle; Z1 and Z2 are the perpendicular depth measure from source and receiver position respectively. In this case, the time (t) is calculated as shown in the following equation: , S1, S2 -slowness in layer 1 and 2 respectively, c -substitution constrain slowness-dependent, x -distance from source to receiver.  (Iwasaki, 2002) Generalizing the equation (1) for a model with multiple layers, we get in a matrix form: The first arrival time calculation using the finite difference method (FDM) was performed using the free software package SeismicUnix (program SUFDMOD2), which works on the Linux operating system (in this work, Fedora 20). This program uses finitedifference method (Stockwell, 1999) to solve the 2D acoustic wave equation: where v(x,z) is the velocity in the acoustic media.
The Laplacian operator can be approximated with central difference operators. Where the velocity, spatial sampling rate, and grid spacing are in consistent units.
The time derivative is calculated by a second order finite difference scheme: Finally, each finite difference scheme has a stability condition (Slawinski et al, 1999). The stability condition for the second order scheme used by SUFMOD2 is shown in the equation 6. This program does not require the grid spacing to be equal in the horizontal and vertical directions. Figure 2 describes the general flow diagram with the list of processes required in SUFDMOD2 for the generation of synthetic seismic records. To create models with complex geometry, we recommend using the package TRIMODEL, which calculates a triangulated velocity model of subsurface. Subsequently, this model is made uniform by using the program TRI2UNI. In this work, three variants of the same model, which differ only in the grid spacing, were studied. In the Table 1 z and x are the spatial interval in the z-axis and x-axis respectively, and nx and nz are the sample numbers in the z-axis and x-axis respectively.
Vmax -maximum velocity in the model (m/s); t -time interval in seconds; x -spatial interval in the x-axis (meters). We used for both modeling methods the general model shown in Fig. 3 with the area of thickness of 2 m simulating a geometry of the completely flat refractor. The zone B (1250 m/s) has a variable thickness 0 -4 m that corresponds to a refractor with complex geometry simulating the advanced processes of weathering on the bedrock indicated in the general model by the area C. We used a sequence of five shot-points with different location. Table 2 shows the position of each shot.

Fig. 4. Uniformly sampled velocity model of layered medium (M3) in binary format used for SUFDMOD2 program
The uniform models (Fig. 4) were built using linear interpolation of the data obtained by seismic tomography. SeismicUnix creates different models using the program UNI2, which generates a 2-D uniformly sampled velocity profile from a layered model. In each layer, velocity is a linear function of position.

Results
The modeling with finite differences for the velocity model with 3 different grid spacing (Table 1) produced a variation in the optimization of stable t (6) and the number of iterations (Table 3).
A seismogram with the traces corresponding to the registration of the P wave to geophone spread was obtained for each model settings, which differ in values nx (number of samples in x-axis) and nz (number of depth-samples).
It was found that the variation in the number of geophones used for seismic acqui- sition has a significant impact on the quality of registration of seismic events. In the seismogram with nx=24 (Fig. 5), the first break of refracted wave is reliably identified, but the accurate arrival time of the direct wave close to position of the source is difficult to define. Using nx=48 (Fig.6), we are able to improve the sampling of the waves between 0.04 and 0.05 seconds, and only seismogram obtained with nx=500 ( Fig. 7) shows clearly the head waves 1 and 2, and direct wave.
The results obtained using TTM and FDM methods are illustrated in Fig. 8. The time of arrival of the P wave on the same model using nx = 24 was calculated with each method. In Fig. 8A, the comparison between the travel-time curves shows a high dispersion of the data for time less than 0.015 s. For times exceeding 0.015 s (first refraction time in Fig. 8C), the evidenced changes in the slopes of the curve corresponding to the different layers of the subsoil are clearly recognized. The highest difference (FD -TTM) about 0.008 s (Fig. 8B) corresponds to the direct wave. To analyze the dispersion of the results, we performed a linear fit between them shown in Fig. 8C.
This curve was divided into two zones (high correlation zone and low correlation zone). The value of calculated coefficient of determination was determined as 0.956866 and the straight fit equation was defined as timeRT=1.07871478*timeFD -0.00284832, where timeRT and timeFD are the times calculated from TTM and FDM data respectively.

Analysis of the seismic survey data
A seismic refraction survey in Barrancabermeja (Colombia) was developed with the 260 Y.P Goyes, S.Khurama, G. Reina, G. Jimenez aim to identify changes in the subsurface and define the depths and types of soils.
The seismic acquisition was performed using a seismograph ES-3000 with 24 vertical 14 Hz geophones separated 2 meters. The total length was 46 meters. The result of the seismic processing is shown in Fig. 9. The results of seismic tomography processing correlate with the data of the geotechnical survey. For optimization the temporal parameters in future acquisition surveys over the same area, we performed a seismic modeling using FDM. Additionally, the modeling allowed us to know the validity of the obtained results by comparing the synthesized traces with the survey traces (Fig. 10).
The events associated with the refraction waves were found between 0 and 0.05 seconds. The seismic survey should decrease the maximum time of the acquisition but increasing, on other hand, the detail of the head waves to improve the process of first arrivals picking.
The arrival times obtained in the seismic record x=0m using FDM and real data, procurement of field (RD) are shown in Fig. 11.
The slope preserves the general trend, but the data differ in terms of lap times. In blue curve shows the difference between two times. The maximum difference is 0.011 seconds.

Conclusions
Modeling using the finite difference method showed precision in the calculation of the first arrival time of P-wave, which was verified using statistical correlations among the times calculated by the TTM method and the times obtained from the seismogram (first arrival time picking using the FDM method). The modeling allowed an analysis and study of seismic survey on models with different geometrical configurations and variation of the seismic velocity with depth. With respect to the number of geophones necessary to record the seismic events involved in seismic refractions, 24 and/or 48 geophones