High-resolution seismic reflection surveying to delineate shallow subsurface geological structures in the karst area of Shenzhen, China

Near-surface seismic reflection surveys can produce high-resolution geological structural images for engineering and hydrological investigations. To help delineate shallow subsurface complex geological structures in a karst area and to better understand limestone cave formation, a high-resolution 2D seismic reflection profile was acquired and processed in the urban area of Shenzhen, China. The stacked images detail subsurface structures down to depths of 80–90 m, including a concave shaped reflection that curves upwards, two thrusts and one normal fault, as well as a hard rock basement reflection at the southern end of the profile which could not be mapped by borehole investigations due to the limited drilling depth. Our interpretations correlate well with borehole data and were supported by synthetic seismic modelling. Limestone caves are mainly found along faults and near rivers in this area. The results provide new insight on the formation mechanism and distribution of the karst caves. Movement along faults and erosion generated fractures and fissures in the limestone that provide channels for rainwater and groundwater to circulate. These waters, rich in carbonic acid, dissolve minerals in the limestone, resulting in the formation of underground cavities.

1 Introduction

Karst landscapes, characterized by their unique geological formations shaped by dissolution of soluble rocks, such as limestone and dolomite, are renowned for their beauty and ecological significance. However, beneath the surface hidden dangers are present. Karst hazards pose challenges for various industries, particularly those dependent on stable ground conditions, such as construction, agriculture, and infrastructure development (Chen et al., 2019; Medici et al., 2024). Losses due to the fast-acting nature of karst encounters can be direct (e.g. human casualties and damage to property), indirect (e.g. interruption to businesses, transport infrastructure and communication networks), or intangible, especially if they occur in areas of high population density (Galve et al., 2012; Bobrowsky, 2013; Intrieri et al., 2015; Sevil et al., 2017; Pazzi et al., 2018).

Figure 1 Regional (a) and local (b) geological map of the Pingshan district, Shenzhen (after Mao et al., 2014 and Shenzhen Geology Writing Group, 2009). Blue rectangle (Fig. 3a) shows study area location.

In southeast China, buried karst, with high fissure water content, high permeability, and variable shapes, is widely distributed (Cui et al., 2015). Shenzhen, a world-scale metropolis, located in the southern part of the Guangdong province (Fig. 1a), has significant areas of Carboniferous rock distributed in the Longgang, Pingshan, and Dapeng districts where karst features are a hazard. During metro, railway, and building construction, karst and karst caves at depths from 2 to 50 m are found in dolomite, dolomitic limestone, dolomitic marble, marble, crystalline limestone, and breccia limestone (Li et al., 2011; He et al., 2020). Its presence challenges the construction of tunnels and use of shield tunnelling machines. For example, the disturbance produced by a shield tunnelling machine may induce the ceiling of a karst cave to collapse, with water present in the cave damaging the machine and constructed tunnel. Furthermore, diaphragm wall collapse, water or mud ingression, ground collapse, and long-term instability are often encountered in karst regions (Cui et al., 2015). Currently, many direct and indirect techniques to detect buried karst and karst caves have been proposed (Lolcama et al., 2002; Hoover, 2003; Kaufmann and Romanov, 2009; Kaufmann, 2014; Samyn et al., 2014; Putiška et al., 2014; Kaufmann and Romanov, 2016; Sevil et al., 2017; Pazzi et al., 2018; Hussain et al., 2020; Wang et al., 2020; Muzirafuti et al., 2020; Stan-Kłeczek et al., 2022; Yordkayhun et al., 2022; Liu et al., 2023). Among the direct methods, drilling and electric cone penetration tests are the most common and useful. On the other hand, indirect techniques can be employed to delineate subsurface karst size and distribution and extrapolate borehole data to a wider area. Geophysical techniques based on a physical contrast between a cave and the surrounding rocks provide an economical and non-invasive alternative, or complement, for mapping geological structures, and are often used in attempts to detect the presence of karst caves and voids below the surface. Methods include seismic reflection/refraction, multichannel analysis of surface waves (MASW), the H $/$ V spectral ratio method, electrical resistivity tomography (ERT), induced polarization (IP), self-potential (SP), ground penetrating radar (GPR), and micro-gravity. In recent years, with a focus on improving resolution, efficiency, and cost-effectiveness (Juhlin et al., 2000; Martinez and Mendoza, 2011; Brodic et al., 2015; Yilmaz, 2021; Sun et al., 2022; Pertuz and Malehmir, 2023), the high-resolution seismic reflection method for shallow exploration and imaging of local subsurface heterogeneities has emerged as a powerful technique to identify and map near surface geological structures with good precision and depth penetration.

In this study, we delineate the near surface geological structures in a karst area with high-density and -resolution seismic reflection data to carry out the research specific objectives as follows: (i) testing of the reflection seismic method over karst terrain; (ii) comparison of reflection seismic results with geotechnical drilling; (iii) integration of physical properties of samples with reflection seismic modelling; and (iv) delineating shallow subsurface complex geological structures in a karst area and understanding limestone cave formation better. Our seismic results are correlated with borehole data and synthetic modelling simulated by the 2D elastic wave equation. The seismic reflection images provide new insight in understanding the formation mechanism and distribution of karst and karst caves.

2 Geological setting and physical properties

Geologically, Shenzhen is located in the South China Block southwest of the Zhenghe-Dapu fault (Fig. 1a). It consists of a multi-stage, complex formation of folds, faults, and intensive metamorphism which took place in the pre-Caledonian, Caledonian, and Hercynian to Indosinian orogenies (Cui et al., 2015), with the area having undergone a complex tectonic evolution from assembly to break-up to re-assembly (Shu et al., 2014). The early activity along the Zhenghe-Dapu fault zone was formed through ductile shear, and the middle and late stages are composed of thrusting with a mid-development detachment. The conspicuous Lianhuashan fault zone in Shenzhen is interpreted as the southern extension of the lithosphere-scale Zhenghe-Dapu fault zone, considered as the boundary between the late Mesozoic Coastal terrane and the early Paleozoic Wuyi-Yunkai orogen (Li et al., 2020). The sedimentary units affected by the Lianhuashan fault zone include Paleozoic and Early to Middle Jurassic rocks. Due to significant magmatic activity, granite and igneous rocks that formed during the Yanshanian period are the dominant rock type. The Yanshanian orogeny involved tectonic movement that mainly included block orogeny and the development of deep faults and wide folds, along with significant intrusions of granitic magma, massive acidic volcanic eruptions, and magmatism-related mineralization (Li et al., 2020).

On a local scale (Fig. 1b), the Paleozoic sequence comprises Silurian-Carboniferous sandstone, shale, slate, limestone, siltstone, gritstone, and coal seams (Shenzhen Geology Writing Group, 2009) that are distributed along the NE-striking and NW-dipping Paotaishan-Hengtougang transpressional fault (PHF) and the Shijingling-Huangzhukeng transpressional fault (SHF). The Lower Carboniferous comprises two units, the Shidengzi formation (C₁s¹) and the Ceshui formation (C₁c¹). C₁s¹ is buried under the Ceshui formation and not exposed in the area. There are shallow marine facies carbonate rocks, consisting of grey and dark grey crystalline limestone which mostly have been metamorphosed into off-white marble and dolomitic marble, with a thickness of more than 340 m. C₁c¹ contains coastal facies carbonaceous clastic rocks, composed of argillaceous shale, siltstone, fine sandstone, etc., with a thickness of more than 1200 m. Holocene alluviums and flood plain deposits affected by the PHF and SHF consist of gravel, silty sand, sandy clay, silty clay, and other unconsolidated deposits. Medium-grained granites are present northwest of the study area with U–Pb zircon crystallization ages of about 108 Ma (Shenzhen Geology Writing Group, 2009). The seismic survey line is located on Holocene flood plain deposits, as shown by the blue rectangular area in Fig. 1b.

Figure 2 Borehole lithology and representative core photos. (a) Borehole section consists of 19 boreholes down to a maximum depth of 70 m, located about 20 m away from the seismic survey line. (b) Core box images from BH62 borehole, 0–6.9 m, fill soil; 6.9–27.3 m, silty clay; 27.3–70.3 m, weathered shale, some weathered into clay and some are still hard rocks.

The geological section is based upon 19 engineering geology boreholes that provide more detailed information on lithology down to 70 m (Fig. 2). There are six main geological units from the top to the bottom of the section, including fill soil, sand, silty clay, weathered sandstone, and shale, and limestone. Limestone caves were found in boreholes BH52 and in BH56 to BH59 in the vicinity of the interface between weathered shale and limestone. The fill soil has a relatively constant thickness of ca. 10 m in all boreholes. The thickness of the sand layer increases from north to south, however, weathered shale shows the opposite trend. Weathered sandstone is only present from BH46 to BH51. Limestone is present from BH46 to BH59 and absent from BH60 to BH64. Thickness variations of the different sedimentary deposits in the horizontal direction suggest that multi-stage tectonic events occurred in this area. The top of the weathered shale has a concave upward shape between BH53 to BH57. Based on borehole sampling, the groundwater table is present at the top of the sand layer. Due to the different compositions of the inner shale, differential weathering is suggested to have occurred at the southern end as indicated in BH60 and shown in Fig. 2b; it is a potential source for seismic reflections.

Compressional and shear velocity data were acquired by velocity logging and density data was obtained by the volumetric cylinder method in the lab using samples from the vicinity of study area. These measurements provide information on the physical properties of the subsurface soils and rocks (Table 1) and can help identify where reflection impedance changes may be located. Fill soil with a density of 1550 kg m⁻³ has the lowest density and limestone with 2650 kg m⁻³ has the highest value; the density of sand with gravel is somewhat higher than silty clay and other weathered rocks. Compressional and shear velocities have similar variational trends, however, compressional velocity is 3–4 times faster than shear velocity. Limestone compressional velocity is up to 3125 m s⁻¹; sand with gravel and silty clay have the lowest compressional velocity. Unconsolidated soils including fill soil, sand with gravel, and silty clay have the highest ratio of compressional to shear velocity; the ratio of compressional and shear velocity of weathered rocks is higher than limestone. A four-layer seismic reflection model is suggested in this area based on the velocity and density data. The first layer consists of fill soil, the second is made up of gravel and silty clay, the third is completely and strongly weathered rocks (including sandstone and shale), and the fourth is limestone.

Table 1 Compressional velocity (V_p), shear velocity (V_s), and density of soil and rock collected from the adjacent area. V_p and V_s are from logging velocity and density from geotechnical testing of borehole samples.

3 Data acquisition and geometry

High-resolution seismic data were acquired in November 2022 in the Pingshan district of Shenzhen, southeast China (Fig. 1). The survey line is located along Longping Road and lies between light rail transit line 1 and the Ciao river (Fig. 3a, c). The surface topography from south to north varies smoothly with an elevation difference of about 2 m. The geological section from the 19 boreholes allows a comparison with the seismic results, however, they were drilled at an elevation ca. 5 m higher than the seismic survey line (Fig. 3c) and 20 m offset from it. A 5 kg sledgehammer with a 4 cm thick metal plate was used as the seismic source along with 5 Hz SmartSolo 3C nodal units for recording. The seismic data were acquired using five segments. Every segment consisted of 148 units with a fixed geometry of 1 m receiver spacing and 2 m source spacing with 74 units overlapping (Fig. 3b). The total length of the survey line is 417 m. A sampling rate of 1 ms was used and 1000 ms of data were retrieved from the nodes for each source point. Table 2 shows the seismic acquisition parameters. Three raw shot gathers recorded at different locations along the survey line show some of the characteristics of the seismic wave field (Fig. 4). Direct waves, surface waves, reflections, and air waves are present in all gathers. Reflections dominate in the time window of 50–150 ms; shallower reflections are masked by surface waves, air waves, and direct waves at the near offsets.

Figure 3 (a) Locations of receivers (cyan triangles), shots (red stars), and boreholes (green dots). (b) Geometry of the seismic acquisition profile, receiver spacing 1 m, shot spacing 2 m, separated into five segments, each segment consisted of 148 fixed units with 74 units overlapping. (c) Field work photo of the seismic survey line (cyan line) and geological section consisting of 19 boreholes (red line); the geological section is 5 m higher than the seismic survey line (the aerial image from © Google Maps).

Table 2 Acquisition parameters for the high-density and -resolution seismic data in the municipality of Shenzhen, southeast China.

Figure 4 Three typical shot gathers from different locations along the seismic survey line: (a) located in the south (shot No. 968), centre (shot No. 1214), and north (shot No. 1374). Reflections from the time window of 50–150 ms are notable, shallower reflections are masked by surface waves, air waves, and direct waves at the near offsets. F, S, R, and A represent first breaks, surface waves, reflections, and air waves, respectively.

4 Data processing

Data processing followed a standard workflow to improve the signal-to-noise ratio and resolution after stacking vertically repeated shot records (twice at each shot location). Tomographic refraction statics were applied to account for travel time variations in the very near surface, while band-pass filtering, spectral equalization, and linear noise suppression were applied to remove different types of noise, including air waves, surface waves, and coherent noise with low frequency. An iterative velocity analysis was performed to obtain the best velocity model and continuity of reflections. The continuity was improved further through surface-consistent residual static corrections. After stacking, band-pass filtering and f−x domain deconvolution were used to further reduce random noise and enhance the coherency of reflections. A finite difference time migration routine was used to move dipping reflectors to their true subsurface positions and collapse diffractions after stacking. Finally, time to depth conversion with the smoothed NMO velocity field was performed to obtain an approximate depth model and to help in interpreting geological structures. Processing steps and parameters are shown in Table 3, while two important steps are discussed in detail below.

Table 3 Processing workflow and parameters for the high-density and -resolution seismic data in the municipality of Shenzhen, southeast China.

4.1 Noise attenuation

Raw shot gather recordings are dominated by low frequency and strong amplitude noise, such as surface waves and coherent noise. A typical power spectrum shows that the data contain frequencies in the range 5–440 Hz (Fig. 5a). To suppress the low frequencies and strong amplitude noise, a band-pass filter with corner frequencies of 50–80–180–220 Hz was applied to the raw data after carrying out different frequency ranges tests (Fig. 5b). Direct waves, refractions, reflections, air waves, and surface waves are recognized with apparent velocities of approximately 1000, 1800, 2200, 340, and 220 m s⁻¹, separately. After band-pass filtering, a spectrum equalization process was used to further reduce noise and enhance the weaker amplitude signals. Compared with Fig. 5b, two sets of reflections are notable in Fig. 5c and marked with arrows. Finally, the air waves, surface waves, and other linear noise were attenuated by median filters at velocities of 340, 250 and 600 m s⁻¹, respectively. Reflections are now clearly observed in Fig. 5d.

Figure 5 Noise attenuation processing steps for the No. 964 shot gather. (a) Raw shot gather and its power spectrum; (b) after a band-pass filter (50–80–180–220 Hz), typical apparent velocities for direct waves, refractions, reflections, air waves, and surface waves are about 1000, 1800, 2200, 340, and 220 m s⁻¹, respectively; (c) after spectrum equalization (50–60–200–220 Hz); (d) after median filtering to remove noise, air wave 340 m s⁻¹, surface wave 250 m s⁻¹, and linear noise 600 m s⁻¹.

4.2 First break travel time tomography

First break picking and travel time tomography were performed to correct for the seismic wave travel time delays in the very near surface low velocity zone and to build a P-wave velocity model. Based on the picked first breaks, a two-layer model was chosen as the initial velocity model, with apparent velocities of about 700 and 1650 m s⁻¹ (Fig. 6a). The ray density plot shows an even distribution in most regions except at about 100–240 m along the profile and at elevations of 20–28 m (Fig. 6b). Root-mean-square errors are approximately 2.4 ms after the fifth iteration (Fig. 6c), indicating a stable result has been attained. The inverted velocity model suggests a two-layer structure over the depth of investigation, the top one with low P-wave velocity is consistent with the fill soil layer; the bottom one with a high P-wave velocity may be correlated to the sand-silty clay layer (Fig. 6d).

Figure 6 Travel time tomography for first breaks. (a) First break travel times for all shots, a two-layer model with apparent velocities of 681 and 1667 m s⁻¹ was taken as the initial velocity model. (b) Ray density model of 20–36 m above sea level. (c) Root mean square (RMS) errors for different iterations. (d) Inverted P-wave model, the black dashed lines represent no rays in the area and velocity is interpolated below it.

5 Synthetic modelling

2D synthetic seismic modelling was carried out using a finite-difference solution of the elastic wave equation (Juhlin, 1995). A four-layer model (Fig. 7a) with fill soil (Layer 1), sand and silty clay (Layer 2), weathered shale (Layer 3), and limestone (Layer 4) was used to simulate the subsurface geological structures. The model parameters for the compressional velocity, shear velocity, and density are based on the values in Table 1 and first break velocities of the compressional and shear waves using the vertical and radial components, respectively. The 417 m length model is the same as for the real data and the depth is to set to 300 m to mitigate reflections from the bottom of the model. The model was discretized to a 0.15 m × 0.15 m grid size. A Ricker wavelet with a central frequency of 70 Hz was used as a source. In order to be consistent with the field data, the geometry and recording parameters were identical with Fig. 3b and Table 2.

Figure 7 Synthetic model and results. (a) the model of P-wave, S-wave, and density based on borehole data and field raw data, Layer 1, Layer 2, Layer 3, and Layer 4 represent fill soil, sand and silty clay, weathered shale, and limestone, respectively, and can be correlated with borehole data. (b) Synthetic shot gather at shot No. 964, the same location as Fig. 5a–d. (c) Pre-stack time migration section.

Figure 7b shows one synthetic shot gather, which has the same location as the real gather in Fig. 5a–d. The R3 reflection is generated from the top of the limestone and matches well the real data gather. Pre-stack reverse time migration was applied using the CREWES MATLAB Toolbox (Fig. 7c). The reflections from the top of sand and silty clay (R1), weathered shale (R2), and limestone (R3) are clearly visible. The result also suggests the presence of a concave shaped structure consistent with the borehole section (Fig. 2a) and real data (Fig. 8). At a distance 0–100 m along the profile and a time of ca. 100 ms, a strong amplitude reflection indicates the limestone which is not mapped by the borehole data due to the limited borehole depth. This reflection is clear in the real data and provides added constraints on the geological structure of the area (Fig. 8).

Figure 8 Seismic reflection migrated image overlaid with the time converted tomography result; an apparent discontinuity of reflections and velocity at horizontal distance 220 m is observed. R1 and R2 represent reflections from the top and bottom of the sand and silty clay layer, respectively; R3 originates from the top of the limestone. At horizontal distance 220 m and time 75 ms, a concave shaped structure exists that correlates with the borehole section.

6 Results and discussion

6.1 Reflection characteristics from the migrated seismic section

Reflection seismic processing steps focused on removing random and coherent noise and imaging high-frequency reflections in order to resolve the near-surface geological structure. The seismic tomography image suggests that the uppermost reflection horizon on the real data can be interpreted as the fill soil – sand and silty clay interface (Fig. 8). Furthermore, the physical property sample analyses indicated that three seismic interfaces can be mapped in the area and we interpret our seismic section accordingly. In addition to the top of the sand and silty clay unit (R1) two other sets of reflections are imaged down to depths of 80–90 m. These reflections, labelled R2 and R3, represent reflections from the weathered shale and the bedrock interface, consistent with the borehole data and synthetic modelling. Between R2 and R3, some reflections are produced by clay from the weathered and unweathered shale, as shown in Fig. 2b. The R1 reflection is distorted at a distance of 220 m. When compared to the near surface velocity model from first break tomography (Fig. 8), there is also a variation in velocity at this location. This suggests that the thickness of fill soil varies due to palaeotopography or fault slipping. However, it is not possible at present to determine which factor is correct. We have interpreted three faults in the middle part of the seismic section, F1, F2 and F3 (Fig. 9). Seismic reflections show an apparent offset at the locations of these suggested faults; together they form the upward concave structure seen in the seismic image and in the seismic modelling. The locations also correlate with where karst caves are found in the borehole section. If correctly interpreted as faults, this indicates a connection between the karst caves and faulting in this area.

Figure 9 Geological interpretation of the depth converted seismic reflection section. Four-layer model as indicated by the seismic data with fill soil, sand and silty clay, weathered shale, and limestone. Two thrusts and one normal fault are shown and labelled by F1, F3, and F2, respectively.

Figure 10 Core box images from BH57 (a) and BH59 (b). Karst cave filled with silty clay was developed at depths of 42.8–48.8 m in the borehole BH57 and unfilled at depths of 52–55.9 m in the borehole BH59.

6.2 Formation mechanism and distribution area of karst caves

Karst distributed in the eastern part of Longgang, Pingshan, and Dapeng districts, Shenzhen, can be classified as covered karst with a buried depth of 4–40 m and buried karst with a buried depth of greater than 40 m (He et al., 2020). Buried karst is developed in dolomitic limestone in the research area. Karst caves, unfilled and filled with silty clay, were revealed from the drill cores of BH52, BH56, BH57, BH58, and BH59 (Fig. 10). The discovery rate of karst caves and the linear karst rate of boreholes were 26.31 % and 1.38 %, respectively. No sinkholes and dolines were developed at the surface. Karst caves and sinkholes are usually localized along distinct lines situated along the lithological and tectonic discontinuities (Putiška et al., 2014). A limestone micro-scale cavity readily develops along the fault planes and joints due to weathering and chemical corrosion (Song et al., 2012). Groundwater flows along discontinuities such as bedding planes, joints, and faults. Preferential structural deformation along faults or bedding can enhance dissolution during subsequent interaction with groundwater, so physical rather than chemical properties of some carbonates can control karstification (Lolcama et al., 2002). Previous studies suggest that at least four stages of regional scale tectonic and magmatic events have taken place in South China (Shu et al., 2014), forming sets of transpressional fault structures, N to NW dipping (Fig. 1a). As shown in Fig. 9, two thrusts labelled as F1 and F3 and one normal fault (F2) are suggested to be present, which may help in understanding the formation and distribution of karst caves in this area. Slipping along these faults and erosion generated fractures and fissures in the limestone that probably provided channels for rainwater and groundwater to enter the subsurface from the Ciao river (Fig. 3c). Carbonic acid in the water dissolved the limestone until it became saturated with carbonate minerals. As saturated groundwater flowed away from the area, unsaturated underground water flowed in and further dissolved the limestone. Finally, underground cavities were formed (Fig. 11a). This may give a reasonable explanation why limestone caves developed along the faults and near the rivers in the Pingshan district, Shenzhen (Fig. 11b), however, this needs to be proved by further work.

Figure 11 Schematic model of karst cave formation in the study area (a) and Karst distribution in the Panshan district, Shenzhen (b). Karst is not only correlated with carbonate rocks, but also faults and rivers (after Shenzhen Geology Writing Group, 2009).

7 Conclusions

A 417 m long high-resolution seismic reflection profile was acquired to image shallow subsurface geological structures in a karst area. Notably, four layers with fill soil, sand and silty clay, weathered shale, and the bedrock are imaged in high resolution on the stacked seismic section. Also, three small-scale faults are delineated which are closely related to karst caves. To validate the seismic image, synthetic modelling using the elastic wave equation was conducted to compare with the real data and help in the geological interpretation. Furthermore, a near surface velocity model from first break tomography, coincident with a borehole section 20 m away from the seismic line, correlates well with the seismic reflection image and confirms the reliability of it. Integrating the geology and geophysics results in the delineation of the fine subsurface geological structure and provides a better understanding of the formation process and spatial distribution of karst caves in the Pingshan district, Shenzhen. Multi-stage regional scale tectonic and magmatic events produced cracks, fractures, and fissures in the limestone that formed channels where rainwater and groundwater with carbonic acid could circulate and form caves and voids. This interpretation provides valuable insight into the development of karst caves and the research contributes to the mitigation of karst hazards.