Application of Geophysical Methods for Geotechnical Parameters Determination at New Borg El-Arab Industrial City , Egypt

Due to a rapid increase in the population during the last few decades, the banks of the Nile River and its delta have reached maximum capacity. As a consequence of this increase, the Egyptian Government has constructed a number of new urban areas and industrial cities outside the Nile Delta. New Borg El-Arab City is one of these new industrial cities. This city is located around 60 km southwest of Alexandria City. This industrial city is proposed to include an airport, a number of factories, worker settlements and heavy truck roads. Therefore, a detailed study of site characterization should be performed before construction being in order to. The main target of this study is to determine the dynamic characteristics and geotechnical parameters at the proposed site using seismic refraction and electrical resistively techniques. Analysis and interpretation of the obtained results reveal that the subsurface consists of three layers with a gentle general slope toward the Mediterranean Sea. The classification of rock material for engineering purposes reveals that the study area is divided into three zones.


Introduction
New Borg El-Arab industrial zone is one of the new industrial cities in Egypt.It lies between latitudes 30.75043221 and 31.04684093N and longitudes 29.44052705 and 29.6723497 E, covering an area of about 90 km 2 (Figure 1).This new industrial city is an urbanization project that is relatively large-scale regional development project and it will be good if the determination the dynamic characteristics and geotechnical parameters being incorporates in the planning stage.So, the main target of the present study is investigating the shallow subsurface structure conditions, the dynamic characteristics and geotechnical parameters of subsurface rocks using shallow seismic refraction profiling and electrical resistivity tomography surveys.The output of this study is very important for solving problems, which associated with the construction of various civil engineering purposes, land use-planning and earthquakes resistant structure design.The integrated interpretation of these techniques classifies the subsurface succession into three layers.

Geological Setting
According to (El Shaazly, 1964), the geological setting of the study area is similar environment consisting of Quaternary deposits (scattered), Oligocene (Sand), Beniolitic Clay, and Middle Miocene as Limestone (from top to bottom), Figure 2. The geomorphological units west of Alexandria showed by (Abd El Mawla, 2010) as in Figure 3.

Methodology
In the present study two from the advanced technuqes have been applied; electrical resistivity tomography and shallow seismic refraction.Electrical resistivity tomography is a geophysical technique for imaging and mapping the vertical and horizontal sub-surface structures from electrical resistivity measurements made at the surface.The shallow seismic refraction technique is considered as one of the most effective geophysical method, which can be applied in the field of engineering seismology i.e. tunnels, dam sites, landslides, quarries, roads, reclemaited lands, caves and cavities.

Electrical Resistivity Tomography
In the ERT method the distribution of the electrical resistivity of the subsoil is obtained by injecting electrical current (by the current electrodes) into the ground and measuring the potential difference (by the potential electrodes) at two determined points of the surface.The method is based on the application of Ohm's law:   where ρ a is the apparent resistivity (Parasnis, 1997), k is a geometric constant that depends only on the reciprocal positions of the current and potential electrodes; ΔV is the measured potential difference, and I is the intensity of the injected current.The apparent resistivity values depend on the true resistivity distribution.
The true resistivity distribution in the investigated medium can be estimated by an inversion procedure based on the minimization of a suitable function.This function is generally the sum of the squared difference between measured and calculated apparent resistivities.The investigated medium is discredited in a 2D (or 3D) grid of cells, where each cell is assigned an initial resistivity value.A finite-difference (Dey & Morrison, 1979a, b) or finite-element (Silvester & Ferrari, 1990) procedure computes the predicted apparent resistivity at the surface.The solution to the problem, as is well known, is not unique.For the same measured data set, there is wide range of models that can give rise to the same calculated apparent resistivity values.To narrow down the range of possible models, normally some assumptions are made concerning the nature of the subsurface (i.e.geology of the subsurface, whether the subsurface bodies are expected to have gradational or sharp boundaries) that can be incorporated into the inversion subroutine.The current method of solution minimizes the difference between measured and calculated apparent resistivities using the smoothness-constrained inversion formulation, which constrains the change in the model resistivity values to become smooth (Loke & Barker, 1996;Loke, 2001).The "smoothness-constrained robust inversion" method (Loke, 2001) has proved to be much more useful when the subsurface bodies have sharp boundaries (Loke, 2001).

ERT Data Acquisition
The ERT measurements are carried out along 56 profiles oriented approximately N-S (Figure 1).For this type of data set, a series of 2D inversions is usually first carried out.The 3D inversion is then used on a combined data set with the data from all the survey lines.Apparent resistivity data were collected using the Syscal R2D instruments in multi-electrode configuration.The distance between profiles was 2 km and a total of 24 electrodes spaced at 5-m intervals are employed.Data are collected using the pole-pole array.The pole-pole array was chosen because, as is well known, it is very sensitive to horizontal changes in resistivity, and is therefore suitable for mapping vertical and horizontal structures (Loke, 2001).Figure 4 shows the principle of ERT technique.

ERT Data Processing and Interpretation
The least-square inversion by Jacobian matrix calculation was used in order to obtain an image of the true resistivity distribution as a function of depth.The least-square inversion method (Loke, 2001) minimizes the absolute changes (l1-norm) in the model resistivity values (Loke, 2001).All pole-pole inverted resistivity sections produced low RMS errors (<13%).The Res2dinv software (Loke, 2001) is used to invert all 2D and 3D data sets.Figure 5 shows an example for inverted resistivity section with tomography data along profiles "1, 17, and 56".
As the measured sections are straightforward to combine all 2D images and produce resistivity slices for different depths.The Res3dinv resistivity inversion software (Loke, 2001) is used to automatically invert the apparent resistivity acquired data and to yield a three-dimensional resistivity model.The slicer-cubic software (Samuil, 2008) is used to automatically draw output interpreted resistivity data from Res3dinv software to acquiesce a three-dimensional resistivity model (Figure 6).From this figure it can be concluded that, the topmost layer presents a larger area with high resistivity variations over short distances (about 2100 Ωm to 1500 Ωm).In comparison, second layer, lying at a depth ranging from 13.35 to 16.75 m, shows more gradual lateral variations in the model resistivity values, this layer presents moderately resistivity values (ranging from 1500 Ωm to 600 Ωm).Third layer lying at a depth varies between 22.25 and 24.60 m.The resistivity values of it are clearly visible, ranging from about 600 Ωm to about 200 Ωm.
Seismic exploration involves generation of seismic waves and recording the arrival times of these waves from the source to the series of geophones (Figure 7).Seismic refraction is used to evaluate the necessary parameters    for constructions, or to solve the problems related to the geological nature of sub-surfaces, mining works, and the environmental conditions overcame in the site (Stumpel et al., 1984).
In the shallow seismic refraction method, the seismic waves, created by artificial sources such as an auto-hammer viberator, propagate through the medium and are refracted at interfaces, where the seismic velocity or density changes.Geophones laid on a single line record the waves returning to the surface after travelling Hole different distances through the ground.By measuring the travel time between the break and the recording of a seismic signal, the seismic velocity in the subsurface and the depth of the interfaces may be inferred.Figure 7 schematic diagram shows the principle of the shallow seismic refraction.
Conventional analysis of seismic refraction data sets makes simplified assumptions about the velocity structure that conflict with observed heterogeneity, lateral discontinuities, and gradients (Parasnis, 1997).Refraction tomography is designed to resolve velocity gradients and lateral velocity changes, enabling it to be applied in settings where traditional techniques fail.
The method used in this paper utilizes non-linear travel time tomography consisting of ray tracing for forward modeling and simultaneous iterative reconstruction technique (SIRT) for inversion.In this method, the velocity model is represented by quadrangle cells (Figure 8).The width of each cell is chosen as the receiver interval.First-arrival travel times and ray paths are calculated by the ray tracing method based on Huygen's principle (Parasnis, 1997).A ray is expressed as a line connecting the nodes arranged on the cell and the travel time between a source and a receiver is defined as the fasted travel time of all ray paths (Figure 8).A model is updated by the SIRT (for more information see (Morey & Schuster, 1999), (Nemeth et al., 1997).

SSRT Data Acquistion
The target to be executed was the determination of dynamic properties of soil and foundation rocks by using seismic refraction method and this implied firstly to throw light upon the end product of seismic refraction technique that was used in this study.Shallow seismic refraction technique was applied in the present study to determine the subsurface layering and its dynamic characteristics at the study area.The seismic refraction data have been acquired along 56 shallow seismic refraction profiles spread over the same area proposed for the R2D resistivity imaging survey (Figure 1).The survey has been accomplished using 24 channels signal enhancement seismograph "GEOMETRICS SMARTSEIS" along 120 m length profiles.Three shots are selected on the seismic line for measuring perpendicular and surface waves: two near shots at the edges of the line, and one shot at the middle of the line.In the surface P and SH wave velocities have been measured.The hammer connects to sensor to provide the time break to the seismograph.The power of the hammer helps to avoid loosing of waves strength that may be caused by Blind layer in some condition.The low pass filter in the recording system has 7 -10 MHz as frequency response, which is suitable for the recording condition, positioned before the analog-digital conversion circuit.Data were stacked at least three times for each source.Figure 9 shows the measured raw data (typical recorded seismic traces) of profiles number 1, 17, and 56 as examples.

SSRT Data Processing and Interpretation
For interpretation of refraction data an initial depth model is obtained using Gardner's method (Gardner, 1939).To perform the processing and interpretation of the seismic refraction tomography data in the current study, (SIPEEDIT (V. 3) 2002) software developed by OYOO Company is used.Therefore, an initial velocity model is estimated.In this case, the initial velocity model is represented by the results obtained from the simple interpretation of refraction data (such as the model shown in Figure 9).The model is represented by 1-1 m cells.Refraction rays are traced through this model to give calculated travel times.A misfit function, consisting of the squared difference between the observed and computed travel times, is calculated.The model is adjusted until the misfit is minimized.The iterations are stopped when the RMS travel time residual (difference between  the calculated travel times for the initial model and the observed ones) is less than the average travel time pick error.

SSRT Velocity
The model generated by seismic refraction tomography shows the variation of the Vp and Vs in the subsoil.The seismic refraction tomography survey points out that the shallow subsoil may be divided into three main layers.The surface layer exhibits a P-wave velocity of about 734 to 772.2 m/s, and a Vs velocity of about 389 to 408 m/s, these velocities correspond to surface sandy weathered zone with thickness ranges from 11 to 13.4 m.This layer corresponds to the area in which the ERT profile shows high resistivity values (Figures 5 and 6).These velocity values are typical of more weathered subsoil (Leucci, 2004).The second layer shows a P-wave velocity ranges from 1000 to 1120 m/s, and a Vs velocity ranges from 521 to 606 m/s, these velocities correspond to calcareous sands weakly cemented with thickness varies between 15 to 16 m.This layer is characterized by the mainly high seismic velocities (Vp and Vs), corresponds to the area about where the ERT profile shows the moderately low resistivity values.The third layer shows a P-wave velocity between 1400 and 1550 m/s, and a Vs velocity between 669 to 781 m/s, those correspond to sandy marl layer strongly cemented with depth vary from 25 and 26 m where the ERT profile shows the moderate resistivity values.The third layer is characterized by the highest seismic velocities (Vp and Vs), corresponds to the area about where the ERT profile shows the lowest resistivity values.The grey zone presents the no-ray area.
The increase of seismic wave velocity could be due to a more compact soil (calcareous sands weakly cemented).According to the variation in Vp and Vs velocities are detected, four cross-sections, over the study area, have been drawn to show the distribution of these layers in the four sides (Figures 10-13).Three dimension maps of the distribution of Vp and Vs velocities for the three layers over the study area are shown in Figures 14  and 15.

Density Gradient (Di)
There is a direct relation between the material competence and the value of the density gradient and the compression velocity (Vp)."The compaction" can be defined as "name of the physical mechanism that converts sediments from their initial state to a progressively denser state under the influence of their own weight or of tectonic movements" (Cordier, 1985).The density gradient (Di) relates to how much consolidation or settlement will take place.Many relations can express it "according to Stumpel et al. (Stumpel et al., 1984)" ( ) The distribution of the density gradient in the subsoil layers has been shown in Figure 16.Generally over the three layers, the relatively low values are observed in the southeastern part of the study area, which indicates the relatively moderate competent materials.The middle values are occupied in the middle part that may be due to difference in sorted materials.While the high values are observed in the northwestern part of the study area, they reflect the fairly moderate competent materials in this part of the study area.Table 1 shows the distribution of the density gradient.

Standard Penetration Test (SPT) [N-Value]
The Standard Penetration Test (SPT) is geotechnically known as "the resistance to penetration by normalized cylindrical bars under standard load".N-value is geophysically evaluated using Imai's (Imai, 1976;Stumpel et al., 1984), which is given by the formula: Vs 89.9N0.341= Table 2 shows the distribution of N-Value.In general for the three layers of the study area, Figure 17 shows the comparatively low values that are experimented in the southeastern part of the study area, which point toward the relatively moderate competent materials.The middle values are inhabited in the middle part that may be suitable to difference in sorted materials.Whereas the high values are observed in the northwestern part of the study area, they reflect the fairly moderate competent materials in this part of the study area.

Kinetic Bulk Modulus (K)
"Sheriff" defined the Kinetic Bulk modulus as "the stress-strain ratio under simple hydrostatic pressure" (Sheriff, 1991).The law can calculate the different values of Bulk modulus to each layer: where ρ is the density.
Table 3 shows the division of Bulk modulus.Over the three subsoil layers, the minimums values of Bulk modulus are observed in the southeastern part of the study area.These values increase to be moderate ones in the center part, while the maximum values are observed in the western and the northwestern corners of the study area (Figure 18).

Ultimate Bearing Capacity (Qult)
The ultimate bearing capacity (Qult) can be defined as "the maximum load required for shear failure or sand li-  quefaction".This capacity is controlled by shear strength factor.The ultimate bearing capacity for cohessionless soils can be calculated by using the standard penetration test (SPT) from Parry's formula (Parry, 1977) as: where "N value" is the resistance to penetration by normalized cylindrical bars under standard load.
Figure 19 shows the overview of the three layers in the study area.The minimums values of ultimate bearing capacity are observed in the southeastern part of the study area.These values increase to be moderate ones in the center part, while the maximum values are observed in the north and the northwestern corners of the study area.Table 4 reveals the distribution of the ultimate bearing capacity (Qult) in the subsoil layers of the study area.

Allowable Bearing Capacity (Qa)
The allowable bearing capacity is the maximum load to be considerable to avoid shear failure or sand liquefaction.It can be termed as allowable bearing pressure too, the foundation materials are affected by both strength and deformation characteristics.
The allowable bearing capacity can be calculated by dividing the ultimate bearing capacity value (Qult) by suitable factor of safety (Abd Elrahman, 1989) as: Qa Qult F = (6) The safety factor (F) equals 2 when the soil is cohessionless and equals 3 when the soil is a cohesive material.Table 5 reveals the distribution of the allowable bearing capacity (Qa) in the subsoil layers of the study area.The minimum values of allowable bearing capacity are observed in the southeastern part of the study area.These values increase to be moderate ones in the center part, while the maximum values are observed in the north and the northwestern corners of the study area (Figure 20).

Rock Material Calssification
This study suggests the classification of foundation rock materials in the study area for engineering purposes based on all the calculated moduli and parameters (Figure 21).This study is divided (the selected study area only) into:

1-Zone (A): (competent material zone):
This zone covers the northwestern corner and the small strip in the middle quarter of the study area.Good competent materials characterize this zone suggesting suitability for any engineering purposes.

2-Zone (B): (Moderated competent material zone):
This zone covers the middle quarter with small strip in the southwestern part of the study area.This zone suggests suitability for moderately and low-load engineering purposes.

3-Zone (C): (non-competent material zone):
This zone covers the eastern part of the study area with some extent in the north.This zone consists of less competent materials not suitable for limited construction purposes.This zone is suggesting non-suitability for any engineering purposes that must be kept away from any constructions or loadness and watery activities.

Discussions and Conclusions
In order to map the Bourg Al-Arab industrial city (Alexandria, Egypt), electrical resistivity tomography (ERT)    Electrical survey, effectively visualized as iso-resistivity surfaces, classified the area into three layers.Shallow seismic refraction tomography correlates well with ERT.The integrated interpretation of seismic refraction and resistivity tomography makes it possible to reduce ambiguity.The proposed method of 3D visualization (iso-resistivity and iso-velocity) confirms the interpretation of the 2D standard horizontal sections.
Measuring Vp and Vs velocities is able to calculate the geotechnical parameters of the study area; these parameters classify the area, according to the compaction, hydrostatic pressure, the resistance to load penetration, and load of shear failure or sand liquefaction, into three zones.The geotechnical parameters categorize the area

Figure 1 .
Figure 1.Location map of the study area and the site of ERT and SSRT stations.

Figure 4 .
Figure 4. (a) Sequence of measurements to build up a pseudo section using a computer controlled multi-electrode survey setup; (b) Arrangement of blocks used in a model together with data points in pseudo section.

Figure 6 .
Figure 6.Three-dimensional resistivity model of the study area.

Figure 7 .
Figure 7. Velocity estimation using the slope method calculation with shot point and geophones array in field.

Figure 8 .
Figure 8. Principle of the ray tracing.

Figure 14 .
Figure 14.Three-dimension velocity map for seismic perpendicular waves.

Figure 16 .
Figure 16.Three-dimension map of density over the study area.

Figure 17 .
Figure 17.Three-dimension map of Standard penetration Test [N-Value] over the study area.

Figure 18 .
Figure 18.The allotment of Kinetic Bulk Modulus (K) in the study area.

Figure 19 .
Figure 19.Three-dimension map of the ultimate bearing capacity over the study area.

Table 4 .
The ultimate bearing capacity (Qult) over the study area.
-refraction tomography (SSRT) methods are used.The results highlight the reliability of the integrated data interpretation based on two physical parameters with different resolution and sensibility.

Figure 20 .
Figure 20.Three-dimension map of the allowable bearing capacity over the study area.

Figure 21 .
Figure 21.Classification of the foundation rock material quality for engineering purposes according to the geotechnical characteristics in the study area. 29

Table 1 .
The distribution of the density gradient of the study area.

Table 2 .
The distribution of (N-Value) over the study area.

Table 3 .
The distribution of Kinetic Bulk modulus (K) over the study area.

Table 5 .
The allowable bearing capacity (Qa) over the study area.