3D Full-Waveform Tomography of Standard Penetration Test-Seismic Wavefields in Karst Florida Limestone

Abstract

The shallow subsurface across large parts of Florida consists of weathered karst limestone, a spatially variable porous stratigraphy with significant influence in the design and load capacity of deep foundations. For reliable load transfer of the super structure into the bearing soil/rock, detection of voids and weak zones is crucial. This study investigated the capabilities of a recently developed standard penetration test (SPT)-seismic testing method for characterizing a large volume of soil/rock properties with a single SPT. The method utilizes a three-dimensional full-waveform inversion (3D FWI) of wavefields induced by SPT blows at depths (in-depth source) to characterize the subsurface around the SPT borehole. A site in Florida that presents shallow, surface karst limestone, was the location of a field experiment that consisted of two 36- × 18-m areas, each with an SPT at the area center. Seismic wavefields induced by SPT blows at depth were recorded by 72 vertical geophones on the ground surface and analyzed by the 3D FWI. The results revealed that subsurface soil/rock properties were characterized in submeter pixels over a large 3D domain of 24 × 36 × 18 m (depth × length × width). Multiple voids at various depths from 5 to 17 m were successfully detected at the site and confirmed by SPT data. The results suggested that the SPT-seismic method is an efficient tool for site investigation, as a bridge pier or pile group could be designed with only one SPT.

Keywords

3D FWI SPT-seismic subsurface exploration 3D subsurface imaging void in karst

Deep foundations (pile groups) are often used to support heavy loads from superstructures such as bridges and buildings. As pile capacities are from skin and tip resistances, solid soil/rock around and below a pile group are required to carry the loads. A void near or below a pile group may cause its failure or even collapse. For foundation design, soil properties are often obtained from the standard penetration test (SPT) or cone penetration test, and rock properties are determined from SPTs, measuring while drilling or coring. However, these invasive tests can only provide material properties within a borehole; offline voids (i.e., away from the borehole) cannot be found. Multiple invasive tests may be needed for the design of large pile groups, particularly in highly variable sites. This study aimed to reduce the number of invasive tests needed for the design.

Surface-based geophysical methods such as seismic tomography, electrical resistivity (ER), and ground penetrating radar (GPR) are often used to increase the volume of test materials. Compared with the ER or GPR methods, seismic methods are more favorable because of their ability to provide highly characterized resolutions with depth, and their engineering properties (e.g., S-wave and P-wave velocities, Young and shear moduli) that are useful for geotechnical analyses. For example, with surface-based seismic testing, full-waveform inversion (FWI) methods ( 1 , 2 ) can provide two- and three-dimensional (2D/3D) S- and P-wave velocity profiles at meter pixels to a depth of 20 m. They have also been used for the detection of shallow voids at depths of less than 10 m ( 1 , 3 ). However, because of limited energy penetrating into rock masses from surface sources and high wave attenuation with depth, the void detectability of these FWI methods is limited to a depth of a few void diameters. A void embedded at a depth greater than 5 void diameters is typically undetectable by these surface-based seismic methods.

Borehole-based seismic methods can be used to increase the characterized resolution and depth of investigation. Several studies (4 –6, 7 –11) have reported simulation and inversion borehole data. Cross-hole seismic FWI methods ( 5 - 7 ) have shown excellent capabilities in characterizing materials between a pair of boreholes. However, they require several boreholes for 3D characterization.

In an effort to reduce the number of boreholes, Mirzanejad et al. presented a novel SPT-seismic method, which couples 3D FWI with the SPT for deep site characterization ( 12 ). Specifically, as the SPT hammer strikes the drill rod and soil sampler at various depths, seismic waveforms are generated and measured on the ground surface by a 2D grid of geophones. Unlike surface-based wavefields dominated by Rayleigh waves, SPT-seismic wavefields are dominated by body waves, emitting from point sources within the rock mass, passing through volumes of deeper materials, and arriving at the ground surface. Inversion of such seismic data provides a higher characterized resolution and greater accuracy at depth than those of surface-based methods ( 12 ). This case study further investigated the method’s capabilities in detecting multiple voids and characterizing soil/rock over a large volume with a single SPT.

Methodology

The 3D SPT-seismic FWI method has been presented in detail ( 12 ), but is briefly summarized in two steps as follows:

Step 1: simulate wave propagation using 3D elastic wave equations as,

ρ \dot{v_{i}} = σ_{i 1, 1} + σ_{i 2, 2} + σ_{i 3, 3}

(1)

ρ \dot{σ_{i j}} = λ (v_{1, 1} + v_{2, 2} + v_{3, 3}) + 2 μ v_{i, j} i f i = j

(2)

ρ \dot{σ_{i j}} = μ (v_{i, j} + v_{j, i}) i f i \neq j

(3)

where

i, j are of 1, 2, 3, representing x, y, and z directions, respectively;

$σ_{ij}$ is ij-th component of stress;

$v_{i}$ is i-th component of particle velocity;

ρ is density;

λ and µ are Lamé’s coefficients, computed from Vs, Vp, and density;

. (the over-dot) is time derivative; and

, (comma) is spatial derivative (e.g., $v_{i, j} = \partial v_{i} / \partial x_{j}$ ).

To solve these equations, the staggered-grid finite-difference method ( 13 ) is used for the grid discretization, together with the image technique ( 14 ) for the top free-surface boundary and the perfectly matched layer ( 15 ) for the other five boundaries. See Nguyen and Tran’s research for a detailed implementation of the forward simulation ( 16 ).

Step 2: update material properties (Vs and Vp) by minimizing the misfit function E(m) via Gauss-Newton optimization as,

E (m) = \frac{1}{2} [F (m) - d]^{t} [F (m) - d],

(4)

m^{n + 1} = m^{n} - α^{n} {[J^{t} J + λ_{1} P^{t} P + λ_{2} I^{t} I]}^{- 1} J^{t} [F (m) - d] .

(5)

In Equation 4, $d$ and $F (m)$ are field (observed) and simulated (estimated) waveform data, respectively. $F (m)$ is computed via the numerical solution of Equations 1 and 3 for model m, which contains Vs and Vp values of cells. The transcript, $t$ , is the matrix or transpose of the vector. In Equation 5, model m is updated iteratively until the convergence criterion is reached. The Jacobian matrix, J, is the derivative of $F (m)$ with respect to model m. Symbol I is the identity matrix, P is a matrix representing the 3D Laplacian operator, and $λ_{1}$ and $λ_{2}$ are constants. Step length $α^{n}$ is fixed at 0.5 during analyses. See the study by Tran et al. for details on the efficient computation of the Jacobian matrix and selection of optimization parameters ( 2 ). Note that only Vs and Vp are updated, whereas mass density is assumed and kept constant during inversion. Our efforts to invert the density were unsuccessful owing to its low sensitivity to wavefields. Thus, a homogenous density of 1,800 kg/m³ (typical soil/rock density) was used for field data analyses in this study.

Field Experiment

The field experiments were conducted at a site in Newberry, FL (Figure 1a). It is a dry retention pond, consisting of fine sand and silt layers underlain by weathered karst limestone. The top of the limestone is typically found at a depth between 2 and 10 m across the site ( 17 ), with the exception of a few locations where limestone can be seen on the surface. The site is prone to sinkhole activity, and several voids have been found at shallow depths of a few meters ( 1 , 3 ).

Figure 1.

Field experiment: (a) two test zones with an SPT at the center of each zone and (b) site photo with the SPT rig.

Data acquisition was performed for two test zones of 36 × 18 m, each with an SPT at the center as shown in Figure 1a. For each zone, 72 vertical geophones of 4.5-Hz resonance frequency were used on the ground surface, covering a 6 × 12 grid of 3-m spacing as shown in Figure 2. This test configuration was selected based on our parametric studies (not shown for brevity). A 2D grid of geophones at 3-m spacing (or less) is required to achieve acceptable accuracy and the desired resolutions at submeter pixels.

Figure 2.

Test configuration used for each test zone, with 72 geophones placed on the ground surface in a 6 × 12 grid of 3-m spacing and SPT at the center of the geophone array.

The SPT was used as the seismic source. For each source location at depth, the 63.5 kg (140 lb) SPT hammer was dropped 760 mm (30 in.), striking the top of the SPT rod (above the ground surface), and imparting energy that propagated to the SPT spoon (attached to the end of the SPT rod), which acted as the in-depth seismic source and induced seismic wavefields at various depths. At each hammer strike (blow), the seismic trigger, mounted on the SPT rod, activated the seismograph for data recording.

Test Zone 1

SPT 1 was conducted to an 18.9-m depth, with intervals of 0.45 m above the 6-m depth and 1.5 m below the 6-m depth. SPT-seismic data were recorded for 5 to 10 blows at each interval, for a total of 146 blows from the ground surface to the bottom of boring. From about a 13- to 17-m depth, the weight of hammer was recorded for the SPT blow counts (N = 0), suggesting a void that could be partially filled with raveled soils. Note that no seismic data were recorded at the void depth.

Shown in Figure 3 is the measured field data for the deepest source at a depth of 18.9-m. The data were filtered through a 10- to 40-Hz bandwidth (Figure 3a) with the corresponding frequency spectra (Figure 3b). The wavefield emitted from this deep source location arrived at all receivers with a consistent propagation pattern in both magnitude and phase. This suggests that the SPT hammer provided sufficient wave energy at the depth required for propagating over large distances (26 m to the farthest geophone) to the ground surface. A similar quality of measured data was obtained for most blows.

Figure 3.

Field experimental data for the deepest source at 18.9-m: (a) measured waveform data and (b) frequency spectra.

For inversion, the analyzed medium of 24 × 36 × 18 m (depth × length × width) was discretized into cells of 0.60 × 0.75 × 0.75 m (depth × length × width). The source locations (i.e., actual depths of the SPT spoon) were assigned to the nearest numerical nodes (0.6 m vertical spacing) for wave simulation. The 146 blows were assigned to 20 depths, and data from blows at the same depth were averaged to enhance signal quality.

Based on the spectral analysis of surface data at the site (Mirzanejad et al. [12]), a one-dimensional initial velocity model (Figure 4) was taken with Vs linearly increasing from 200 m/s for soil on the surface to 500 m/s for limestone at the model bottom (24-m depth). The initial Vp was assumed to be twice that of the initial Vs. The inversion started on the initial model (Figure 4) using waveform data was filtered through a 10- to 40-Hz bandwidth. The inversion stopped at a preset maximum number of iterations (100 iterations) and took about 26 h on a desktop computer (32 cores of 3.4 GHz, 512 GB of RAM). The high computing demand was the main limitation of this FWI approach.

Figure 4.

Field experiment: initial Vs used for the inversion of data from both test zones, and initial Vp as twice that of Vs (not shown).

Figure 5 presents the comparison of waveform data for two source locations at an 11-m depth (above the void) and 18.9-m depth (i.e., the deepest source). The estimated and observed (measured) data matched well at the end of the analysis. No cycle skipping was observed, confirming that inversion converged to the global minimum. Figure 6 shows the normalized error for all 100 iterations. The error reduced continuously from 1.0 at the beginning to 0.67 at the end of the analysis. The remaining error (waveform misfit) was a result of the field conditions (noise, damping, soil/rock properties) not being fully captured by the forward simulation and model update. Nevertheless, the final estimated data mimicked the field data (Figure 5).

Figure 5.

Test zone 1: waveform comparison between observed data and estimated data: (a) source at 11.0-m depth (right above the void) and (b) source at 18.9-m depth (deepest source). Poor data from two channels were removed from the analysis.

Figure 6.

Test Zone 1: normalized least-squares errors of all 100 iterations.

Shown in Figure 7 are the final inverted results. Vs profile (Figure 7a) contains two distinct layers: 1) a soft layer with Vs∼200 m/s from 0 to about a 7-m depth and 2) a stiff layer with Vs∼500 m/s from a 7- to 24-m depth. They agreed with the material samples from SPT with soil underlain by limestone. Importantly, a large low-velocity zone (Vs∼100 m/s) was found at the middle of the medium from about a 13- to 17-m depth. The Vp profile (Figure 7b) was similar to the Vs profile, including the existence of the low-velocity zone and soil/limestone layers.

Figure 7.

Test Zone 1: 2D cross-section view of (a) inverted Vs and (b) inverted Vp.

The 3D renderings of inverted Vs and Vp are shown in Figure 8. Both profiles exhibit the shallow soil layer underlain by the limestone layer. The top of the limestone varied from 6 to 8 m deep below the test area of 36 × 18 m. In Figure 8a, two separate voids are imaged in the middle of the medium near a depth of 15 m. One void intersected with the SPT borehole; the other void was offline and was not identified by the SPT. The Vp image (Figure 8b) is consistent with the Vs image. Compared with the 1D SPT, the seismic results provided much more information about the subsurface such as locations and 3D dimensions of voids and 3D soil/rock variation.

Figure 8.

Test Zone 1: 3D rendering of (a) inverted Vs and (b) inverted Vp.

Test Zone 2

Similar to Zone 1, SPT2 was conducted to a depth of 18.10 m at intervals of 0.45 m above the 6-m depth and 1.5 m below the 6-m depth. Seismic data were recorded for 102 blows from the ground surface to the maximum SPT depth, with 5 to 10 blows at each SPT interval. From about 5- to 9-m depth, the weight of the rod or weight of the hammer were recorded (SPT N-values = 0), suggesting that there was a void filled with air and some raveled soils. In addition, a complete loss of drilling fluid circulation happened at a depth of 17.1 m, suggesting another void.

Data analysis was performed in the same manner as that of Test Zone 1. The analyzed medium of 24 × 36 × 18 m (depth × length × width) with cells of 0.60 × 0.75 × 0.75 m (depth × length × width) were used. The 102 blows were assigned to 19 depths, and data from blows at the same depth were averaged to reduce noise and improve signal quality. The inversion started on the same initial model (Figure 4), used 10- to 40-Hz data, and stopped after 100 iterations.

The final inverted results are shown in Figure 9. The Vs profile (Figure 9a) consists of a soft layer with Vs∼200 m/s from 0 to about a 5-m depth, and a limestone layer below the 5-m depth. It contained a large, low-velocity anomaly (Vs < 100 m/s) at shallow depths (5 to 10 m) and a deeper low-velocity anomaly (Vs < 150 m/s) below a depth of 15 m. The Vp profile (Figure 9b) was similar to that of Vs, in relation to the soil/rock layers and low-velocity anomalies. For a better visualization, Figure 10 shows a 3D rendering of inverted Vs and -Vp. Both profiles exhibited the shallow soil layer underlain by the limestone layer. The large low-velocity anomaly is clearly imaged at a depth between 5 and 10 m, and the deep anomaly is also shown at about 17 m deep.

Figure 9.

Test Zone 2: 2D cross-section view of (a) inverted Vs and (b) inverted Vp.

Figure 10.

Test Zone 2: 3D rendering of (a) inverted Vs and (b) inverted Vp.

Comparison of S-wave velocity to SPT N-values

The comparison of inverted Vs to the SPT N-values at the two SPT locations is shown in Figure 11. For Test Zone 1, Figure 11a shows the initial and inverted Vs profiles at the SPT location, together with the SPT N-values (blow counts) of SPT 1. The Vs values appeared to change significantly during the inversion. At the void location, at a depth of 13 to 17 m, the cell Vs changed from the initial value of about 400 m/s to the final value of about 120 m/s. The final Vs profiles generally agreed with the SPT N-values, including a high-velocity zone about the void at a 10- to 12-m depth, and the low-velocity zone at the void location (13 to 17 m).

Figure 11.

Comparison of inverted and initial Vs with SPT N-values for (a) SPT 1 and (b) SPT 2.

For Test Zone 2, the initial and inverted Vs profiles at the SPT location are shown in Figure 11b, together with the SPT N-values of SPT 2. The good agreement between the inverted Vs and SPT N-values was observed from the ground surface to a 13-m depth. Both profiles showed low values from a depth of 5 to 10 m at the void location, and high values at about 3- and 13-m depths. Although there were discrepancies between Vs and SPT N-values below 13-m, the complete loss of drilling fluid circulation at a 17.1-m depth suggests that the SPT was next to the offline void in rock (shown as the deep anomaly in Figure 10).

However, the SPT N-values were more erratic than the Vs for both SPTs. This was owing to the SPT N-values representing the local soil properties at the device’s tip, whereas the Vs values represented the average material properties of a 0.60 × 0.75 × 0.75 m cell. In addition, the regularization used in Equation 5 tends to tie a cell to its six adjacent cells (i.e., top, bottom, left, right, front, and back), leading to smooth inverted models. Nevertheless, the SPT-seismic FWI was able to detect buried voids and characterize the subsurface soil/rock properties over large 3D volumes.

Finally, although the method has been only applied to the Florida geology comprising fine sand, silt, and karst limestone, it is expected to be applicable to any geomaterials, both above and below water tables. The water content only affects Vp values, not Vs values. In addition, using 3D geophones would potentially improve characterized results, as SPT-induced wavefields consist of both vertical and horizontal components.

Conclusion

Field experiments were conducted at a site in Florida that presents weathered, karst limestone, to investigate the capabilities of the SPT-seismic method in imaging voids in rock. The method is based on 3D FWI of SPT-induced wavefields with in-depth sources and surface receivers. The experiments were performed for two surface areas of 36 × 18 m, using a 2D grid of 72 vertical geophones to record data induced by an SPT at the area center of each area. The seismic results suggested that by using SPT-induced wavefields emitted within rock mass and dominated by body wave components, the method provided new imaging capabilities of subsurface structures in terms of accuracy and resolution with depth. Soil/rock properties and voids were characterized at submeter pixels over a large 3D domain of 24 × 36 × 18 m (depth × length × width) or up to 18 m away from the SPT location. Multiple voids at various depths from 5 to 17 m were detected with locations and dimensions consistent with the SPT data. The method’s ability of characterizing soil/rock in a 3D domain around a SPT borehole is significant for the efficient design of deep foundations. It allows use of a single SPT at the center of a pile group to characterize the entire volume of materials supporting all piles in the group. The method could help reduce the cost of invasive tests, as well as minimize the risk of pile failure/collapse as a consequence of unexpected site conditions.

Footnotes

Acknowledgements

We thank the Florida Department of Transportation state material office for their assistance with the SPT-seismic data.

Author Contributions

The authors confirm contribution to the paper as follows: study conception and design: K. T. Tran, M. Mirzanejad; data collection: K. T. Tran, S. J. Wasman; analysis and interpretation of results: K. T. Tran, M. Mirzanejad, D. Horhota; draft manuscript preparation: K. T. Tran, D. Horhota, S. J. Wasman. All authors reviewed the results and approved the final version of the manuscript.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was supported by the Florida Department of Transportation: grant no. BED31-977-01.

ORCID iDs

Majid Mirzanejad

Scott J. Wasman

References

Tran

K. T.

McVay

Faraone

Horhota

Sinkhole Detection Using 2D Full Seismic Waveform Tomography. Geophysics, Vol. 78, No. 5, 2013, pp. R175–R183. https://doi.org/10.1190/geo2013-0063.1.

Tran

K. T.

Mirzanejad

McVay

Horhota

3-D Time-Domain Gauss–Newton Full Waveform Inversion for Near-Surface Site Characterization. Geophysical Journal International, Vol. 217, No. 1, 2019, pp. 206–218. https://doi.org/10.1093/gji/ggz020.

Mirzanejad

Tran

K. T.

McVay

Horhota

Wasman

Sinkhole Detection with 3D Full Seismic Waveform Tomography. Geophysics, Vol. 85, No. 5, 2020, pp. B147–B157.

Pratt

R. G.

Goulty

N. R.

Combining Wave-Equation Imaging With Traveltime Tomography to Form High-Resolution Images From Crosshole Data. Geophysics, Vol. 56, No. 2, 1991, pp. 208–224. https://doi.org/10.1190/1.1443033.

Pratt

R. G.

Shipp

R. M.

Seismic Waveform Inversion in the Frequency Domain, Part 2: Fault Delineation in Sediments Using Crosshole Data. Geophysics, Vol. 64, No. 3, 1999, pp. 902–914. https://doi.org/10.1190/1.1444598.

Zhou

Greenhalgh

S. A.

Crosshole Seismic Inversion With Normalized Full-Waveform Amplitude Data. Geophysics, Vol. 68, No. 4, 2003, pp. 1320–1330. https://doi.org/10.1190/1.1598125.

Wang

Rao

Crosshole Seismic Waveform Tomography–I. Strategy for Real Data Application. Geophysical Journal International, Vol. 166, No. 3, 2006, pp. 1224–1236. https://doi.org/10.1111/j.1365-246X.2006.03030.x.

Cheng

Liu

Wang

Zong

Multi-Hole Seismic Modeling in 3-D Space and Cross-Hole Seismic Tomography Analysis for Boulder Detection. Journal of Applied Geophysics, Vol. 134, 2016, pp. 246–252. https://doi.org/10.1016/j.jappgeo.2016.09.014.

Rao

Wang

Chen

Wang

Crosshole Seismic Tomography With Cross-Firing Geometry. Geophysics, Vol. 81, No. 4, 2016, pp. R139–R146. https://doi.org/10.1190/geo2015-0677.1.

10.

Butchibabu

Sandeep

Sivaram

Y. V.

Jha

P. C.

Khan

P. K.

Bridge Pier Foundation Evaluation Using Cross-Hole Seismic Tomographic Imaging. Journal of Applied Geophysics, Vol. 144, 2017, pp. 104–114. https://doi.org/10.1016/j.jappgeo.2017.07.008.

11.

Y. C. H.

Danovan

The Potential of Seismic Cross-Hole Tomography for Geotechnical Site Investigation. In E3S Web of Conferences ( Tarantino

Ibraim

, eds.), EDP Sciences, 2019, p. 18006. https://doi.org/10.1051/e3sconf/20199218006.

12.

Mirzanejad

Tran

K. T.

McVay

Horhota

Wasman

Coupling of SPT and 3D Full Waveform Inversion for Deep Site Characterization. Soil Dynamics and Earthquake Engineering, Vol. 136, 2020, p. 106196. https://doi.org/10.1016/j.soildyn.2020.106196.

13.

Virieux

P-SV Wave Propagation in Heterogeneous Media: Velocity-Stress Finite-Difference Method. Geophysics, Vol. 51, 1986, pp. 889–901. https://doi.org/10.1190/1.1442147.

14.

Robertsson

J. O.

A Numerical Free-Surface Condition for Elastic/Viscoelastic Finite-Difference Modeling in the Presence of Topography. Geophysics, Vol. 61, 1996, pp. 1921–1934. https://doi.org/10.1190/1.1826736.

15.

Komatitsch

Martin

An Unsplit Convolutional Perfectly Matched Layer Improved at Grazing Incidence for the Seismic Wave Equation. Geophysics, Vol. 72, No. 5, 2007, pp. SM155–SM167. https://doi.org/10.1190/1.2757586.

16.

Nguyen

T. D.

Tran

K. T.

Site Characterization With 3D Elastic Full-Waveform Tomography. Geophysics, Vol. 83, No. 5, 2018, pp. R389–R400. https://doi.org/10.1190/geo2017-0571.1.

17.

Tran

K. T.

Hiltunen

D. R.

Inversion of Combined Surface and Borehole First-Arrival Time. Journal of Geotechnical and Geoenvironmental Engineering, Vol. 138, No. 3, pp. 272–280. https://doi.org/10.1061/(ASCE)GT.1943-5606.0000587.