Characterization of karst hazards from the perspective of the doline triangle using GPR — Examples from Central Ebro Basin (Spain)

Engineering Geology 108 (2009) 225–236

Contents lists available at ScienceDirect

Engineering Geology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / e n g g e o

Characterization of karst hazards from the perspective of the doline triangle using GPR — Examples from Central Ebro Basin (Spain) Ó. Pueyo-Anchuela ⁎, A. Pocoví Juan, M.A. Soriano, A.M. Casas-Sainz Geotransfer Research Group, Departamento de Ciencias de la Tierra, Universidad de Zaragoza, C/Pedro Cerbuna, no. 12, CP.50.009 Zaragoza, Spain

a r t i c l e

i n f o

Article history: Received 9 August 2008 Received in revised form 21 June 2009 Accepted 29 June 2009 Available online 9 July 2009 Keywords: GPR Doline Sinkhole Karst Hazard 50 MHz antenna

A

b s t r a c t

Ground Penetrating Radar (GPR) is a widespread technique to locate and characterise karst hazards (cavities and paleocollapses). In this paper, we show that the internal structure of sediments obtained from GPR surveys can also be used as an indicator of active karst processes. A classification of karst hazard problems by means of GPR in the GPR-doline triangle is proposed, with three end-members: cavities, evidences of subsidence and paleocollapses. These end members show particular signatures in the GPR-profiles. The field examples shown in this paper indicate that GPR survey is a geophysical technique that offers a very high resolution and provides structural and sedimentological information of the subsoil. The use of grid maps, elaborated from GPR data, is an efficient way to determine anomalous sectors (lateral changes of electromagnetic properties, different reflectivity and qualitative penetration depth variations or velocity wave propagation changes). Their relationships with the structural features obtained from GPR-profiles (onlap geometries, zones with depressed reflectors, folded reflectors or laterally abrupt structural limits) can be used as indicators of karstic processes. © 2009 Elsevier B.V. All rights reserved.

1. Introduction The geophysical approach has been widely used in addressing karst hazards including different techniques such as: gravimetry (Bishop et al., 1997; Rybakov et al., 2001; Closson et al., 2003), seismic reflection, refraction and SASW (Steeples et al., 1986; Jansen et al., 1993; Lund and Korhonen, 1989; Debeglia et al., 2006), magnetometry (Rybakov et al., 2005; Mochales et al., 2006), electric tomography (Schoor, 2002; Zhou et al., 2002), seismic tomography (McDowell and Hope, 1994), and electromagnetic surveys including GPR (Collins et al., 1990, 1994; Robert and de Bosset (1994); Grandjean and Gourry, 1996; Freeland et al., 1998; Doolittle and Collins, 1998; Chang and Basnett, 1999; Chamberlain et al., 2000; Al-fares et al., 2002; Cunningham, 2004; Kofman et al., 2006). GPR (Ground Penetrating Radar) has been widely applied to locate and characterise karst problems during the past twenty years. However, most of these works have been focused in the location of underground cavities and the characterisation of paleocollapses. These can be understood as the main hazard sources in shallow karst areas. Nevertheless, the presence of active subsidence related to cavities below the GPR-penetration depth or dissolution processes not directly related with the development of cavities, can also represent an important hazard source. The presence of active subsidence processes without development of fully-grown cavities can be infer⁎ Corresponding author. E-mail address: [email protected] (Ó. Pueyo-Anchuela). 0013-7952/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.enggeo.2009.06.022

red by the characterization of the soil structure by the GPR survey. GPR shows several advantages over other geophysical methods, particularly in urban settings, since ubiquitous structural steel and vibrations from roadway or air traffic make it difficult to use gravimetry, terrain conductivity, magnetic, or seismic methods with an acceptable reliability (Mellett, 1995). The main goals for GPR surveys used in this work are to locate and characterise karst hazard zones in the areas surrounding the city of Zaragoza. Non-disturbed areas are represented by the general horizontal bedding of the Miocene and Quaternary deposits, and therefore, local thinning and thickening of the Quaternary deposits can indicate the presence of series sedimented within active karst areas. The initial hypothesis is that the internal structure of deposits can be used as an indicator of karstic processes, and the GPR-profiles can show evidences of active subsidence, presence of cavities and paleocollapses. In this work we present a new systematic methodology to be applied in GPR studies in potential or actual karst areas. The examples presented are taken from the Ebro Basin area in NE Spain, where karst risks are mainly associated with gypsum dissolution processes. 2. Geological setting The Ebro Basin is in the north-eastern part of Spain, bounded by the Pyrenees, the Iberian Range and the Catalonian coastal range. The Tertiary evolution of the basin was controlled by the development of the Pyrenean Range. During the Late Eocene the connection with the

226

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

Atlantic Ocean was closed and in the central sector of the Ebro Basin the Oligocene–Miocene series were related to alluvial fan systems, evaporitic and carbonate lacustrine systems (Pardo et al., 2004) that generated detrital, evaporite and carbonate rocks (Fig. 1). The Neogene rocks show near horizontal bedding with very gentle folds (Arlegui and Simón, 2001). During the Neogene, this endorheic basin opened to the Mediterranean Sea (Pérez-Rivarés et al. 2002; GarcíaCastellanos et al. 2003), thus changing its base level with subsequent dominating erosion processes and short sedimentation episodes. The Quaternary deposits of the Central Ebro Basin (gravels, sands and silt beds) were generated by fluvial (terraces) and alluvial systems (mantled pediments) mainly related with the Ebro river and its tributaries, and partially covering the Tertiary rocks. The Tertiary materials that outcrop in the surroundings of Zaragoza are mainly composed by anhydrite, secondary gypsum with marls, shales and mudstones. Borehole data in this area also indicate the presence of soluble salts as glauberite and halite. Karst processes in the surroundings of Zaragoza are related to the dissolution of Miocene evaporite materials partially covered by Quaternary materials. Karst activity is evidenced by different types of structures in most of the Quaternary fluvial levels, and is still active at present in the lower terraces and the alluvial plane. Dissolution processes cause intrastral dissolution, cavities, subsidence and collapses. The differences between surface morphologies resulting from karstic processes are mainly related to the cohesive behaviour of the Quaternary materials (Soriano and Simón, 1995). A cavity can form in the subsoil where a cap-rock permits to individualise the evolution of the dissolution zone from the uppermost materials. These cap-rocks can be limestones or gypsum Miocene levels and carbonate crusts (the so-called mallacán) in the Quaternary deposits. When these cap-rocks exist the presence of active karst processes is only observed after their sudden collapse. Other

situations occur when there is not a cohesive level and continuous subsidence is the dominant process. Intermediate cases appear when a thick alluvial level is present, the cavity has no expression at surface and deposits accommodate to a deep cavity (alluvial dolines). Intrastral dissolution produces accommodation and sagging of the Quaternary deposits. The surficial expression of intrastral or deep collapses are usually difficult to separate, the most common cases being intermediate cases between end-members (see for example, Pueyo Anchuela et al., 2009a, b). 3. Methodology Ground Penetrating Radar (GPR) is a non-destructive geophysical method that produces a cross-sectional image of the shallow subsurface – the resulting image (radargram) being very similar in style to seismic reflection profiles. GPR acquisition is based upon the propagation, reflection and scattering of high frequency electromagnetic waves in the soil (generally ranging from 10 to 1000 MHz). Reflections are produced by boundaries between elements with contrasting dielectrical properties, depending on grain size distribution (sorting, clay content), porosity, water content and the electrical properties of the particles themselves (Knight and Nur,1987; Davis and Annan, 1989; Annan et al., 1991). Consequently, subsurface features such as sedimentary structures, lithology changes and the water table can all generate primary reflections (Neal et al., 2002). Therefore, GPR results depend on the soil type and its conditions (saturation degree, compaction, mineralogy, etc.) and also on the frequency of the antennas used. In this work, a Ramac (Mala Geosciences, Sweden) GPR system with unshielded bistatic antennas and central frequency of 50 MHz was used (distance between antennas of 2 m). The trigger interval used was fired with time and the survey was aided with metric

Fig. 1. Geological map of Northeastern Spain with the location of the Ebro basin and the studied sites (1: Paleozoic, 2: Mesozoic, 3: Paleogene, 4: Neogene, 5: Quaternary).

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

measurements every 0.30 m and ticks in the profiles at fixed horizontal distances. The use of a 50 MHz equipment represents an acceptable resolution, in the range of 1 m, and a near surface penetration needed to usual conditions of alluvial karst hazards. The penetration depth obtained from the used equipment depends upon the conditions and characteristics of the subsoil, the different surveys carried out in this work present more penetration depth than the used TWT in the profiles, ranging from 15 to 30 m. On the other hand, the use of unshielded antennas makes that the presence of high conductivity elements or active sources of EM waves at/or near the surface, can have an important influence on the raw data. The surveys have been done in areas where none of these disturbations where present, or minimising their effect by making the survey parallel to walls and carefully eliminating all metallic elements present at surface before the survey. The time-depth conversion along the vertical axis (two-way time in ns) was performed using a constant propagation velocity estimated from the modelling of subsurface scattering diffractions (fitting hyperbolae) appearing on the raw data in each grid. The estimated velocities by hyperbola fitting are between 65 and 78 mm/ns. These velocities are representative for agricultural soil (Reynolds, 1997), what coincides with local stratigraphic descriptions of soil logs. The profile separation in each grid was 5 m, and the survey profiles followed two perpendicular directions, in order to obtain a three dimensional survey of the subsoil. This method allows for the lateral continuity of the anomalous sectors to be defined. The vertical resolution obtained by means of the wavelength and the velocity propagation of the waves is comprised between 0.3 and 0.4 m. The horizontal discrimination resolution by means of the first Fresnel zone gives values between 1.7 and 1.9 m (Reynolds, 1997) for a target located 1 m deep. These values are a mathematical approximation, and in the case of First Fresnel Zone they represent the distance at which two different objects with similar electromagnetic properties could be distinguished, the identification resolution being lower than this distance, although it also depends on the dielectric constant of the target and its depth. Profile processing included sequential zero-time corrections, running average filter, exponential gain function, band-pass filter and plotting. Interpretation was based upon the semi-quantitative evaluation of profiles, since in GPR profiles each individual measurement is represented by all the EM properties of the subsoil (conductivity and permeability, electric and magnetic); for example, the changes in the structure of the subsoil that can be easily observed in the profiles are related to changes in the lithology, mineralogy, state, clay content, iron oxides…. In the GPR profiles these changes can be identified from the change in color (in this case black and white) that is related to the amplitude of the recorded wave in each point. In this way, and not only in the profiles, the study of the changes in the wave-amplitude along one profile and different profiles are a semiquantitative assessment of the EM properties of the subsoil. In the same way that sedimentary structures are observed in the profiles, they can be observed by changes in the amplitude of waves, when these ones are analysed. Moreover, changes observed in the profiles sometimes are not evident in different profiles, because of the relationship between the structure observed and the surveying direction , or because the change in properties is not so sharp along certain directions. In this sense, the use of grid maps obtained from the lateral interpolation of the wave amplitude can reveal the lateral correlation between profiles of certain features that cannot be evident from the direct study of the profiles. Grid-maps can be drawn by lateral interpolation of profiles with the same direction in each grid. The lateral interpolation is made by the individualization of time and distance intervals in the profiles, analyzing the strength of the waveamplitude and comparing these results within the same profile or between different profiles. The discrete intervals defined depend on the expected resolution and the size of the element. In the same way

227

different depths (TWT intervals) can be analysed or residual maps can be done at high depths. The residual maps are elaborated where the strength of the wave is clearly attenuated and where the main observed changes are related with echoes of more surficial elements, making this maps without signal gain. The development of grid maps can be done directly from the raw data of GPR profiles in a xyz format (x: position in the profile, y: distance between profiles and z amplitude strength of the wave). Along the survey direction the interpolation depends more of the dimension of the looked element along the survey direction than on the needed degree of accuracy, whereas for the interpolation between different profiles the degree of definition depends of the mesh configuration. A usual practice in this kind of maps since the GPR wave has positive and negative values of the amplitude strength, is to use the square of the data (or absolute value) to avoid the usual horizontal disposition that can be observed in the profiles, whereas for map view, the change of the wave phase (more that the value itself) can be an indicator of change of the propagation velocity of the waves. These maps can be drawn either for the total vertical profile or drawing different grid-maps at different time intervals. Both types of maps were used to maximise anomalies linked to problematic sectors. The depth intervals to draw the grid maps were obtained from the GPR-profiles. The amplitude maps, or grid-maps, represent qualitative changes between different triggers and profiles, contours representing arbitrary units (depending on the processing and the magnitude of the changes). These maps were first applied systematically in Pueyo Anchuela et al. (2007) and have been used with reasonable success in karst environments (for example Chamberlain et al., 2000; Mochales et al. 2008; Pueyo Anchuela et al., 2009a, b). Correct profile processing is necessary to characterise the structure of the subsoil, and the previous routine has demonstrated an acceptable success. In the case of topographic irregularities corrections are needed. However, this correction was only necessary in several of the studied profiles, because slopes were usually less than 5°. Topographic correction should be considered in regions usually with surface gradients N10%, and slope angles N6° (Lehmann and Green, 2000). The stacking processing was only applied in profiles where vegetation growth does not allow a stable displacement. The grid map interpolation was done before and after the processing routine. The use of this kind of maps to locate anomalous sectors is based on the different apparent penetration (velocity propagation or reflectivity). Therefore, using these maps before applying gain filters maximizes the accuracy in locating anomalous sectors. In some cases, as in the case of profiles, the stacking (running average) filtering was applied to the grid maps. 4. Classification of karst hazards from evidences obtained from GPR Previous studies of karst realised in the Central Ebro Basin have been focused on their cross-section geometry, temporal evolution and structural damages (Soriano, 1990; Gutiérrez-Elorza and GutiérrezSantolalla, 1998; Soriano and Simón, 1995; Benito et al., 2000; Soriano and Simón, 2002; Guerrero et al., 2004; Soriano et al., 2004; Gutiérrez-Santolalla et al., 2005a and b and references therein). In the last years, geophysical methods have been applied to this subject for locating hazardous areas (Mochales et al., 2006). The main features of existing models are related to the presence of voids, the subsidence related to the absence of a cap-rock, the superficial collapse of a cavity, or the subsurface collapse of a cavity related to the Miocene materials or the Quaternary carbonate crusts with the subsequent genesis of an alluvial doline at surface. These features can be included within several groups of processes, all of them related to evaporite dissolution and underground lateral and vertical mass movements. Evidences for these processes can be simplified into three different groups: subsidence, cavities and collapses. These three

228

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

main terms or morphologies can be represented in the vertices of a triangle. They can be considered as the end-members of the doline triangle (Fig. 2): since the indicators used to classify the karst hazards are related to evidences inferred from the GPR, the triangle will be accordingly named the GPR-doline triangle. The end members of this classification are: Type 1. Cavities characterised by abrupt lateral contrasts in the electromagnetic properties of the materials. These features are represented by hyperbolic anomalies and are related to high resistivity zones. Depending on the shape and dimensions of these cavities it is possible to find either an isolated hyperbolic anomaly, or clustering of hyperbolic anomalies. The latter feature is related to wide cavities or roof irregularities of the cavity. The migrated profiles show zones with a change of the aspect of the radargram (usually higher apparent reflectivity). The presence of natural cavities is always found with partial filling collapse or associated with subsidence indicators. This can be related with the low competence of the cap-rock, that undergoes accommodation before collapsing. In Fig. 2 two usual cases in alluvial karst in the Zaragoza area have been represented. One of them is related to the development of cavities in the Miocene substratum and the second one in an area where the cohesion of carbonate crust allows the development of cavities in unconsolidated sediments. Type 2. Although subsidence itself cannot be measured by means of GPR-profiles, several indicators allow subsidence to be

inferred. These indicators are related to the internal structure of the soil. Presence of on-laps or wedge-like geometries of reflectors, progressive lateral variations of the electromagnetic properties of the subsoil and presence of depressed areas (inferred from reflector geometries associated with onlap reflectors in the upper part) can be used as indicators of subsidence processes. These features must be identified after applying topographic correction and migration of the profiles, and, in some cases, these features can be linked to sedimentary structures. In areas with active subsidence, GPRprofiles show two different geometrical features: a) The presence of local on-lap reflector patterns is the only indicator of active karst; in these cases, a closed envelope of accommodation features can be identified. The filling of the subsidence zone can be related to differential electromagnetic behaviour. The filling usually shows an increase of clay content with respect to the surrounding areas and therefore higher attenuation of the waves (negative anomalies in the grid maps) can be expected. Where subsidence zones are filled with urban debris, multiple reflection anomalies could be identified and positive anomalies in the grid maps were obtained. b) In other cases, there are more indirect indicators of subsidence processes, such as concave geometry of reflectors (deflected reflectors), linked to accommodation reflector geometries in the upper part of the profiles. The marginal zones of this on-lap reflector geometry show variations in the electromagnetic properties and are related to hyperbolic anomalies. These lateral variations of the electromagnetic properties and the presence

Fig. 2. GPR-doline triangle with the position of the different examples presented in this work.

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

of hyperbolic anomalies located at the margins of the subsidence zone are related to the presence of a sharp contact between the filling of the subsident zone and the host soil, and could be interpreted as the result of reverse faulting. In general, the identification of anomalous sectors related to subsidence can be done by interpretation of penetration variations in the grid maps (EM wave velocity changes or reflectivity variations) when the filling of the subsidence zone consists of materials different from the host rocks (different composition, grain size, structure, porosity, etc.). These features are usually related to on-lap geometry, folded reflectors and/or depressed zones in the GPR-profiles. The interpretation of the different subsident indicators is not straightforward, since they can be also present in non karst related environments. The geometrical relationships between reflectors in the side of a paleochannel for example, can be very similar to the expected adaptation features in the margin of a subsident zone. On the other hand, in the case of a paleochannel the EM characteristics of the filling can be also different from the host-sediment. In these cases, understanding of the 3D structure and its map view (lateral correlation of isolated indicators) can permit the correct interpretation of the geometrical relations identified in the GPR-profiles (for example closed envelopes, circular or elliptical, in map view). Type 3. Paleocollapses are one of the karst features that can be interpreted straightforward from GPR-surveys. Paleocollapses are characterised by sharp lateral contacts in the conductivity of the subsoil. In many cases, it is possible to identify the rupture surface as a reflector, when the properties of the collapse filling are different from the host rock. The EM wave behaviour in collapsed zones is enhanced when the filling consists of different materials, especially when sinkholes are filled with urban debris. Evidences associated with these structures are: i) lateral changes in EM behaviour, ii) presence of reflectors representing the marginal faults, iii) the loss of lateral reflector continuity (that can

229

indicate laterally truncated synclines) related to the collapse and iv) the structure of the filling (folded reflectors, depressed zones, reflectors with on-lap geometry,…). These three types of the GPR-doline triangle are end-members of a continuous phenomenon, common karstic features being intermediate cases. Therefore some aspects that can be ambiguous in interpretation of GPR profiles can be solved by considering the presence of more than one karst feature. 5. Case studies in the Zaragoza area In this section we present some examples of GPR surveys applied to locate sectors of anomalous behaviour of the EM waves, and its characterisation by the GPR-profiles in several karst hazard zones around Zaragoza. All the case studies come from test-sites, where some knowledge of the actual kart hazard was recorded by means of oral information, surface data or geological–geomorphological techniques. In some cases, structures inferred from GPR profiles were verified through direct observations in quarry outcrops. 5.1. Examples of cavities Isolated end-members of the GPR-doline triangle are not usually found. In natural environments, in nearly all the studied cases related to cavities, gentle subsidence could be identified. An endmember cavity is shown in Fig. 3, where two hyperbolic anomalies, whose apexes are located at depths of 5 and 8 m, can be observed. The horizontal or shallowly-dipping reflector is related to the contact between the Quaternary deposits and the Miocene limestones where the anomalies are located. These two hyperbolic anomalies are each related to cavity in the subsoil. The presence of cavities in the subsoil can be interpreted also from the grid-maps. The presence of a cavity can be inferred from an anomalous zone in the EM properties, but its interpretation as a cavity is not straightforward. In this case the cavities have a metric size and can be studied directly in

Fig. 3. (A) GPR profile showing two underground cavities. In the lower figures the hyperbolic anomalies have been interpreted (B) Anomaly correlation map (145 b t b 260 ns) with the location of the anomalous sector (in grey) with two example profiles from a cavity with irregular roof in the subsoil.

230

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

a near gully; these cavities have been historically used as dwelling places. A cavity shows EM properties different from the host-rock, as higher reflectivity, conductivity in some cases, or higher attenuation of the waves, and therefore the grid maps show a closed envelope of differential properties of the EM waves. The study of the internal structure by means of the GPR-profiles, can show hyperbolae anomalies or a clustering of hyperbolae anomalies, and the migrated profiles show zones with differential electromagnetic behaviour. In the profiles, the internal structure of the upper part is characterised by on-lap reflector geometries over the zone of clustering of hyperbolic anomalies (Fig. 3B). The interpretation of this anomalous zone is ambiguous: it can be related to a paleocollapse, a cavity with irregular roof shape or an intermediate case. The presence of a zone with different EM behaviour can be related to a cavity or an anomalous filling, and the existence of a hyperbolic anomaly is related to variations of the dielectric constant (cavity or materials different from the host rock). Cavities detected in other GPR-profiles (see for example Beres et al., 2001) show similar behaviour to those presented here. The presence of a cavity seems the most feasible interpretation and the on-lap reflectors geometry can be related to subsidence processes. On the other hand the microgravimetric survey over this section shows a decrease in the intensity of the residual anomaly whose modelling can be ascribed to a shallow cavity of metric scale (Mochales et al., 2008). 5.2. Examples of continuous subsidence In some cases the only evidence of internal structure of the sediments inferred from the GPR-profiles is the on-lap reflector geometry that indicates a process of accommodation to subsidence

processes. When this geometry shows a closed envelope and is related to a different electromagnetic behaviour of the materials in the enveloped zone, it can be interpreted as karst-related. In Fig. 4A, the GPR profile shows on-lap reflector geometries as the only indicator of internal structures different from horizontal layering (parallel reflectors). These on-lap geometries show a near-circular envelope in plan view. The grid-map shows a different penetration depth of the EM waves in the central zone of the grid. The interpolation map (Fig. 4A) shows no data in the central zone of the grid. Differences in the penetration depth occur near the on-lap reflector geometries, showing a plan-view geometry consistent with the location of the adaptation geometries. The presence of on-lap reflector geometries defining a circular structure is here interpreted as an indicator of differential subsidence (Fig. 4A). The different penetration depth can be explained by the differential electromagnetic properties between the filling of the subsident zone and the host rock. In this case the filling represents a wave-attenuation zone that can be interpreted to contain a higher proportion of clay. The survey grid was done along an agricultural field where the farmer indicated that he usually had to fill the surveyed zone with new soil every several years to level the topographic depression. In some cases, where there is the possibility to study the natural cross-sections, more evidences of karst can be observed. In the case shown in Fig. 4B, the filling of a depressed zone with a 10 m deep collapse in its central part can be observed in outcrop. The survey surface had been refilled and compacted by bulldozers. The presence of a collapse underneath can be identified in the outcrop, but not by GPR because of its limited penetration depth. The clay filling and its compaction can be the responsible of the high wave attenuation observed in this zone. In this case, the subsidence processes related with the collapse, but not the collapse itself, can be identified. In the

Fig. 4. (A) Anomaly correlation map (230 b t b 260 ns) with two example profiles of the subsident zone, the grey scale represents the change rate of the amplitude strength of the GPR wave. (B) GPR profiles with several on-lap reflector geometries related to the filling of a subsidence zone and direct comparison with an outcrop photograph where the survey was done (C) GPR profile with a central depressed zone associated with on-lap geometry reflectors and hyperbolic anomalies related with the boundaries of the on-lap reflector geometries.

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

profiles it is possible to identify accommodation features with opposite dips that could be correlated with the on-lap geometries observed in the outcrop. The marginal areas of the subsidence zones can also show hyperbolic anomalies. These anomalies can be related to the lateral contrast of EM properties between the filling and the host rock. In many cases, these marginal areas of subsidence zones can show nearvertical normal faults. In Fig. 4C there is a sector in the profiles with concave-upwards geometry (deflected reflectors), whose curvature decrease upwards and the upper reflectors show on-lap reflector geometries. The marginal zones of this deflected zone show a variation of the electromagnetic properties and are related to hyperbolic anomalies. They are best represented towards the right side of the profile. These lateral changes of the electromagnetic properties and the presence of hyperbolic anomalies located towards the end of subsidence zone, are related to the presence of a sharp contact between the materials of the subsiding zone and the host rock. Penetration tests in these zones indicated a cavity of metric scale in the left hyperbolic anomaly from Fig. 4C. 5.3. Examples of paleocollapses Paleocollapses can be directly detected by the very sharp limits of the collapse, the presence of deflected and folded reflectors, the loss in continuity of the reflectors and the variation of electromagnetic properties of the deposits. In grid-maps, a paleocollapse can be observed as an end-member of a subsident zone, where the propagation property variations between the filled zone and the host sediment are maximised. When the paleocollapse is reached by the GPR-survey, the complexity of its internal structure allows to classify structures and wave behaviour of different zones in the profile. In Fig. 5A three kinds

231

of geometry can be distinguished: (i) on-lap reflectors, (ii) folded reflectors without lateral continuity and (iii) a zone with high wave attenuation or homogeneus behaviour. Zone (iii) is surrounded by sharp boundaries covered by reflectors with on-lap geometry. The grid map shows an elliptical anomalous sector with a positive anomaly (Fig. 5B). This structure is related to a historical collapse filled with urban debris that can be identified from surface conditions (it is still active as can be assessed from the continuous filling by urban debris and the opening of ground cracks around the sinkhole). The zone with homogeneous behaviour has higher reflectivity and a qualitativelyestimated more penetration depth. This collapse is associated with folded reflectors cut by the collapse, subsidence and accommodation processes pre-dating the uppermost filling. In cases where the collapse happens suddenly, without any other evidence of karst activity, the only criterion of the collapse is the filling of the sinkhole, with an enhanced signal in GPR surveys when it consists of urban debris. In Fig. 6 the results of a GPR-survey over a filled historical collapse are presented. The collapse happened suddenly without any indicator at surface pre-dating the event and the pipe was subsequently filled several times with urban debris. The collapse shows a vertical pipe geometry 8 to 10 m deep with sharp, near-vertical walls and its base is more than 8 m deep. The GPR profiles show (Fig. 6) an anomalous sector defined by a higher penetration of the EM waves in the ground. The grid map shows a zone with a higher positive gradient that is related to multiple reflection anomalies and higher reflectivity targets in the profiles. The filling of the sinkhole consists of urban debris such as ceramic elements, bricks, iron electrical appliances, concrete, etc. These filling materials show in general high reflectivities. The distinct filling of the collapsed cavity is used to locate the paleocollapse and to characterise its filling from the structural point of view. The internal

Fig. 5. (A) GPR profile and its interpretation based on the defined classification. (B) Arial photograph with the studied historical sinkhole, grid maps obtained by the lateral correlation of parallel profiles (contour map and 3D shadowed relief maps of the contour map).

232

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

Fig. 6. Anomaly correlation map for Grid 1 (0 b t b 260 ns), with the location of the anomalous sector (in dashed line) and two example profiles from the grid. This survey is centered over a recent paleocollapse filled by urban debris, the photograph of 2005 was done before its filling and the GPR-survey.

structure of the filling shows on-lap reflector geometry that is related to its compaction. 5.4. Intermediate cases The most common case in karst areas, as pointed out in the preceding sections, is intermediate between the end-members of the GPR-doline triangle. In all the previously described cases there is a dominant end-member, whereas in the following examples mixed cases of the three end-members will be described. The location of the different case-studies in the GPR-doline triangle depends on the identification of the karst features from the GPR profiles and their contribution to the whole survey grid. In one of the quarries the GPR-profiles can be directly compared with outcrop photographs (Fig. 7). Several zones can be individualised in the profiles, according to their structure and electromagnetic behaviour (Fig. 7A and B). These zones are i) a homogeneous zone, ii) a zone with on-lap reflector geometry, iii) a zone without significant reflectors that could be associated with a high attenuation of the EM signal or a homogeneous structure, and iv) a zone with folded reflectors lacking continuity. This zone classification could be applied in the same sense in the case of the outcrop photograph (Fig. 7). Zones with homogeneous behaviour are related to horizontal structure, the zone with on-lap reflector geometry is related to the filling of a depressed area, with onlap geometries in sediments towards its borders. The zone with folded reflectors or without reflector continuity is related to stratigraphic levels that are folded and cut by collapse-related structures associated with the most subsident area. The zone with high wave attenuation is related to the higher clay content of the collapse filling. In this case, it is possible to establish a direct correlation between reflectors and stratigraphic levels (Fig. 7); although the lateral collapse contact cannot be identified as a reflector, probably because the rupture surface is filled by detrital materials, it can be located because of the different dip of reflectors. For other reflectors, it is possible to identify more detailed features in the GPR-profile than in the outcrop because GPR is very sensitive to the

water content in different levels (as a consequence of grain size variations, porosity, clay content, organic matter, iron oxides,…). These features cannot be directly observed in a photograph, where the layering can only be related to thin clay levels between sands and gravels. Folded reflectors can be interpreted as the result of bending before the final collapse. In some cases, the folding process is amplified after the collapse. In this situation, when reflectors showing on-lap geometries are not folded, pre-collapse folding can be interpreted. In this case the presence of a paleocollapse with indicators of subsidence before and after the collapse, and the absence of cavities or cavity-related anomalies, allow to locate it along the horizontal axis of the GPR-doline triangle (Fig. 2). The same zone classification can be applied in cases where no direct control of the outcrop is possible. A GPR survey was done in a doline field, where collapses and subsidence areas are present and a recent partial collapse was recorded. In the case shown in Fig. 8 three different anomalous sectors with circular to elliptical geometries and negative anomalies were identified. One of the anomalous sectors is related to the partial collapse and the others show differential vegetation growth and in some cases open cracks at surface in semicircular to semi-elliptical map view (see photographs in Fig. 8B). The collapses, with a diameter of about 2 m, permit to observe the shallow subsurface stratigraphy in this sector: the uppermost unit consists of gravels, usually with a carbonate crust, and without evidences of internal structure different from graded bedding. The GPR profile (Fig. 8) shows a surficial zone with homogeneous behaviour and on-lap reflector geometries in the central part of the profiles. The on-lap zones are related to penetration depth variations and inward-dipping reflectors. The central part of the GPR profile shows a zone without reflections or with homogeneous behaviour. In this case, all the anomalous sectors are related to shallower penetration of the survey that can be interpreted according to the presence of more attenuating material in these zones. There are not evidences of cavities in the surveyed profiles, what can be explained either because they do not exist or because they are located below the maximum penetration depth of the GPR. These anomalous sectors can be interpreted as three aligned alluvial dolines, with accommodation synclines. Wave

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

233

Fig. 7. (A) GPR profile with the structural interpretation based on the classification of structural zones (B) Outcrop photograph with the same structural zone classification applied in (A). (C) Outcrop photograph with its interpretation and GPR profile of the same zone.

attenuation in the central part of the profiles, related to the collapsed zone, can be interpreted as the filling of a natural collapse with higher clay content than the host rock. Therefore, in the grid maps, anomalies are negative, in contrast with the positive anomalies associated with urban debris filling. The open cavity may represent the vertical propagation of a cavity involving the Miocene substratum, what is compatible with the absence of cavities in the profiles within the surveyed depth. The on-lap geometry filling does not show folded geometries (accommodation signatures), what means that the subsidence happened before the beginning of the filling. 6. Discussion The electromagnetic properties of dolines or subsidence zone filling parameters allow for anomalous sectors to be located. The study of these anomalous sectors within the GPR-profiles can help to improve and reduce the ambiguity of some interpretations obtained from other methods. The integrated study of grid maps and profiles, and the particular features of each GPR-doline type, is the clue to the GPRsurvey success.

The simplification of karstic hazard structures in the GPR-doline triangle, summarised in cavities, evidences of subsidence and paleocollapses, allows for most of the studied examples to be classified according to their properties. Each category shows distinct signatures in the GPR-profiles. The grid maps show differences in the penetration depth of the GPRsurvey. These differences can be related to penetration depth variations, changes in the EM wave propagation velocity or changes in the reflectivity of the subsoil. These parameters allow the three end-members to be characterised. Their main features can be identified from the study of the subsoil internal structure inferred from the GPR-profiles. Several features that could be interpreted as karstic (in the sense used in this work) can be ascribed to sedimentary structures not related to karst development. The study of geophysical techniques together with outcrop study and the combination of several karstic features will contribute to the success in identification of possible hazards. The presence of on-lap reflector geometries is ambiguous, since they can be either the result of subsidence processes (long-term or short-term) or they can just be related to depositional processes. The association with closed depressed areas, folded reflectors and the loss

234

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

Fig. 8. (A) Anomaly correlation map with location of profile A, anomalous sector found and location of the EM grid. (B) 3D shadowed relief map of the grid and photographs from a surficial collapse that correlates with the most Southern anomaly and open cracks identified surrounding the central anomaly. (C) Plot of apparent conductivity and susceptibility measured over the Profile A (Gem-2 multifrequency EM device) and map view from the EM-Grid over the central anomaly identified in the GPR-grid maps and interpretation of profile A based on the defined classification shown in Fig. 7A. See profile location in panel A.

in continuity of reflectors are unequivocal features of karst processes. These features bring about analogous identification problems in the grid maps, when only one of these criteria is present. The presence of an anomalous sector with a positive penetration anomaly in a collapsed zone is related to urban debris fillings. Negative anomalies are associated with natural fillings with higher clay content. The geometry of these anomalous sectors can help to its interpretation, but without structural indicators, their understanding in terms of karst hazard is not straightforward. In the case of cavities, similar ambiguities can be defined. The presence of hyperbolic anomalies is related to dielectric variations, and high permitivities are related to cavities. However, dielectric variations are only qualitative indicators that can be related to other natural targets. The collapse end-member (type 3) is possibly the easiest recognizable karstic feature, because of the presence of several indicators in

the same place (several structural indicators and anomalous sectors in the grid maps). 7. Conclusions Karstic features detected by GPR can be summarised in the conceptual GPR-doline triangle with three end members: cavities, subsidence and paleocollapses. In karst areas common features reveal intermediate cases between the three end-members of the GPRdoline triangle, as shown in this work. The examples presented here are mainly intermediate cases in the GPR-doline triangle, and the integration of different features of the same process allows to locate, characterise and identify the karst problems. The interpretation of the GPR features shows that the GPR surveying represents a good and proper geophysical technique to be applied in problematic karst zones, even where cavities are below the

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

maximum penetration depth of the survey. One of the advantages of GPR surveys against other methods is that GPR does not need modelling of the karst structure to determine the origin of anomalies, because the GPR-profiles show directly the complicated internal structure of the karstic subsoil. The use of grid maps represents a fast way to determine anomalous sectors (lateral changes of electromagnetic properties, different reflectivity, and qualitative penetration depth variation or velocity wave propagation) and their relationship with the structural features obtained from the GPR-profiles (on-lap geometries, depressed reflectorszones, folded reflectors or lateral abrupt structural limits) that can be considered direct evidences of karst activity. The application of this methodology to the surroundings of Zaragoza is favoured by the relatively simple structure of the deposits affected by karst: widespread bending of the cover materials adapting to the underground cavities or the usual practice of filling the collapsed and the subsident zones with urban debris materials. These materials have a high reflectivity and therefore can produce multiple reflection anomalies. They represent qualitative changes in the wave velocity propagation and penetration depth that allow for anomalous sectors in urban areas to be easily defined. The integration of GPR-surveys with classical detection techniques such as the study of aerial photographs and farmer interviewing, can help to date karst activity and therefore to establish the relative chronology between processes inferred from the GPR survey. Acknowledgements This work has been financed by the Geotransfer Research Group (Zaragoza University) and a predoctoral grant to the first author (Ministerio de Educación). The authors want to thank the work of Dr. R. J. Shlemon (editor) and three different anonymous reviewers that have improved significantly this work. References Al-fares, W., Bakalowicz, M., Guérin, R., Dukhan, M., 2002. Analysis of the karst aquifer structure of the Lamalou area (Hérault, France) with ground penetrating radar. Journal of Applied Geophysics 51 (2–4), 97–106. Annual International Meeting Program with AbstractsAnnan, A.P., Cosway, S.W., Redman, J.D., 1991. Water table detection with ground-penetrating radar. Society of exploration geophysicists 494–497. Arlegui, L.E., Simón, J.L., 2001. Geometry and distribution of regional joint sets in a nonhomogeneus stress field: case study in the Ebro basin (Spain). Journal of structural geology 23, 297–313. Benito, G., Gutiérrez, F., Pérez-González, A., Machado, M.J., 2000. Geomorphological and sedimentological features in Quaternary fluvial systems affected by solution-induced subsidence (Ebro Basin, NE-Spain). Geomorphology 33, 209–224. Beres, M., Luetscher, M., Olivier, R., 2001. Integration of ground-penetrating radar and microgravimetric methods to map shallow caves. Journal of Applied Geophysics 46, 249–262. Bishop, I., Styles, P., Emsley, S.J., Ferguson, N.S., 1997. The detection of cavities using the microgravity technique: case histories from mining and karstic environments. Geol Soc Eng Geol Spec Publ 12, 153–166. Chamberlain, A.T., Sellers, W., Proctor, C., Coard, R., 2000. Cave detection in limestone using ground penetrating radar. Journal of Archaeological Science 27 (10), 957–964. Chang, K.R., Basnett, C.,1999. Delineation of sinkhole boundary using Dutch cone soundings. Engineering Geology 52 (1–2), 113–120. Closson, D., Karaki, N.A., Hussein, M.J., Al-Fugha, H., Ozer, A., Mubarak, A., 2003. Subsidence et effondrements le long du littoral jordanien de la mer Morte : apports de la gravimétrie et de l'interférométrie radar différentielle: Subsidence and sinkholes along the Jordanian coast of the Dead Sea: contribution of gravimetry and radar differential interferometry. Comptes Rendus Geosciences 335 (12), 869–879. Collins, E., Cum, M., Hanninen, P., 1994. Using ground-penetrating radar to investigate a subsurface karst landscape in north-central Florida. Geoderma 61 (1–2), 1–15. Collins, M.E., Puckett, W.E., Schellentrager, G.W., Yust, N.A., 1990. Using GPR for microanalyses of soils and karst features on the Chiefland Limestone Plain in Florida. Geoderma 47 (1–2), 159–170. Cunningham, K.J., 2004. Application of ground-penetrating radar, digital optical borehole images, and cores for characterization of porosity hydraulic conductivity and paleokarst in the Biscayne aquifer, southeastern Florida, USA. Journal of Applied Geophysics 55 (1–2), 61–76. Davis, J.L., Annan, A.P., 1989. Ground penetrating radar for high resolution mapping of soil and rock stratigraphy. Geophysical prospecting 37, 531–551.

235

Debeglia, N., Bitri, A., Thierry, P., 2006. Karst investigations using microgravity and SASW: application to Orléans, France. Near Surface Geophysics 4, 215–226. Doolittle, J.A., Collins, M.E., 1998. A comparison of EM induction and GPR methods in areas of karst. Geoderma 85 (1), 83–102. Freeland, R.S., Yoder, R.E., Ammons, J.T., 1998. Mapping shallow underground features that influence site-specific agricultural production. Journal of Applied Geophysics 40 (1–3), 19–27. García-Castellanos, D, Vergés, J, Gaspar-Escribano, J.M., Cloething, S., 2003. Interplay between tectonics, climate and fluvial transport during the Cenozoic evolution of the Ebro Basin (NE Iberia. Journal of Geophysical Research 108 (B7), 2347. Grandjean, G., Gourry, J.C., 1996. GPR data processing for 3D fracture mapping in a marble quarry (Thassos, Greece. Journal of Applied Geophysics 36 (1), 19–30. Guerrero, J., Gutiérrez, F., Lucha, P., 2004. Paleosubsidence and active subsidence due to evaporite dissolution in the Zaragoza area (Huerva River valley, NE Spain): processes, spatial distribution and protection measures for transport routes. Engineering Geology 72 (2004), 309–329. Gutiérrez-Elorza, M., Gutiérrez-Santolalla, F., 1998. Geomorphology of the Tertiary gypsum formations in the Ebro Depression (Spain). Geoderma 87, 1–29. Gutiérrez-Santolalla, F., Gutiérrez-Elorza, M., Marín, C., Desir, G., Maldonado, C., 2005a. Spatial distribution, morphometry and activity of La Puebla de Alfindén sinkhole field in the Ebro river valley (NE Spain): applied aspects for hazard zonation. Environmental Geology 48 (3), 360–369. Gutiérrez-Santolalla, F., Gutiérrez-Elorza, M., Marín, C., Maldonado, C., Younger, P.L., 2005b. Subsidence hazard avoidance based on geomorphological mapping in the Ebro River valley mantled evaporite karst terrain (NE Spain). Environmental Geology 48 (3), 370–383. Jansen, J., Anklam, J., Goodwin, C., Roof, A., 1993. Electromagnetic induction and seismic refraction surveys to detect bedrock pinnacles: Proc 4th Multidisciplinary Conference on Sinkholes and the Engineering and Environmental Impacts of Karst, Panama City, Florida, P115-122. Publ Rotterdam: A A Balkema. International Journal of Rock Mechanics and Mining Science and Geomechanics Abstracts 30 (6), 362 December. Knight, R.J, Nur, A.,1987. The dielectric constant of sandstones, 60 Khz to 4 Mhz. Geophysics 52, 644–654. Kofman, L., Ronen, A., Frydman, S., 2006. Detection of model voids by identifying reverberation phenomena in GPR records. Journal of Applied Geophysics 59 (4), 284–299. Lehmann, F., Green, A.G., 2000. Topographic migration of georadar data: implications for acquisition and processing. Geophysics 65 (3), 836–848. MarchLund, C., Korhonen, R., 1989. Karst cavity detection by means of seismic tomography to help assess possible dam site in limestone terrain Proc Symposium on Rock Mechanics and Power Plants, Madrid, 12-16 Sept 1988 V1, P25-31. Publ Rotterdam. International Journal of Rock Mechanics and Mining Science and Geomechanics Abstracts 26 (2), 72. McDowell, P.W., Hope, V., 1994. Location and delineation of karst and dissolution collapse features by acoustic tomography : Proc 4th Multidisciplinary Conference on Sinkholes and the Engineering and Environmental Impacts of Karst, Panama City, Florida P123-129. Publ Rotterdam: A A Balkema. International Journal of Rock Mechanics and Mining Science and Geomechanics Abstracts 31 (1), A28 February. Mellett, J.S., 1995. Ground penetrating radar applications in engineering, environmental management, and geology. Journal of Applied Geophysics 33, 157–166. Mochales, T., Casas, A.M., Pueyo, E.L., Pueyo, O., Román, M.T., Pocoví, A., Soriano, M.A., Ansón, D., 2008. Detection of underground cavities by combining gravity, magnetic and ground penetrating radar surveys: a case study from the Zaragoza area, NE Spain. Environmental Geology 53 (5), 1067–1077. Mochales, T., Pueyo, E.L., Casas, A.M., Soriano, M.A. 2006. Magnetic prospection as an efficient tool for doline detection. A case study in the central Ebro Basin (Northern Spain). Geological Society of London. Special Publications: Natural and Anthropogenic Hazards in Karst: Recognition, Analysis and Mitigation. Neal, A., Pontee, N.I., Pye, K., Richards, J., 2002. Internal structure of mixed-sand-andgravel beach deposits revealed using ground-pentrating radar. Sedimentology 49, 789–804. Pardo, G., Arenas, C., González, A., Luzón, A., Pérez, A., Pérez-Rivarés, F.J., Vázquez-Urbez, M., Villena, J., 2004. La cuenca del Ebro. In: Vera, J.A. (Ed.), Geología de España. SGEIGME, pp. 533–543. Pérez-Rivarés, F.J., Garcés, M., Arenas, C., Pardo, G., 2002. Magnetocronología de la sucesión miocena de la Sierra de Alcubierre (sector central de la Cuenca del Ebro). Revista de la Sociedad Geológica de España. 15 (3–4), 217–231. Pueyo Anchuela, Ó., Pocoví Juan, A., Soriano Jiménez, M.A., Casas Sainz, A.M., Mochales López, T., 2007. GPR as a tool to detect and characterize subsidence and collapse associated with shallow karst. Examples from the Central Ebro Basin (Spain. Geophysical Research Abstracts 9, 08911. Pueyo Anchuela, Ó., Casas-Sainz, A.M., Soriano, M.A., Pocoví Juan, A., 2009a. Mapping subsurface karst features with GPR: results and limitations. Environmental Geology. doi:10.1007/s00254-008-1603-7. Pueyo Anchuela, Ó., Pocoví Juan, A., Casas-Sainz, A.M., Soriano, M.A., Ansón López, D., 2009b. Radarfacies and sinkhole hazards. Geophysical Research Abstracts, Volume 11, 8116. Reynolds, J.M., 1997. An Introduction to Applied and Environmental Geophysics. John Wiley and Sons, Chichester. Robert, A., de Bosset, C., 1994. Application du géoradar à la localisation de cavités, de nids de gravier et de zones karstiques. Journal of Applied Geophysics 31 (1–4), 197–204. Rybakov, M., Goldshmidt, V., Fleischer, L., Rotstein, Y., 2001. Cave detection and 4-D monitoring: a microgravity case history near the Dead Sea. Leading Edge (Society of Exploration Geophysicists) 20 (8), 896–900. Rybakov, M., Rotstein, Y., Shirman, B., Al-Zoubi, A., 2005. Cave detection near the Dead Sea—a micromagnetic feasibility study. The Leading Edge (Society of Exploration Geophysicists) 24 (6), 585–590.

236

Ó. Pueyo-Anchuela et al. / Engineering Geology 108 (2009) 225–236

Schoor, M., 2002. Detection of sinkholes using 2D electrical resistivity imaging. Journal of Applied Geophysics 50 (4), 393–399. Soriano, M.A., 1990. Geomorfología del sector centro-meridional de la Depresión del Ebro. Diputación Provincial de Zaragoza. Soriano, M.A., Simón, J.L., 1995. Alluvial dolines in the central Ebro Basin, Spain: a spatial and developmental hazard analysis. Geomorphology 11, 295–309. Soriano, M.A., Simón, J.L., 2002. Subsidence rates and urban damages in alluvial dolines of the central Ebro Basin (NE Spain). Environmental Geology 42 (5), 476–484. Soriano, M. A., Simón, J. L., Arlegui, L. E., Liesa, C. L., Pocoví, A., 2004. Problemas causados por el karst aluvial en el centro de la cuenca del Ebro (España). Stato Dell'arte sullo

studio dei fenomeni di sinkholes e ruolo delle amministrazioni statali e locali nel governo del territorio, 521–532. Steeples, D.W., Knapp, RW., McElwee, CD., 1986. Seismic reflection investigations of sinkholes beneath interstate Highway 70 in Kansas Geophysics V51, N2, P295–301. International Journal of Rock Mechanics and Mining Science and Geomechanics Abstracts 23 (5), 188. Zhou, W., Beck, B.F., Adams, A.L., 2002. Effective electrode array in mapping karst hazards in electrical resistivity tomography. Environmental Geology 42, 922–928.