*Article* **Longitudinal Monostatic Acoustic Effective Bulk Modulus and Effective Density Evaluation of Underground Soil Quality: A Numerical Approach**

**Yuqi Jin 1,2, Tae-Youl Choi <sup>2</sup> and Arup Neogi 1,3,\***


**Abstract:** In this study, we introduce a novel method using longitudinal sound to detect underground soil voids to inspect underwater bed property in terms of effective bulk modulus and effective density of the material properties. The model was simulated in terms of layered material within a monostatic detection configuration. The numerical model demonstrates the feasibility of detecting an underground air void with a spatial resolution of about 0.5 *λ* and can differentiate a soil firmness of about 5%. The proposed technique can overcome limitations imposed by conventional techniques that use spacing-consuming sonar devices and suffer from low penetration depth and leakage of the transverse sound wave propagating in an underground fluid environment.

**Keywords:** ultrasonic elastography; underground detection; soil inspection; underwater acoustics

#### **1. Introduction**

According to the US Geological Survey in 2014, the average cost of karst collapses in the United States over the past 15 years is more than \$300 million per year. The subsidence from sinkhole collapse is especially highest in Florida, Texas, Alabama, Missouri, Kentucky, Tennessee, and Pennsylvania. It is impossible to know when a catastrophic sinkhole collapse occurs. However, it is possible to predict the occurrence of such likely events. Sinkholes in karst terrain occur naturally and from anthropogenic activity, e.g., groundwater development, oil and gas drilling, surface loading, and urban expansion into previously undeveloped sinkhole-prone areas and drought or precipitation extremes [1,2]. Most states with substantial damage attributed to karst sinkholes have public resources documenting sinkholes and sinkhole density locations, except for Texas [3].

Hence, the appropriate geophysical methods to provide subsurface information are crucial for the migration of catastrophic disasters due to subsidence or sinkholes. Noninstructive tests or non-destructive (NDT) such as ground-penetrating radar (GPR) [4,5], spectral analysis of surface waves (SASW) [6,7], multi-channel analysis of surface waves (MASW) [8,9], and micro-tremor array measurement (MAM) [10,11] are useful methodologies to detect underground voids. They can provide 2D or 3D subsurface stiffness profiles from the measurements at the ground surface. Each method has advantages and limitations. For example, GPR can identify layering sites, but it cannot resolve material properties. However, seismic methods SASW and MASW can resolve layer thickness and stiffness of materials. These methods can be inserted into boreholes and can be used to measure subsurface characteristics from the inside of a borehole. Accurate voids detection can be appreciable. However, the softening process that occurs before the air voids formation was not usually involved in the checklist for inspection. The density and mechanical properties

**Citation:** Jin, Y.; Choi, T.-Y.; Neogi, A. Longitudinal Monostatic Acoustic Effective Bulk Modulus and Effective Density Evaluation of Underground Soil Quality: A Numerical Approach. *Appl. Sci.* **2021**, *11*, 146. https:// dx.doi.org/10.3390/app11010146

Received: 10 December 2020 Accepted: 23 December 2020 Published: 25 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/ licenses/by/4.0/).

undergo a clear decrease when the soil became soft. Unfortunately, those methods are limited to provide more in-depth engineering information, such as soil type, strength, stability, and so on [12].

Furthermore, these techniques do not resolve subsurface layering in the presence of certain anomalies. Impedance contrasts, moisture, and cavities can affect different tests in different ways. Ultrasonic techniques have been recently used to study the underwater distribution in soil [13] and soil properties. Recently pulsed velocity ultrasound has been used to detect hard objects in farmland [14]. However, most of these techniques require close contact of the transducer with the soil.

Therefore, it is necessary to find an appropriate non-contact method for mapping subsurface voids and monitoring the soil's healthy in terms of mechanical properties and density while providing material properties through the evaluation of geophysical methods. Electromagnetic and seismic monitoring systems are the most commonly used techniques to detect voids on land. To overcome the limitations in these methods, such as compactness, low penetrate depth, and leakage of the transverse sound wave propagating in underground fluid, we introduce a recently developed elastographic mapping technique. The effective bulk modulus and effective density detection (EBME) [15] have been applied to underground soil health monitoring and void detection in a compact monostatic setup.

#### **2. Numerical Experiment Design**

As Figure 1 shows, the typical underground is formatted in layers structure modeled as ambient air, soil layer with and without voids, and an underlying rock layer. The basic principle of the model employed in this work involves using low-frequency acoustic waves to detect the soil's effective density. It is based on the amplitude ratio of the reflected wave between the soil layer and the underlying layer due to acoustic impedance mismatch. The effective density can be presented in terms of an absolute value or a relative scale estimated from the recently invented non-invasive imaging technique: effective bulk modulus elastography (EBME) [15,16]. The previous studies used this technique to distinguish different materials such as hard and soft materials and similar tissue phantoms in terms of effective bulk modulus and effective density [15]. The application of EBME showed that the unique technique could differentiate the various regions of 3D printed plastic differing in density due to the air porosity introduced during the printing process under various conditions of printing. One of the fabricated samples had five density zones varying from 100% to 60%, which is similar to the varying packing density of porosity in soil due to environmental conditions. The effective density imaging technique remotely evaluated the absolute elastic values of the various regions with a maximum 6% error in absolute density values [17]. The technique can be applied to the underlying layers using acoustic radiation force to estimate the effective density in both lateral and axial directions [17]. This work motivated us to use this technique to study the void formation and the packing density in soil using a remote and rapid scanning technique applied to characterize other material systems.

## **Figure 1.** Material model of the layers used for simulations.

#### *2.1. Effective Bulk Modulus and Effective Density Calculation*

Effective bulk modulus and Effective density in scanned imaging were calculated as [15]:

*ρ* = *c*−<sup>1</sup> *Z*<sup>0</sup> ⎛ ⎜⎜⎝ <sup>−</sup><sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> − <sup>4</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> + 1 *p*1 *pe* − *p*<sup>0</sup> − 2 ⎞ ⎟⎟⎠, *Z Z*0 > 1, *ρ* = *c*−<sup>1</sup> *Z*<sup>0</sup> ⎛ ⎜⎜⎝ <sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> + <sup>1</sup> <sup>−</sup> <sup>4</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> *p*1 *pe* − *p*<sup>0</sup> + 2 ⎞ ⎟⎟⎠, 1 3 < *Z*1 *Z*0 < 1. *ρ* = *c*−<sup>1</sup> *Z*<sup>0</sup> ⎛ ⎜⎜⎝ <sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> − <sup>1</sup> <sup>−</sup> <sup>4</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> *p*1 *pe* − *p*<sup>0</sup> + 2 ⎞ ⎟⎟⎠ , 0 < *Z*1 *Z*0 ≤ 1 3 . (1) *K* = *c Z*<sup>0</sup> ⎛ ⎜⎜⎝ <sup>−</sup><sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> − <sup>4</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> + 1 *p*1 *pe* − *p*<sup>0</sup> − 2 ⎞ ⎟⎟⎠, *Z Z*0 > 1, *K* = *c Z*<sup>0</sup> ⎛ ⎜⎜⎝ <sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> + <sup>1</sup> <sup>−</sup> <sup>4</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> *p*1 *pe* − *p*<sup>0</sup> + 2 ⎞ ⎟⎟⎠, 1 3 < *Z*1 *Z*0 < 1. *K* = *c Z*<sup>0</sup> ⎛ ⎜⎜⎝ <sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> − <sup>1</sup> <sup>−</sup> <sup>4</sup> *<sup>p</sup>*<sup>1</sup> *pe* − *p*<sup>0</sup> *p*1 *pe* − *p*<sup>0</sup> + 2 ⎞ ⎟⎟⎠ , 0 < *Z*1 *Z*0 ≤ 1 3 . (2)

In the above equations, *pe*, *p*0, and *p*<sup>1</sup> are, respectively, the highest amplitude of the source pulse from the probe, the reflection of the wavefront back from the front interface of the sample layer, and the second echo back from the interface between the target sample layer and next material layer separately. *pe* was the maximum amplitude value of the emission source. It was set in the software to 1 μPa on the absolute scale. *p*<sup>0</sup> and *p*<sup>1</sup> were the maximum absolute values of the detected reflection signal amplitudes obtained from the probe's upper surface. The values were averaged from the ten linear distributed arrays on the probe line laterally. *c* is indicated the sound velocity from the time of flight in the sample layer at the measured location, described as *c* = 2*d*/(*t*<sup>1</sup> − *t*0), where *t*<sup>1</sup> and *t*0, are the first peak of the first and second echoes. *d* is the thickness of the target layer. Effective density *ρ* is calculated from *ρ* = *Z*/*c*, where *Z* is the acoustic impedance of the sample at the scanned location. The baseline impedance value was referred to in the previous layer is *Z*<sup>0</sup> = *ρ*0*c*0.

For the multiple layers of various materials:

$$p\_k = \frac{Z\_n}{Z\_0} (p\_\varepsilon - \text{sig}(Z\_1 - Z\_0) |p\_0|) \cdot \left(\prod\_{i=2}^k t\_{i-1,i}\right) r\_{k-1,k} \left(\prod\_{i=1}^k t\_{i,i-1}\right), k = 1, 2, \dots, n. \tag{3}$$

where the transmission and the reflection coefficients are *ti*−1,*<sup>i</sup>* = (2*Zn*)/(*Zn*−<sup>1</sup> + *Zn*) and *rk*−1,*<sup>k</sup>* = (*Zn* − *Zn*−1)/(*Zn* + *Zn*−1) the reflection coefficient at the interface between layer (*n* − 1) and *n*. The reflection coefficient of the interface between the last layer and ambient material is expressed as *rk*−1,0 = (*Z*<sup>0</sup> − *Zn*−1)/(*Zn*−<sup>1</sup> + *<sup>Z</sup>*0). *<sup>n</sup>* numbers of *Zn* values are obtained by solving *n* numbers of Equation (2) for the *n* numbers of layers in the samples. The effective density values are expressed as *ρ<sup>n</sup>* = *Znc*−<sup>1</sup> *<sup>n</sup>* .

#### *2.2. Numerical Modeling*

The numerical simulations were performed using COMSOL Multiphysics. The geometry was designed in two-dimension to reduce computational time. The whole detected region was eight meters in length and 4 m wide, including a 1.5 m layer of air thickness between the probe and soil. We also consider 2 m of the rocky layer under the thick layer of soil. The physical properties of the regular soil layer were defined as *c* = 800 m/s and *ρ* = 2000 kg/m3, the rock layer was presumed to have *c* = 2000 m/s and *ρ* = 3000 kg/m3. The room temperature speed of sound in the air was considered to be *c* = 342 m/s and density as *ρ* = 1.225 kg/m3. The physical properties of air and rock layers were provided by the built-in materials library in COMSOL Multiphysics software. The parameters related to soil properties were used from the literature [18,19]. The soft-soil was defined to exhibit a 5% reduction in the speed of sound and its density compared to regular soil. We also considered that the softer soil had a 5% decrease in the speed of sound and density from the soft soil. The time-dependent wave equation was simulated with a general pulse form with its pulse function expressed as sin(*ω*0*t*)*e*−*f*0(*t*−3*T*0) 2 , where *ω*<sup>0</sup> was the angular frequency of the pulse at the operating frequency *f*<sup>0</sup> (2000 Hz), and *T*<sup>0</sup> = 1/*f*<sup>0</sup> is the time period. *t* was the time interval over which the event was simulated. The time window used for the estimation of the wave propagation was 480*T*0.

Each of the simulated model illustrated in Figure 2 shows the geometrical configuration used in this study, which were considered to be 2.05 m tall (in the vertical direction) and 1 m wide (in the horizontal direction). The top probe was 0.05 m-thick and 0.5 m-wide. The rest of the 2 m region was generally separated into three major zones. From top to bottom, the air ambient layer was 0.15 m thick. The major soil zone thickness was 1.6 m, which had a 0.25 m rock layer under it. For the cases of soft soil, softer soil, and air void existing in Figure 2B–D, the center of the anomalous regions located at the center of the soil layer. Figure 2B,C, show that the soft soil and softer soil regions were 0.8 m wide and 0.5 m thick. In Figure 2D, the air void was a circle with 0.05 m in diameter.

**Figure 2.** Sound wave propagation through air, soil and rocky layers in four different cases. The red and blue color scale indicated the positive and negative amplitude of the sound wave pressure. Case (**A**): healthy soil. Case (**B**): healthy soil with 5% properties reduction for the soft soil region. Case (**C**): healthy soil with a softer soil region with 5% lower values compared to the soft soil in (**B**). (**D**): healthy soil with a small internal air void. The size of the air void is smaller than the sound wavelength.

#### **3. Results and Discussion**

A finite element analysis-based numerical simulation was performed to simulate the feasibility of using effective density detection in determining soil voids and the overall density of the porous soil. The simulation was used to visualize the transient sound wave propagation in the four soil with different conditions. As Figure 2 shows, the models of the initial study were categorized into healthy soil (A), soil with a region of soft soil Figure 2B, soil with a region of softer soil Figure 2C, and healthy soil with an internal small air gap (of about 0.5 *λ* in size) Figure 2D. In the axial propagation of the wave in the detection setup, the low-frequency sound wave pulse has a small amplitude reflection at the interface between Air and soil layers. In the soil layers, the sound wave propagation was delayed in the case of Figure 2B,C comparing with healthy soil (A). Moreover, in the case of Figure 2D, the small air void caused a scattering effect on the propagating wave without any clear temporal delay. The sound wave was reflected back into the transducer with a larger amplitude at the interface between the soil layer and the rock layer under it. The reflected wave propagated in the opposite direction of the wave emission direction. The backward trip of the wave undergoes another temporal delay in the case of Figure 2B,C, and once scattering effect in case (D). By measuring the reflected signal pressure over a roundtrip propagation of the wave, a significant difference was found between the four cases, as shown in Figure 3.

**Figure 3.** (**A**) Temporal reflected signal collected at the surface of the probe in all four cases mentioned in Figure 2. "No cavity" was the healthy soil condition in case (**A**). "One cavity" indicated the soil with a small air void in case (**D**). "Soft soil" was the case (**B**), which has a soft soil zone in the soil layer, which has similar properties. "Softer soil" was the case (**C**), which has a soft soil zone in the soil layer which has dissimilar properties. The time window width was 6 ms. The wave completed a roundtrip in the simulated model in the calculated time length. (**B**) was the zoomed-in view between 0.4 ms and 1.8 ms. (**C**) was the zoomed-in view between 2 ms and 4 ms. (**D**) was the zoomed-in view between 4.125 ms and 5 ms.

In Figure 3, Figure 3A shows the entire range of the temporal response of the propagating wave. Figure 3B–D shows the zoomed-in response from the various temporal regions plotted in Figure 3A. Figure 3B shows that the echoes from the two reflections can be expected from the interface between the air and soil and its second roundtrip envelope. Figure 3C depicts the zoomed-in time window, and the propagation of the wave occurs within the soil layer's internal region. No reflection pulse can be observed from the healthy soil line. The first and third echo on the green line occurred from the front and back interface between softer soil and healthy soil. The second and fourth pulses were the second roundtrip of the first and third reflection that occurred inside the softer soil cavity. The only two visible reflections on the soft soil were the echo from the two mismatch interface around the anomalous zone and occurred along the direction of the propagating wave. The red line has two lower amplitude echoes occurring from the air cavity's front and back interface on the wave propagation direction. In Figure 3D, there were three major reflection envelopes on each line. The first echo was reflected from the rocky layer's front interface, and the second echo was from the back boundary of the rock layer. The third reflection envelope occurred from the back surface and underwent a roundtrip internally within the rocky layer.

Figure 3A showed the temporal measurement of the reflected sound wave at the probe surface regarding the normalized sound pressure for all the four proposed cases. Before 2 ms, the sound wave has experienced the same media, and the first two echo were reflected at the interface between air and soil layers due to acoustic impedance mismatch. Between 2 ms and 4 ms, the temporal signal's zoomed-in view shows the reflections from the soft soil, softer soil, and air void in case Figure 3B–D compared to the normal soil ("No cavity"). The larger echo amplitude in case Figure 3C was due to the larger difference between the density and speed of sound properties in softer soil and healthy soil. The difference between the density and speed of sound in soft soil in case (Figure 3B) and healthy soil is smaller comparing with case Figure 3C (softer soil), which led to a smaller amplitude reflected echoes. The red dotted trace showed the two small symmetric reflection envelopes from the air void, as shown in Figure 3D. The air void's axial size could be estimated from the temporal distance between the highest peak of the two small symmetric reflection envelopes.

Figure 3D shows the signal between 4.125 ms and 5 ms, which indicated the reflected wave from the interface between the soil and rock layers. As the sound wave in the four cases experiences different soil conditions during the two roundtrip propagation, the reflection from the interface between the soil and rock layers shows a clear time difference between the cases. Comparing Figure 3A, which is for healthy soil, and Figure 3D with one small air void, the reflected wave does not have any time delay for the arrival pulse. The attenuated amplitude of the signal in case Figure 3D was affected by the scattering due to the small air void. Due to slightly lower density and speed of sound in the soil's temperate region, the reflected wave in the case, Figure 3B, had a small amount of time delay compared to the healthy soil. When the decreased density and sound velocity values occurred more in the soft soil zone, the time delay between the soil with and without the soft region was larger, as illustrated in Figure 3C.

From the temporal waveform in terms of pressure, as shown in Figure 3, the calculated relative density and effective bulk modulus values in healthy soil were normalized to 1. In the case of Figure 3B, the effective bulk modulus decreased to 0.923, and the effective density value of the soil layer reduced to 0.984. In the case of C, the effective bulk modulus decreased to 0.907, and the effective density value of the soil layer was lowered to 0.962. In the case of D, the effective bulk modulus decreased to 0.985, and the effective density value of the soil went down to 0.985. As the listed equations indicated, the amount of reflection time delay back from the rock layer provides a larger effect on the estimated effective bulk modulus and density values than the echo amplitude difference since sound speed is a square term in the equation. However, in Figure 3D, since the rock layer interface echo does not have a clear time delay in the measurement, the effective bulk modulus

and effective density value decreased linearly by a small amount based on a decrease of the echo amplitude due to the scattering effect. In the numerical simulation, the defined media was not dispersive. In reality, the dispersion [20,21] may be introduced by the nonuniform soil status. In a highly dispersive medium, the frequency components in a broadband acoustic pulse have different speed of sound (phase velocity). This dispersion would elongate the pulse width and modify the velocity of the acoustic pulses. For an accurate estimation of the effective bulk modulus and effective density values, the single frequency component phase velocity should be considered and used in place of the speed of sound value from the pulse envelope.

Unlike the existing seismic methods of soil void detecting techniques, EBME uses a monostatic low-frequency sound wave to estimate the effective bulk modulus and effective density over the entire soil layer's depth in the effective area as large as the sound probe beyond the long-wavelength limit. The technique can monitor the soil's health in terms of relative content of air, water (effective elasticity), and foreign objects (Appendix A). The significant advantage of the EBME void detection is the requirement of compact equipment set up on the ground. Almost all the existing sonic soil void detecting techniques use a bistatic setup. To place the emission source and receiver on the ground required more space. Some techniques require either sound wave emission source array or detector array on the ground, which the preparation is time-consuming. In addition, the alignment of the array is important to obtain accurate results. The non-flat ground condition would introduce nonnegligible uncertainty to the detection. Since the EBME technique uses a high penetrating low-frequency sound wave, the probe is not necessary to contact the ground. The air layer between the probe surface and soil does not decrease the signal-to-noise ratio a lot. In this way, the limitation from the non-flat ground surface condition could be overcome. A sound wave array system would be preferred as it provides a 2D surface scan. Using EBME, the 2D surface scan could be carried out through a raster scan since it is independent of the ground surface condition and does not require any contact with the ground surface. Most seismic methods apply transversal mode sound wave or radiational stress, which approaches a limitation of wave propagation in the fluid such as underground water. Since most fluids do not have a shear modulus to transmit transverse mode vibration, the existing seismic methods have difficulty determining the underground structure and properties once the underground water layer exists. The EBME technique uses longitudinal mode wave, which can propagate in both solid and fluid media to overcome the limitation of the underground fluid propagation issue compared to existing seismic methods. Ground-penetrating radar is another non-destructive method in soil void detection. The radar emitted an MHz microwave into the soil and collecting reflection. Since the electromagnetic wave has a shorter wavelength, the resolution of ground penetrating radar is usually appreciable. However, the increased resolution by small wavelength introduces large noise as a tradeoff. The wavelength of the low-frequency EBME technique can minimize noise since it is employing a long-wavelength sound. Instead of using a long-wavelength sound wave to detect the comparable size of soil void, the EBME technique uses a relatively scaled effective bulk modulus and effective density to estimate the volume fraction of the void inside the soil by using a rapid non-contact raster scan. Once the porous soil's target area was determined, a higher frequency sound wave detector [22], acoustic lens/collimator [23–26], or electromagnetic radar could be applied in the region to find detailed information of the specific voids for further interest. This work showed the initial feasibility of underground acoustic detection. Real experiments will be further performed as future works in the lab scaled-down condition and further in real condition with the influence factors studies such as vegetation on the ground surface and unknown solid objects in the soil layer.

#### **4. Conclusions**

This study proposed a novel method to detect underground soil voids and monitor the soil healthy in terms of effective bulk modulus and effective density demonstrated by numerical simulation. The technique can detect about 0.5 *λ* size air void in soil and

5% reduction in soil density due to the decrease in the measured effective bulk modulus and density compared to a healthy reference soil. The proposed technique would provide better penetration depth than electromagnetic methods, more compactness than the multi-detectors array systems, and better resolution than conventional sonic techniques. Compared to surface wave and shear wave techniques, this study's novel method can overcome the limitation of non-guided propagation of the transverse wave in underground water or other fluid.

**Author Contributions:** Conceptualization, A.N. and T.-Y.C.; methodology, T.-Y.C. and Y.J.; software, Y.J.; validation, Y.J., T.-Y.C. and A.N.; formal analysis, T.-Y.C.; investigation, Y.J.; resources, A.N.; writing—original draft preparation, Y.J., T.-Y.C. and A.N.; writing—review and editing, A.N.; visualization, A.N.; supervision, A.N.; project administration, A.N.; funding acquisition, A.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by an Emerging Frontiers in Research and Innovation (EFRI) grant from the National Science Foundation (NSF) Grant No. 1741677. The support from the infrastructure and support of the Center for Agile & Adaptive and Additive Manufacturing (CAAAM) funded through State of Texas Appropriation #190405-105-805008-220 is also acknowledged.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data availability from corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Appendix A**

**Figure A1.** (**A**–**D**) Various time points of an additional case with steel block embedded into the soil layer, which could be clearly recognized. (**E**) Time of flight signal of the case.

#### **References**

