3-D GPR visualization technique integrated with electric resistivity tomography for characterizing near-surface fractures and cavities in limestone

Abstract

Near-surface geophysical techniques, (i.e. ground-penetrating radar (GPR) and electric resistivity tomography (ERT)), continue to enable great possibilities for 3-D visualization of fractures and cavities. An integrated survey, including GPR and ERT, was used for mapping and evaluating a complex subsurface fracture network in Eocene limestone at Al-Mokattam site, Egypt. While only GPR was available to investigate an engineering site of about 200 × 300 m, located in Riyadh city, Saudi Arabia, for cavity detection. At the Egyptian site, GPR data were acquired with a GSSI SIR2000 system using a 400 MHz shielded antenna with 0.2 m in-line spacing. GPR data visualized using an isosurface rendering combined with inverted ERT allowed us to identify and follow the fracture system, which was mainly filled with low-resistivity sediment; such fractured limestone may cause landslides. GPR data at the Riyadh site were acquired using the RAMAC system with a 100 MHz shielded antenna attached to a GPS with real-time kinematic total station equipment; GPR data confirmed the possible presence of cavities at different depths, which may cause severe damage during building construction. In this paper, we demonstrated the efficiency of each geophysical technique available for solving subsurface engineering problems, leading to future improvements in site safety. Additionally, we indicated the efficiency of an advanced 3-D GPR visualization technique (isosurface rendering) for improving fracture and cavity detection, characterization and interpretation.

KEYWORDS:

• GPR
• ERT
• isosurfacerendering
• fracture
• cavitydetection
• engineering-geophysics

1. Introduction

Ground-penetrating radar (GPR) is a high-resolution, non-destructive technique used for accurately probing near-surface structure. It is based on the principle that high-frequency electromagnetic (EM) waves may be reflected or refracted at subsurface features characterized by changes in their electrical properties [1]. The GPR technique can reveal interfaces between different materials, if there is a sufficient contrast between the dielectric properties of these materials. The advantages of GPR surveys, especially for mapping and imaging subsurface geological features are well documented [2–7]. The high horizontal and vertical resolution required for the detection of fractures and cavities makes GPR an effective technique for large-area surveys [8–10]. Therefore, GPR is best suited for high-resolution geophysical and geological engineering investigations, particularly if it is combined with another geophysical technique such as electrical resistivity tomography (ERT). Defects in the subsurface rocks, such as fractures or cavities, create a challenge for imaging, especially fracture networks can produce an abundance of reflections, displacements, and hyperbolic diffractions in GPR profiles (Grasmueck et al. [8]). GPR and ERT have been widely applied in the last decades for fracture and cavity detection (e.g. [11–18]), while 3-D GPR visualization has been widely used for active fault imaging (e.g. [4,19–21]).

GPR anomalies can be easily identified and correlated throughout different adjacent profiles in the acquisition grid, from which we can follow the spatial extension of the target of interest. Conventional GPR used in subsurface exploration is mainly based on relatively widely spaced, parallel 2-D survey profiles with in-line spacing ranging from 0.5 to 1 m (0.5 m is the common in-line spacing in many case studies), with extensive interpolation used to fill in the shortage of data between adjacent profiles. The resulting subsurface images can be displayed as fence diagrams or highly interpolated pseudo 3-D image cubes [22,23].

Conventional time-slice presentation is the most popular and the most rapid method used to display 3-D GPR images from the relatively wide-spaced series of contiguous 2-D parallel profiles, in which amplitude reflections are mapped at a constant two-way travel time [24]. However, time-slice imaging can be considered as a background filter, because GPR recorded traces can have significant horizontal banding noise caused by equipment design as well as reverberations from the ground surface [25].

Horizontal time-slice, in many cases, may not be the suitable choice for subsurface a visualization tool in cases of great subsurface complexity such as structures buried or extended at different depths, and/or anomalies with an irregular morphology [26]. To improve the efficiency and quality of subsurface target images, it is very important to use more advanced 3-D visualization techniques (e.g. the isosurface rendering technique), which has been considered a successful tool for such purpose.

Accurate subsurface imaging of the spatial variability fractures and cavities is essential when investigating land suitability for construction, and/or land use management based on limited subsurface soil information. GPR data analysis combined with ERT is a crucial technique to efficiently prioritize areas with cost-effective optimization of construction times, to efficiently implement hazard mitigation measures, and to design projects with appropriate and reliable foundation systems to compensate for geotechnical risks (e.g. [27–29]). However, existing geotechnical investigations (traditional engineering techniques) that are widely performed to understand the potential problems of soils can be time-consuming, expensive, and limited. The limitations may include required depth of investigation, subsurface soil materials, etc.

Al-Mokattam plateau (Figure 1(a)) represents a serious danger problem, especially at the borders of the plateau, where at the southern part of the Al-Mokattam area there have been many landslides due to the fractured limestone. The landslides have caused damage to many houses and hotels around the area. The fractured limestone is underlain by a clay layer that causes most of the engineering problems because of volume changes in swelling clays resulting from human activities that modify the local environments [30]. In the past few years, several rock failures have damaged houses, roads, and bridges, and caused loss of human life. For example, in September, 2008, a rockslide hits the nearby houses, killing 24 people and injuring over 300 persons [31].

Figure 1. Google earth location maps for the geophysical survey sites. (a) Al-Mokattam survey site in Egypt, where ERT and GPR carried out, while (b) the location for Riyadh investigated site, KSA, where GPR was the only available technique used in the site.

The Riyadh survey site (Figure 1(b)), is mainly characterized by limestone bedrock of the Sulaiy Formation that is considered as a critical hazard for construction safety because of fractures and cavities that occur in this limestone rock. These cavities developed in limestone as a result of the chemical leaching of the carbonate and evaporite formations by infiltrating water from nearby urban areas [16,32]. Section #2 will discuss in detail the geologic and geomorphologic setting of these two sites.

An objective of this study was to evaluate the potentiality of near-surface geophysical techniques available (GPR and ERT) to image out the fracture system running in the limestone plateau at the Al-Mokattam site, Cairo, Egypt, and to apply GPR survey to detect and define the shape and size of cavities in the limestone construction site in the Riyadh area, KSA. This research article introduces an application of 3-D visualization technique using isosurface rendering that is presented along with a comparison to conventional time-slice presentation. The results of this comparison have revealed an implication for easily imaging and tracking of subsurface fractures and cavities comparing with conventional methods used; herein we outline the main advantages of these isosurface techniques for target identification and characterization. We believe that results obtained in this paper will help to solve many future engineering geologic problems in complex situations or in challenging case studies like subsurface fracture and/or cavities.

2. Geological and geomorphological setting of the studied areas

2.1. Al-Mokattam site geological setting

The Al-Mokattam survey site is located between latitudes 29° 59′ 20″ and 30° 00′ 50″ and longitudes 31° 16′ 30″ and 31° 19′ 00″. The surface geology of the study area was described by the Geological Survey of Egypt [33,34] and the National Authority of Remote Sensing and Space Sciences [35]. The area consists of two Upper Eocene formations. The first is the Maadi Formation, which consists of soft clastic rocks (clay, silt and sand) and hard dolomitic limestone; the second is the Giushi Formation, which consists of white fossiliferous limestone with marl intercalation.

Geomorphologically, the Al-Mokattam area consists of three plateaus. The upper plateau lies in the northern part and has high elevations ranging from 205 m in the south to about 170 m in the north; it is separated from the middle plateau by a steep cliff. The middle plateau has an elevation ranging from 110 to 150 m. The lower plateau has an elevation ranging from 50 to 80 m [35]. The slopes of the Middle Eocene rocks that form the boundaries of the middle plateau of Gabal Al-Mokattam are unstable because of dissected fractured limestone with clay intercalation and represent potential areas of rock failure in many places [36] (Figure 2).

Figure 2. Diversity of fractures running in Al-Mokattam limestone Formation at Jabal Al-Mokattam Area.

The foundation rock of the Al-Mokattam upper plateau is highly jointed, and the slopes of the south and south-west part have been affected by slope failures where several buildings were subjected to severe damage. In February, 2019 blocks of limestone fell down on the urban area at the northern part of the plateau (the Manshiyat Naser area just below Al-Mokattam hills), and more than 100 persons were injured according to local newspapers. GPR and ERT surveys were used to investigate the subsurface situation and to image fracture continuity, to avoid future rock collapse and to help preserve people’s life and property.

2.2. Riyadh site geological setting

The Riyadh area is located on the Najd plateau, in the central part of the Kingdom of Saudi Arabia, between longitude 24° 38' North and latitude 46° 43' East with an elevation of about 600 m above the mean sea level. Most of Riyadh city was built over a subsurface of carbonate rocks of the Jubaila formation from the Upper Jurassic era [37]. The geology of the Riyadh site is mainly characterized by two different limestone bedrock units, both of which have high potential of developing karstic cavities, sinkholes, and open fractures [38]. The first limestone unit is a NW-SW limestone belt on the eastern side of the Riyadh site and is known as the Sulaiy Limestone Formation of the Cretaceous age [39]. The Sulaiy Formation is typically composed of compacted limestone with a few thin calcareous beds. On outcrops, this formation exhibits slumping features in its lower beds identical to the features found in the Arab Formation. The higher beds, however, are unaffected by slumping and are moderately strong, forming erosion-resistant, well-defined steep scarp slopes. Subsurface cavities and sinkholes are most frequently created in the lower limestone beds of the Sulaiy Formation at the contact with the Arab Formation rather than in the upper beds as they are (both the Sulaiy Formation and the Arab Formation) formed mainly of compacted limestone. The second karst unit belongs to the Jurassic limestone of the Arab Formation. This formation outcrops as an NW-SE trending of limestone belt at the western part of Riyadh, where our survey site is located (Figure 1(b)). This formation is mainly composed of limestone and evaporites, and it is considered as the most susceptible bedrock to develop soluble cavities in carbonate rocks. This could be attributed to either the dissolution of evaporites or to the nature of the limestone, which is probably composed of limestone boulders within an extremely soluble matrix of softer dolomitic limestone [38]. In some locations, especially at the eastern part of Riyadh city, the highly fractured Sulaiy limestone formation is directly overlaid by the Arab Formation, due to the dissolution of the anhydrites of the Hith Formation [40]. Solution-collapse of soft and porous breccia in the slumping carbonate beds is the main result of dissolution in the near-surface anhydrites [41,42] (Figure 3).

Figure 3. (a) Typical subsurface cavities developed in limestone at Riyadh construction site, KSA. (b) The size and dimension of cavities in limestone bedrock of Sulaiy formation.

Such subsurface geological conditions at the Riyadh site, along with surface human activities, affect the ability of limestone to withstand construction and leave it significantly hazardous that may cause land subsidence, differential settlements and/or the collapse of buildings [43].

3. Geophysical data acquisition

GPR data at the Al-Mokattam site was acquired with a SIR-2000 system from Geophysical Survey Systems Inc. (GSSI) using the 400 MHz shielded antenna (Figure 4). Data were obtained using an odometer survey wheel with 0.2 m in-line spacing and 0.02 m cross-line spacing. GPR scan traces were recorded at every 2 cm with 512 samples/scan using a sampling rate interval of 0.2344 ns producing a two-way travel time with good reflection and penetration to about 120 ns corresponding to a depth of about 5 m based on a subsurface velocity of 0.1 m/ns. The scan area was about 8 m x 21 m at the cliff of a fractured limestone wall exposure (Figure 4(b)).

Figure 4. (a) Location map of the survey site. (b) The 400 MHz shielded antenna from GSSI used for acquisition of the relatively dense GPR data set at Al-Mokattam survey site in Egypt. String was used as guidelines for antenna movement, while measurement tape was used to place GPR antenna at each profile. (c) ERT measurement using Syscal-R2 multi-electrode system.

At the Al-Mokattam site, 2-D geoelectric resistivity survey was available and was used to investigate the subsurface resistivity distribution and geological setting in the site, to help interpretation of 3-D GPR data and image subsurface soil moisture content in the soil (Figure 4). Electrical resistivity tomography (ERT) data were acquired using a SYSCAL-R2 multi-electrode system attached with 80 electrodes (connected to a network of intelligent nodes, which was used as an automatic multi-electrode switching system) using a Wenner-Schlumberger array, with 1.5 m electrode spacing, providing a depth of about 23 m with full subsurface data coverage to reach the expected target of interest for each study area. Each ERT measurement cycle required about 1.7 s, where the delay time was about 0.3 s. The intensity of the current injection during data acquisition was between about 150 and 230 mA (Figure 4(c)).

GPR data at the Riyadh site were acquired using a RAMAC system from MALA Geoscience with 1a 00 MHz shielded antenna attached to GPS real-time kinematic (RTK) total station equipment for position calculations (Figure 5). The trace spacing was 0.05 m, which was used to determine the horizontal resolution. The total recorded time was about 311 ns, where this time is used for depth calculation, corresponding to the depth of about 15 m based on a subsurface velocity of 0.1 m/ns. Positioning data with RTK was acquired every 1 m using high-accuracy differential GPS (TRIMBLE GPS 5800 with TSC), where the x-coordinate and y-coordinate precision is about 1 cm and the z-coordinate precision is about 2 cm. Table 1 presents the acquisition parameters for the RAMAC-100 MHz and the GSSI-400 MHz shielded antenna, GPR systems, respectively.

Figure 5. GPR data acquisition with 100 MHz shielded antenna using RAMAC system from MALA Geoscience. RTK was used to locate antenna positioning with high-accuracy differential GPS (TRIMBLE GPS 5800 with TSC).

Table 1. Acquisition parameters for the RAMAC-100 MHz and the GSSI-400 MHz shielded antennas, GPR systems.

Nominal antenna frequency (MHz)	RAMAC 100 MHz	GSSI 400 MHz
Antenna configuration	Monostatic mode	Monostatic mode
Number of samples per scan	360	512
Sampling frequency (ns)	1083.92	–
Range (trace length in ns or time window)	311 ns	120ns
Trace interval	0.05 m	0.02 m
Minimum stacks per trace	8	8 stack
Sampling interval	0.9225 ns	0.2344 ns
System run mode	GPS positioning system with RTK (TRIMBLE GPS 5800 with TSC)	Odometer survey wheel

Unfortunately, an ERT survey was not available at the Riyadh site. Therefore, we used only GPR as an essential tool for subsurface investigation.

4. Geophysical data processing

The low conductivity of the subsurface limestone in both survey sites (Al-Mokattam & Riyadh sites) provides an exceptional GPR application for fractures and cavities detection. The GPR reflection data that are displayed directly after data acquisition in the field (Figures 6(a) and 7(a)) usually contain noise or interference, which make data interpretation difficult [44]. To clean up the noise and correct the horizontal and vertical scales of the raw field data, they must be subjected to processing prior to interpretation.

All GPR data acquired with the 400 MHz and the 100 MHz shielded antennas were processed in exactly the same way (Figures 6 and 7) for each individual profile using the REFLEXW program [45]. Table 2 indicates the sequence of processing parameters applied to GPR data acquired with the RAMAC-100 MHz and GSSI-400 MHz shielded antennas. The processing steps were as following: (1) “Subtract DC removal” and “time-shift correction” compensates for zero time drift due to temperature changes affecting the GPR electronics, within a predefined time range the mean is calculated for each trace and is subsequently subtracted from all samples of each trace. In this processing step, we also remove zero time shifts (t0) between individually recorded data sets by aligning all first breaks to one common level throughout the GPR profile and also between repeat profiles. (2) The “dewow” [46] step removes very low-frequency components of the data associated with either inductive phenomena or analog/digital conversion. (3) Automatic gain control (AGC): RAMAC profiles acquired with the 100 MHz antenna at the Riyadh site were subjected to AGC as one additional processing step. A fixed gain curve is normally applied in the field with the GSSI system but can be applied later during data processing for the RAMAC GPR system [47]. (4) Background removal (average trace removal): This processing step is used to remove or decrease the effect of horizontal reflectors by calculating an averaged trace from the selected time/distance range of the actual GPR profile and subtracting it from every trace in the section. (5) Bandpass frequency: The filter bands are specified by the setting of four frequency values represented by a cosine window. Amplitude spectrum for every antenna in the two different systems was plotted and used to determine the signal frequency bandwidth for the 400 MHz antenna and the 100 MHz antenna, respectively. Therefore, frequency bands chosen for the 400 MHz shielded antenna were bands of about 230–295 to 550–615 MHz, while frequency filter bands of about 60–110 to 320–3370 MHz, were chosen for the 100 MHz antenna. The frequency spectrum below the low and above the high cut-off frequencies was set to be zero. (6) Running trace averaging was used to enhance the data in the horizontal direction and make the features more visible by taking the mean of a predefined number of traces in a window and subtract it from each individual trace in sequence for the entire GPR section [48]. (7) Velocity estimation for time-depth calculations of the GPR data was based on diffraction hyperbola analysis with REFLEXW software [45]. In addition, RAMAC GPR system with 100 MHz unshielded antenna was used in this part of the study (Riyadh site) to acquire common mid-point (CMP) radar data (Figure 8). The transmitting and the receiving antenna for the 100 MHz unshielded antenna are separated and the CMP measurement can be carried out. Figure 8(a) shows the data of CMP profile along the survey site, where Figure 8(b) represents the velocity spectrum analysis derived from the CMP data. We can find that the root mean square (RMS) velocity is about 0.1 m/ns. Figure 9 shows two different examples of 2-D GPR processed profiles at the Al-Mokattam site. Interpretation indicates vertical fractures and irregular displacement features appearing as strong reflection patterns in the 2-D GPR profiles. Figure 10 indicates three different 2-D GPR profile examples at Riyadh site. The relatively large hyperbolic reflection signal was probably caused by relatively large cavity structure in the subsurface limestone. (8) The processed GPR profiles (Figures 9 and 10) were used to generate three-dimensional horizontal time-slices by spatially averaging the squared wave amplitudes of the recorded radar reflections over the time window (Figures 11 and 12). The interpolation process creates interpolated time-slices, which are normalized to 8 bit following the colour changes between different levels and not actual reflection values [49]. The number of generated time-slices depends on the selected time window length, the slice thickness and the overlay between slices. Horizontal time-slice thickness is usually set to be at least one or two wavelengths of the radar pulse transmitted into the ground [49]. Moreover, the thickness of the generated time-slice should not be thicker than the size of the smallest expected subsurface target (Figures 11 and 12), and in this case the selected window was about 25 ns. GPR-Slices v7.0 [49] processing software was used to create the 3-D visualization with isosurface rendering, to improve interpretation of the GPR results and produce more readable images that can accurately define the dimensions of the target and its continuity (Figure 13). Moreover, static topographic correction was not applied during data processing at either surveyed sites (Al-Mokattam & Riyadh areas), as the ground surface was smooth and flat and no elevation slope could be observed.

Figure 6. An example of GPR field data acquired with 400 MHz shielded antenna: (a) raw data, (b) data after dewowing, (c) data after signal processing, while (d) is the same profile with time cut after 80 ns to make fracture anomalies in subsurface more visible.

Figure 7. An example of GPR field data acquired with 100 MHz shielded antenna (RAMAC system): (a) raw data, (b) data after dewowing and AGC, (c) data after signal processing with background removal and bandpass filtering, (d) is the same profile with running average trace to smooth things up and make cavities features in subsurface more visible.

Figure 8. Common mid-point radar data acquired along Riyadh survey site using 100 MHz unshielded antenna: (a) CMP data collected and (b) velocity spectrum analysis obtained from (a) indicates that subsurface root mean square velocity is about 0.1 m/ns.

Figure 9. Two different examples of 2-D GPR processed profiles at Al-Mokattam site. The basic signal processing includes zero time-shift; background removal, bandpass filtering and running average trace. Interpretation indicates vertical fractures and irregular displacement features appearing as strong reflection patterns in 2-D GPR section.

Figure 10. Three different examples of 2-D GPR processed profiles at Riyadh site. The basic signal processing includes zero time-shift; AGC, background removal, bandpass filtering and running average trace. Strong hyperbolic reflection signal probably caused by relatively large cavities structure running in subsurface rock.

Figure 11. A composite set of time-slice images from surface to about 108 ns (corresponding to a depth of about 4.6 m). The strong reflection anomaly is related to subsurface fracture in limestone (Al-Mokattam site), where complexity and continuity of the fracture is not clearly visible in time-slice images.

Figure 12. A composite set of time-slice images from surface to about 310 ns (corresponding to a depth of about 15 m). The strong reflection anomaly is related to subsurface cavities in limestone (Riyadh site).

Figure 13. 3-D visualization by means of isosurface rendering technique of the trace analysis. (a) Al-Mokattam isosurface map illustrates the lateral continuity and intersection pattern of the steep fractures, where complex internal and geometry of fracture further enhanced in isosurface visualization. (b) Riyadh isosurface map, where expected subsurface cavities in limestone system cloud could be easily detected. The label (C1–C7) is given for the expected cavities location as the same in Table 3.

Table 2. Processing sequence applied to GPR data acquired with the RAMAC-100 MHz and the GSSI-400 MHz shielded antennas.

	Processing step	Description	Parameters used for 100 MHz antenna	Parameters used for 400 MHz antenna
1	Time-shift correction	Zero time correction is adjusted by delaying the first signal break by the time it takes the direct airwave to travel from transmitter to receiver.	−18 ns	−4 ns
2	The “dewow” step	Removes the very low-frequency components related to system-dependent noise.	Running mean value on the whole trace (0 ns to ∼300 ns)	Running mean value on the whole trace (0 ns to 120 ns)
	Automatic gain control (AGC):	AGC	40	–
3	Background removal	Remove the effect of horizontal reflectors and coherent constant noise	From 0 ns to ∼300 ns	From 0 ns to ∼120 ns
4	Bandpass filtering	Remove low and high frequency components	*LCO/LP/ UP/UCO: 60/110/320/370 respectively	*LCO/LP/ UP/UCO: 230/295/550/615 respectively
5	Running trace averaging	Used to smooth up the data and make features more visible.	7 traces	4 traces
6	Subsurface velocity calculation	Time-depth conversion, by using hyperbola fitting analysis from subsurface cavity target.	0.1 m/ns	0.1 m/ns
	Export data to GPR-SLICE		Exported as DZT format	Exported as DZT format
	Producing data in 3-D view	This processing step was used to show surfaces of equal amplitude in 3-D volume.	Generation of 3-D data cubes include the time-slice images and isosurface maps	Generation of 3-D data cubes include the time-slice images and isosurface maps
	*LCO: lower cut-off, LP: lower plateau, UP: upper plateau, UCO: upper cut-off

A geoelectric field dataset, in the form of 2-D electrical resistivity imaging, was processed using RES2DINV computer software [50] (Figure 14). The software generates a calculated model of the pseudo-section from the inverted model; the root mean square (RMS) error between the calculated and the measured pseudo-section is also computed during inversion sequences. A least-squares algorithm was used to reduce the RMS error between the measured and the calculated apparent resistivity in an iterative mode. Two inversion methods are available in the software package: the Gauss–Newton smoothness constrained least-squares [51] and the Gauss–Newton robust model constrained [52]; as a general rule, features depicted by both methods can be considered real [53]. The Gauss–Newton smoothness constrained least-squares method was selected to reduce the RMS error between the measured and the inverted apparent resistivity [53]. On the other hand, the non-linear least square optimization technique was used for the inversion routine [50]. The number of iterations for data inversion was selected to be 5 iterations; this iterative process aims to reduce and minimize the difference between the measured ERT pseudo-section and the inverted pseudo-section based on the starting model.

Figure 14. 2-D electrical resistivity section acquired with Wenner-Schlumberger configuration (Cairo site). A set of 80 electrodes with 1.5 m electrode spacing generate a 120 m profile length with about 23 m depth of complete data coverage. (a) and (b) The measured and inverted model of the electric resistivity profile, (c) its interpretation are based on resistivity values and observed fracture in the adjacent outcrop.

Subsurface interpretations for the inverted 2-D model of the electric resistivity profile are based on resistivity values and excellent limestone wall exposure, which facilitate the correlation between surface outcrops and ERT data as shown in Figure 2.

5. Results and interpretation

Instead of generating a large set of slices, an overlay between slices was applied to connect subsurface reflections that might have been recorded at slightly different depths as well as to improve image transitions between depths. The number of time-slices generated from GPR-Slice (Figures 11 and 12) depends mainly on the length of the time window selected, the slice thickness and the overlay between slices. The thickness of horizontal slices is often set to at least one or two wavelengths of the radar pulse that is sent into the ground [49]. This processing step figured out a set of 30 horizontal slices at the Al-Mokattam site, using a 6 ns thick interval over a 120 ns time window (Figure 11) and produced 26 horizontal slices for the Riyadh site using a 25 ns thick interval over a total two-way travel time of about 300 ns (Figure 12).

Rather than following the conventional procedures of interpreting 3-D GPR subsurface reflection anomalies from multiple subsurface horizontal time-slices, we decided to detect the difference in object reflections on the basis of their 3-D isosurface rendering. Illumination is frequently added to the isosurface rendering presentation to highlight the features contained in the data volume (Figure 13). This three-dimensional display of GPR data simulates the transmission and absorption of light through points inside the volume. Light rays are cast through the volume, where equal amplitudes within the volume may reflect or absorb light. The colour of an individual pixel on the screen is computed by combining contributions from each particle intersecting the ray. Colour-coding of the resulting structures significantly enhances the recognition of the relative position of the fractures and cavities and permits the viewing of homogeneity and heterogeneity within the objects with suitable adjustment of the opacity. As a final visualization step, the isosurface rendering technique was applied to extract a volumetric body that represents the shape of the target [47]. Therefore, interpretation of the GPR data was based on the concept of these generated isosurface rendering volume.

Figures 9 and 10 demonstrate typical 2-D GPR vertical cross sections extracted from the two survey sites after data processing. Fracture in the subsurface is probably correlated with the high-amplitude reflections and vertical displacement in subsurface reflectors (Figure 9), while cavities in limestone indicate relatively large hyperbolic features in the 2-D profiles (Figure 10).

5.1. Joint interpretation for GPR and ERT results at Al-Mokattam site

The characterization of fractures with GPR in subsurface rocks depends on the type and size of the rock and the subsurface soil moisture content of the investigated materials. If the subsurface fractures are quite large and filled with materials with high electric conductivity and/or high electric permittivity (clay- or shale-filled fractures or high-moisture content filling), then a high amount of reflected radar energy will be scattered [54]. Most subsurface fractures are too thin to be detected by GPR 2-D sections alone; therefore, we applied an isosurface rendering technique to visualize the continuity of the fracture and extract the volumetric distribution of the high-moisture content of the fracture filling material. The fractures could be observed in 2-D GPR sections by their semi-vertical displacement with openings reaching more than 0.5 m causing numerous relatively small diffraction hyperbolas in the area where they occurred (Figure 9). Using these visualization approaches, a subsurface fracture network was created in the 3-D volume view (Figure 13(a)). Some of these fractures are confirmed by the surface exposures (Figure 2). The 3-D isosurface visualization technique of volume rendering (Figure 13(a)) shows only the strongest amplitude reflection by setting the weak signals transparent. By using this visualization approach, a subsurface fracture network can be interpreted and followed in the 3-D GPR cubic volume. The 3-D GPR results in the form of isosurface rendering volume image approximately located the 3-D subsurface filling material of the fractures in limestone, helped to easily trace the complexity of the subsurface fractures, and simplified the location for a suitable future engineering solution. These sediment-filled fractures and displacements indicate high-amplitude features in the GPR isosurface rendering volume.

ERT was used at this site to investigate and reduce the ambiguity in GPR data interpretation. The ERT inverted image (Figure 14) shows relatively low-resistivity values (from 36 Ωm to about 68 Ωm) interpreted as near-surface fractured limestone with relatively high-moisture content and/or shale and clay intercalation; various lateral changes in resistivity values within this layer could be related to water seepage from the nearby urban areas and/or heterogeneous solid materials filling the fractures in the limestone. At the lower part of the ERT section (Figure 14), we can observe a highly resistive layer (between 120 Ωm to values greater than 449 Ωm) that may correspond to dolomitic limestone. The depth to the top surface of this layer is interpreted to vary between 7 m in some location to more than 10 m in other locations. This dolomitic limestone exhibits in some parts of the section a relatively low resistivity suggesting shale and clay intercalation in limestone related to the Maadi formation [30]. In the ERT inverted model (Figure 14(c)), black dashed lines are assigned to the near-surface observed fractures running in limestone, while the white dashed lines at the lower part of the section indicate expected fractures in the dolomitic limestone layers. The ERT and 3-D GPR results are shown in Figure 15; the overlapping ERT section and isosurface rendering volume in the first 24 m highlights the coherency between these two techniques used to image subsurface fracture.

Figure 15. Overlay view between 2D electrical resistivity section and the 3-D isosurface view for GPR data demonstrate the efficiency of integration between the two techniques to image subsurface fracture at Al-Mokattam site.

The low-resistivity zones between 1.5 and 4, 9 m and 12 m, and at about 24 m (Figure 15) agree with the high-amplitude reflections in the GPR isosurface image at the same distance, implying a fractured limestone with high-moisture content filling visible in some parts of the survey site. This observation suggests that fractured and cracked rocks may be inter-connected to make the subsurface rock matrix more permeable; therefore, we put forth a plausible water infiltration model under a severely cracked and fractured subsurface rock condition [55–57], which may lead relatively low-resistivity values in the study area.

5.2. 3-D GPR results at Riyadh site

In this dataset Figure 12, the 3-D images produced in the form of slice maps and isosurface rendering (Figure 13(b)) have been generated using GPR-Slice software [49]; analysis of the multiple-reflection profiles indicates that the possible cavities of interest may be located in the subsurface limestone, and recorded in depth from 1.4 m to about 7 m in the surveyed region (Figure 3). Therefore, suitable threshold values were chosen to improve reflection from the cavities and to remove the common GPR reflections produced by the surrounding homogenous limestone (reflections from host soil limestone are not included in the final isosurface images as they show only a relatively low amplitude reflection compared with reflection from the cavities) (Figure 16). Integrated analysis of the individual GPR reflection profiles enables us to understand how the final 3-D images appear. Essential questions about GPR and what is being imaged in the subsurface should be raised during data analysis and processing. The most basic question is: What has produced the reflection? This can often only be determined from viewing reflections in profiles and understanding what produces radar reflections [58]. Subsurface limestone at the Riyadh site shows a relatively weak reflection to radar wave, as it is composed of very homogenous white chalky aphanitic limestone of the Sulaiy formation [59] (Figure 3). Therefore, any subsurface interfaces between the cavities and surrounding limestone materials will reflect the radar waves transmitted from the surface and appear in the GPR sections as diffraction hyperbolas (Figure 10). By analyzing the individual GPR sections, the reflective nature of the buried cavities of interest was distinguished and located in the standard amplitude slice map, where it becomes easy to follow and observe in the isosurface rendering image. At this site, isosurface rendering technique was able to image the detailed shape of subsurface cavities (Figure 16); the shape of the cavities indicates a circular structure. Furthermore, combining conventional time-slice images with 3-D isosurface rendering maps (Figure 13(b)) shows us a relatively high-amplitude reflection, in which interpretation may be due to expected high-moisture content in the subsurface limestone near and around the expected position of cavities. This high-moisture content in the subsurface limestone is probably related to water seepage from nearby site urban areas due to bad sanitary systems. Groundwater seepage is unpredictably moving under many areas in the city of Riyadh, causing cavity formation in the carbonate soil, infrastructure damage, and/or landslides. Therefore, detection of these subsurface cavities is essential to prevent future building damage and loss of human lives.

Figure 16. An isosurface rendering of the 70% strongest reflecting surface in the 3-D data cube, allowing detailed interpretation for expected cavities location at Riyadh site.

In this study area, drilled construction trenches confirmed the existence of subsurface cavities at various depths (Figure 3), where complicated networks of subsurface cavities have developed as a result of limestone dissolution (Figure 16). We have summarized in Table 3, the location and depth for the expected subsurface cavities anomalies at the Riyadh site. The large-amplitude reflection in the 3-D section (Figure 13(b)) around the cavities is probably correlated to limestone with high-moisture content due to water seepage from sanitary systems in nearby urban areas. The karst geomorphology is a distinctive terrain that occur in soluble rocks such as carbonates and evaporates limestone whose landforms are affected by underground seepage and drainage.

Table 3. Summary of expected location and depth for the subsurface cavities and anomalies at the Riyadh site, KSA.

Anomaly #	X-coordinates (latitude)	Y-coordinates (longitude)	Time (ns)	Depth (m)
Anomaly #01	2,737,528.500	666,051.310	42.827	2.141 m
Anomaly #02	2,737,525.750	666,025.440	72.270	3.613 m
Anomaly #03	2,737,549.250	666,017.750	96.360	4.820 m
Anomaly #04	2,737,537.750	665,983.250	42.190	2.109 m
Anomaly #05	2,737,460.062	666,053.690	199.412	9.971 m
Anomaly #06	2,737,461.340	666,027.310	52.195	2.610 m
Anomaly #07	2,737,560.750	666,041.250	49.023	2.450 m

6. Conclusions

In spite of the large fractured limestone wall exposure at our survey site, the subsurface fracture networks and especially the 3-D fracture continuity, in many cases are difficult to evaluate only from just the surface outcrop. Important information about the internal structures of the complex fractures and its continuity in the subsurface limestone rock recorded in 2-D profiles is not clearly visible in the conventional time-slices. Therefore, combining 3-D GPR in the form of isosurface rendering and 2-D ERT over such a fractured limestone is a major contribution toward tracking the fracture running in the Al-Mokattam limestone. This study shows that it is possible to map the soil filling the fracture, to indicate the lateral discontinuities and to identify the expected locations of future zones of collapse. Consequently, GPR in the form of isosurface rendering offers an unparalleled opportunity to determine sub-meter scale detail of fracture distributions, and possible locations of subsidence when it is integrated with the 2-D ERT, providing a well-constrained method for building fracture filling maps.

The study also shows that cavities developed along the Riyadh site are probably due to or associated with vertical water seepage from nearby urban areas. This supports the hypothesis that cavities develop in areas where groundwater is accumulated or temporarily retained. Moreover, in the Riyadh area as construction projects increase, this case study presents an appropriate and time-reducing way to plan the design and pre-evaluation of subsurface geo-engineering imaging for future projects running over formations characterized by cavities.

The proposed 3-D GPR visualization technique has shown its potential as a modern tool for GPR data interpretation. These techniques provide detailed information from different depths and facilitate the volumetric assessment of the underground fractures and cavities for engineering purposes.

Compared with the hazard and potential costs of damages caused by subsidence or rockslides in urban areas characterizing by unknown cavities or fractures, the cost and speed of execution of the discussed approaches (3-D GPR with advanced visualization techniques and ERT) for the geophysical prospection is definitely worth the effort. Moreover, we have shown that advanced 3D visualization by means of isosurface rendering procedures goes, beyond conventional time-slice imaging presentations, allowing for a comprehensive and more detailed interpretation of complex targets within the 3D GPR data volume. This 3-D visualization technique produced understandable images for fractured rock and cavity structures in great detail and provides engineers with information to guide future solution engineering developments technologies.

Acknowledgements

The authors would like to thank colleagues at NRIAG, Egypt for their help and support during field work. Moreover, we want to express our gratitude to Dr. Dean Goodman for helping to process the GPR data using GPR-Slice software.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The authors would like to thank deanship of scientific research, Majmaah University, Kingdom of Saudi Arabia, for funding this work [project number 38-66].