Application of seismic refraction tomography for tunnel design in Santa Clara Mountain, San Juan, Argentina

Armando Luis Imhof¹, Manuel Sánchez², Carlos Calvo³ and Adriana Martín⁴

¹ Instituto Geofísico Sismológico Volponi / Facultad de Ciencias Exactas Físicas y Naturales. UNSJ. Ignacio de la Roza y Meglioli. Rivadavia. C.P. 5400, San Juan, República Argentina. Fax: (54) 264 4234980. Tel.: 54 264 4945015. E-mail: aimhof@unsj.edu.ar

² Instituto de Mecánica Aplicada. Facultad de Ingeniería. UNSJ. E-mail: msanchez@unsj.edu.ar

³ Departamento de Matemática Aplicada. Facultad de Ingeniería. UNSJ. E-mail: ccalvo@unsj.edu.ar

⁴ Departamento de Informática. Facultad de Cs. Exactas, Fcas y Naturales. UNSJ

Record

Manuscript received: 15/12/2010 Accepted for publications: 30/11/2011

ABSTRACT

A geophysical survey involving seismic refraction tomography (SRT) for mapping 'P' waves was carried out in Sierra Santa Clara, San Juan Province, Argentina in July 2009. The purpose of the geophysical survey was to determine the degree of fracturing and the rigidity of the rock mass through which it is planned to build a 290 m long road tunnel traversing the mountain almost perpendicular to the axis thereof, at around 100 m depth from the summit.

Several difficulties arose from the operational point of view which made it almost impossible to conduct fieldwork in normal circumstances. Firstly, the topography had almost 45° slopes and 100 m research depths which would have involved having had to use explosives to generate seismic waves reaching sensors which had sufficient signal-to-noise ratio for distinguishing them. Legal restrictions regarding the use of explosives on the one hand and insufficient power when using hammer blows on the other made it necessary to design and build a gas-powered gun to achieve the minimum energy (2 kJ) required for detecting seismic signals.

Secondly, using conventional interpretation methods involving layered models was inoperable in such geological structures; seismic tomography methods were thus used which make use of the velocity gradient concept (both lateral and in-depth). This allowed mapping subsurface velocity variations in the form of velocity contour lines.

The methodology used with the new seismic waves' source generator, as well as SRT application in this type of geological structure, demonstrated that satisfactory results could be obtained for this kind of geophysical study for geotechnical purposes.

Keywords: seismic tomography, P wave, seismic refraction, compressional waves, signal processing.

RESUMEN

Se llevó a cabo un estudio geofísico de tomografía sísmica por refracción (SRT) con ondas compresionales en la Sierra de Santa Clara, provincia de San Juan, República Argentina en julio de 2009. El propósito del estudio geofísico fue determinar las variaciones de velocidad de ondas P en el interior del macizo rocoso con el fin de inferir el grado de fracturación y rigidez del mismo, y proyectar la construcción de un túnel vial de 290 m de longitud, que atravesará la montaña de forma casi perpendicular al eje de la misma, a 100 m de profundidad desde la cumbre.

Varias dificultades se presentaron desde el punto de vista operativo, que tornaron altamente compleja la realización del trabajo de campo en circunstancias normales. En primer lugar las condiciones topográficas, con pendientes de 45° y profundidades de investigación de 100 m, que implicaba el uso de explosivos para generar ondas sísmicas que alcanzaran los sensores con suficiente relación señal-ruido para distinguirlas. Restricciones legales respecto al uso de explosivos, por una parte, y la potencia insuficiente en caso de uso de golpe de maza por otra, hicieron necesario el diseño y construcción de un cañón impulsado por gas en busca de generar la energía mínima (2 kJ) necesaria para la detección de las señales sísmicas.

En segundo lugar, no fue posible la utilización de los métodos convencionales de interpretación que emplean modelos de capas, inoperable en este tipo de estructuras geológicas; por ello resultó imprescindible utilizar el método de tomografìa sísmica que hace uso del concepto de gradientes de velocidades, tanto laterales como en profundidad. Ello posibilitó mapear las variaciones de velocidad en subsuelo en forma de curvas de nivel.

Palabras claves: tomografia sísmica — sísmica refractiva — ondas de compresión — procesamiento de señales.

La metodología empleada con el nuevo sistema de generación de ondas sísmicas, así como la aplicación en este tipo de estructuras geológicas de la SRT, demostraron la posibilidad de obtener resultados satisfactorios para este tipo de estudios geofísicos con fines geotécnicos.

Introduction

Conventional refraction inversion methods use a layered model approach. The subsurface is divided into a number of continuous constant velocity layers having velocities and thicknesses which are modified through interactive forward modeling in an effort to match traveltimes determined from the field data.

Several interpretation methods have been based on this premise (see among others Jakosky, 1940; Heiland, 1940; Palmer, 1986; Dobrin, 1975).

These conventional methods require sections of traveltime curves to be mapped to refractors, a task that can be difficult in several geological situations. For example, the presence of extremely hard topographic conditions in mountains surely implies that there can be large and sudden changes in the shape of the bedrock. There may also be localised feature? such as voids that contradict the assumption of continuous constant velocity layers.

The progress of computer technology has simplified calculus procedures, thereby leading to a great development in geophysical data processing methodologies, particularly for application to soil studies regarding civil engineering infrastructure. The amazing increase in microcomputer processing velocity during the last few years, added to their huge storage capability, has enabled the development of algorithms from different theoretical backgrounds, giving birth to refraction tomography (Schuster et al 1993; Sethian, 2007; Zhang et al., 1998).

Refraction tomography

Unlike conventional refraction methods, seismic refraction tomography (SRT) does not require that the model be broken into continuous layers having constant velocity. Instead, the model is made up of a large number of small constant velocity grid cells or nodes. Inversion is performed by an automated procedure involving raytracing through an initial mode and comparing the modelled traveltimes to the field data, and adjusting the model grid-by grid to match the modeled traveltimes to these. Thif process is iteratively repeated until a preset number of iterations has been reached. Because there is no assumption of continuous constant velocity layers, SRT can model localised velocity anomalies.

Another difficulty related to traditional seismic refraction methods for data interpretation is that they are only applicable to relatively gentle topography, as mentioned in previous paragraphs.

Several difficulties appear regarding geophysical studies in mountain areas, namely:

• The terrain is extremely rugged; and

• Such material has great lateral variation, making layered model theory inapplicable.

Survey area location and general geology

The study area was located in Calingasta, Sarmiento, in the southern sector of the province of San Juan, Argentina, and the district of Las Heras in the north of the province of Mendoza and about 150 km southwards of the capital city of San Juan, on the eastern slope of the front range and western edge of the pre-cordillera at an average height of 2,500 m above sea level (masl) (Figure 1).

The region's climate can be described as semi-arid to arid, being of the mesothermal type.

The present study was conducted in a mountainous area having an ancient marine environment, with the presence of high-angle folding sedimentary rocks, general lithology consisting of sandstones interbedded with siltstones.

The area is characterised by the presence of two major mountain bodies, the west one belonging to the Front Range and the eastern one belonging to the pre-cordillera with a central mountain body, corresponding to the central and western pre-cordillera subprovinces (Figure 2)

The Pedernal-Los Pozos hills form an outstanding feature of the eastern mountain range, extending 25 km NNE-SSW and 8 km W-E. They are mainly made up of Paleozoic age calcareous rocks, being more than 2,000 masl high.

Figure 2 shows the prevalent geological formations in the study area. The main formations are as follows:

• The Marquesado group (Bordonaro, 1980), which includes La Laja, Zonda and La Flecha formations and is composed of limestone, dolomitic limestone, marl, shale and chert;

• The Mariño formation, made up of sandstones and tuffs presenting strongly structured (high tilt) intercalated materials; and

• The Uspallata Triassic formation, which is composed of medium to coarse-grained olive green metasandstone alternating with olive green shale (Cortés et al., 1997).

There is abundant folding and fracturing (jointing). The formations dip at a high angle towards the NW, where there are greater shale levels.

Basic theory

Refraction interpretation methods leads to seismic velocity determination (typically calculated from the first break pick of the P-wave) as part of the interpretation results.

Travel time inversion is nonlinear because a perturbation in the velocity field results in changes in the ray paths through it. Linearised inversions involve local linearisation of such nonlinear problem. This makes them initial-model dependent. A poor choice of starting model could cause an inversion to converge to an incorrect solution at a local minimum, due to most algorithms looking for a solution close to the starting model. Obviously, this hampers interpretation.

Global Optimisation

Several investigators have developed algorithms that search for global minimum, instead of local ones. The algorithm used in this work for seismic tomography data interpretation used this type of scheme; it is known as simulated annealing optimisation and it was developed by Kirkpatrick et al., (1983) and Bohachevsky et al., (1986). Although the original theory was developed for statistical mechanics related to temperature, it soon became adapted to geophysical problems (e.g. Chunduru et al., 1996). It follows the basic procedure given below:

1. Travel times were computed through an initial model. There are several ways of doing this, basically grouped into two methods: straight or curved ray theory (Cerveny, 2001) or numerically resolving the wave equation with some constraints, such us the eikonal equation (Vidale, 1988). The latter option has two main advantages: it accounts for curved rays from all types of primary arrivals (i.e. direct ones, refractions or diffractions) and a ray tracing algorithm does not need to be used.

There are several schemes for the numerical solution of the wave equation, most being grouped into two methods: finite differences and finite elements methods. Here, a fast finite difference algorithm was used. The domain was divided into rectangular blocks or cells having a lateral extension generally half the distance between geophones and depth was defined by the interpreter.

2. The least-square error, E_g, was calculated among the measured times (tmeas) and the first theoretical predicted ones (tpred). The least-square error was defined for any iteration it, which is named the objective function, as:

where n is the number of observations and j denotes each observation.

2. The velocity model was then altered by repeating step 1 and the new least-squares error El was computed.

3. The new model was unconditionally accepted if its error was small (El <= Eg ), and conditionally accepted with probability:

when E1>E0 . In the above equation, AE = E0 - E1, Twas temperature, q was an empirical parameter and Emin global minimum objective function value (ideally, equal to zero). Conditional acceptance of models having larger least square errors enabled the algorithm to escape from local minima in its search for the global minimum. (Emin-El)q made the probability of accepting a detrimental step towards local minima tend to zero as the inversion approaches the global minimum.

4. The procedure was repeated until the optimisation converged. The software performed thousands of iterations to obtain the best matches between:

i) Actual arrival times and distances between sources and receivers; and

ii) Modelled, calculated times based on ray paths through the finite-difference mesh.

Each finite-difference array element forming part of ray paths was characterised with a velocity value throughout the optimised interpretation. The extension of the area being researched was the part of the finite-difference mesh having optimised velocities. The interpretation results were typically presented in colour bar form, different colours representing different velocities. Interpretation results could also be converted into spreadsheet files for further analysis.

Methodology

Figure 3 a-c shows a top view and cross-section, as well as an overview of the study area where source points T1 to T3 were located at the coordinates shown in Table 1.

Figure 3-b shows that the trace of the projected tunnel is about 100 m deep from the apex of the mountain, which made logistical tasks become very difficult.

A 24-digital-channel seismic device was configured from 12 physical-channel recording equipment having a three pointed shot: a central one (T2) and two at both ends of the profile (T1 and T3). Figure 3.b represents the seismic survey location along the hillside.

Due to the seismograph having 12 channels, shooting had to be repeated at each emission point (in addition to particular shotpoint's corresponding stacks), moving the line of geophones each time to cover all 24 channels. The total layout length was 280 m with 5 m to 15 m spacing between stations (Figure 4a). Figure 4b shows the distribution of in-depth information content (coverage) regarding the rays with this three-shot device according to the topography in question.

Source wave system

A new system for generating body waves (P and S), called Tamper (Imhof, 2009), was designed and implemented with the aim of generating waves at the corresponding shooting points (T1-T3, Figure 3 a-b and 4). Normally, a 5 kg to 10 kg hammer is used as P wave source, making it possible to generate power up to 0.2 kJ (Fernandez, 2000). This power scarcely allows reaching 20 m research depths, which was insufficient for this study. Another option was to use explosives, which was then discarded as the work also required S wave generation for other kinds of measurements (not described here, see Imhof, op.cit.).

Ordinary seismic wave generation systems use various energy sources, such as elastomeric bands (Damper), gunpowder or compressed air buffalo-gun type (Pullan et al., 1987), to name but a few. The excellent articles written by Miller et al., (1986) and Herbst et al., (1998) have compared several of them (see also Figure 5). These systems are not very portable, and basically they only allow wave trains to be generated (since they emit vertically). The Tamper works like a cannon loaded with nitrogen gas at high pressure (150 bars) as a source of energy that fires a 3 kg bullet at high velocity. The release of energy is variable and depends on the bullet's trajectory length and gas volume in the chamber, being able to reach up to 20 kJ. It can be placed horizontally to generate SH waves, as well. Figure 5.b shows an image of this device.

2.5 kJ were enough for this work, which allowed reaching a 100 m depth research (Herbst et al., 1998) (Figure 5.a) with a seismic device like that shown in Figure 6.

Results

Time-distance graphs (dromochrones)

Six 500 msec. recording time records were obtained for P-waves (see Figure 7.a; a typical record for T1). Contiguous records for the same shot were assembled in pairs to build one of the 24 channels which allowed observing consistency among the first arrivals, thus facilitating choosing them.

The tasks of extracting information from records and transferring it to the time distance graphs were very complicated due to the poor quality of the seismic signals, resulting from prevailing weather conditions (wind) in the study area combined with the terrain's extreme unevenness which hampered geophone deployment.

Pre-processing

Records in SEG-1 format were transformed into ASCII with SEG1-TO-ASCII software prepared by the working group; each signal was then processed separately with bandpass filters. The whole record was reconstructed to enable choosing them more accurately and achieve a more adjusted dromochrone. Filtering algorithms (Fourier transforms) were prepared in Mathcad and Matlab platforms.

After pre-processing and choosing all records, data were represented in the form of time-distance graphs (dromochrones; Figure 7.b).

Input information required by the seismic modelling software consisted of the device's geometry (Figure 4.a) and the first arrivals (i.e. dromochrones).

The grid size for finite difference modelling depended on the tomogram's desired resolution; the highest was used here, so minimum grid size was considered, applying the expressions (Pullammanappallil et al, 1994):

Figure 8 shows contour lines that represented boundaries between compressional wave velocity gradients.

Discussion

The area being investigated consists of coarse-grained sandstones intermixed with dark shale (the latter most noticeable on the NW side), having Bz = 65° NW and Az = 100°. The rock is extremely fractured (altered).

Regarding shale, it should be noted that it has a wide range of seismic wave velocity variation, according to the degree of fracturing. For example, velocity can reach up to 1,500 m/s when weathered, increasing significantly up to 5,000 m/s if not having undergone external alteration.

The seismically surveyed area could be divided into two sections: the first one from T1 position up to about 120 m NW where the SE mountain side is made up of coarse-grained sandstones which, although somewhat fractured, are nevertheless less altered from around 15 m deep downwards.

The situation changed towards the NW side where the presence of altered black shale caused lower V_P values. Nevertheless, at projected tunnel depths, velocities and hence material strength tended to increase (about 1,600 m/s in the extreme NW; Figure 8, under T3). The interblended shales were not very thick.

Conclusions

Geophysical refraction tomography was an important tool for evaluating rock rigidity, adding extra information regarding lithological characteristics. P velocity variation was determined at depths greater than 100 m in the survey area. This is the first time it has been done in Argentina, bearing the depth of investigation plus hard topographical conditions and refraction tomography in mind.

The Tamper system would appear to be a valuable tool for countless applications in the future, since it allows P and S waves to be generated at greater depths than the traditional hammer (its maximum range not having been studied yet) without the risk of using explosives. Moreover, its portability allows it to be carried into high mountain terrain, allowing the complete characterisation of rock masses, something which could not be performed by means of surface geophysical methods only, having to employ cross-hole studies which require expensive drilling and also provide data having much more bounded distances.

Acknowledgements

We would like to thank Juan Marcet, Director of "Escuela de Ingeniería de Caminos de Montaña" - at the Universidad Nacional de San Juan and his team for their tremendous cooperation in logistical tasks and positive attitude towards developing new geophysical methods, as well as their valuable contributions so that this daunting task could be carried out.

We would also like to thank Manuel Sánchez PhD, who was in charge of designing the Tamper emission system, for his continued willingness to cooperate with the project.

References

Bohachevsky, I. O.; Johnson, M. E. and Stein, M. L. (1986). Generalized simulated annealing for function optimization, Technometrics, 28, 209-217        [ Links ]

Bordonaro, O.L. (1980). El Cámbrico en la Quebrada de Zonda, Provincia de San Juan, Revista de la Asociación Geológica Argentina. 35, 26-40         [ Links ]

Cerveny, V. (2001) Seismic Ray Theory. Cambridge University Press, 700pp.        [ Links ]

Chunduru, R. K.; Sen M. K. and Stoffa, P. L. (1996) 2-D resistivity inversion using spline parameterization and simulatd annealing. Geophysics, 61, 151-161.        [ Links ]

Cortés, J; González Bonorino, G; Kourkharsky, M; Brodtkorb, A y Pereyra, F. (1997). Memoria y Mapa Geológico a escala 1:100.000 de la Hoja 3369-03 - Yalguaraz, provincia de Mendoza, Dir. Nacional de Servicio Geológico. Informe Inédito. Buenos Aires.        [ Links ]

Dobrin, M. (1975). Introducción a la Prospección Geofísica. Ed. Omega. España.        [ Links ]

Fernandez, A.L. (2000). Tomographic imaging the state of stress. PhD thesis, Georgia Institute of Technology. Atlanta. USA.        [ Links ]

González Díaz, E & Fauqué, L (1993) Geología y Recursos Naturales de Mendoza. V.A. Ramos (Editor). Relatorio XII Congreso Geológico Argentino. II Congreso Exploración de Hidrocarburos. Buenos Aires I-17:217-234        [ Links ]

Heiland, C.A. (1940). Geophysical Exploration, Prentice Hall. NY. 1013 pp.        [ Links ]

Herbst, R; Kapp, I; Krummel, H & Luck, E. (1998). Seismic sources for shallow investigations: A field comparison from northern Germany, Journal of Applied Geophysics, 38, 301-317.        [ Links ]

Imhof, A. L. (2009) Estudio Geofísico Ruta Nacional N°153. San Juan-Argentina. Informe Técnico. Escuela de Ingeniería de Caminos de Montaña (EICAM). Universidad Nacional de San Juan.        [ Links ]

Kirkpatrick, S., Gelatt Jr., C.D., Vecchi, M.P. (1983). Optimization by simulated annealing. Science 220, 671—680.        [ Links ]

Miller, R.; Pullan, S. E.; Waldner, J. S. & Haeni, F. P. (1986). Field comparison of shallow seismic sources, Geophysics 51, 2067-2092.        [ Links ]

Palmer, D. (1986). Refraction Seismics: The lateral Resolution of Structure and Seismic Velocity. Ed. Geophysical Press Ltd. London. ISBN: 0-946631-13-1. 269 pp.        [ Links ]

Pullammanappallil S.K. & Louie, J.N. (1994) A Generalized Simulated-Annealing Optimization for Inversion of First-Arrival Times, Bulletin of the Seismological Society of America, 84-5, 1397-1409        [ Links ]

Pullan, S.E.; MacAulay, H.A. (1987). An in-hole shotgun source for engineering seismic surveys, Geophysics, 52, 985-996.        [ Links ]

Schuster, G. T. & Quintus-Bosz, A. (1993). Wavepath eikonal traveltime inversion: Theory. Geophysics, 58, 1314-1323.        [ Links ]

Sethian, J. A. (2007). Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge University Press, 367 pp. ISBN 0521645573        [ Links ]

Vidale, J. E. (1988). Finite difference calculation of travel times, Bulletin of the Seismological Society of America, 78, 2062-2076.        [ Links ]

Zhang, J. & Toksoz, M.N. (1998). Nonlinear Refraction Traveltime Tomography, Geophysics, 63, no.5, 1726-1737.        [ Links ]