AbstractThis work studied the feasibility of imaging a coupled fluid–solid system using the elastodynamic and acoustic waves initiated from the top surface of a computational domain. A one-dimensional system, in which a fluid layer was surrounded by two solid layers, was considered. The bottom solid layer was truncated using a wave-absorbing boundary condition (WABC). The wave responses were measured by a sensor located on the top surface, and the measured signal contained information about the underlying physical system. Using the measured wave responses, the elastic moduli of the solid layers and the depths of the interfaces between the solid and fluid layers were identified. A multilevel genetic algorithm (GA) combined with a frequency-continuation scheme to invert the values of sought-for parameters was employed. The numerical results showed that (1) the depths of solid–fluid interfaces and elastic moduli can be reconstructed by the presented method, (2) the frequency-continuation scheme improves the convergence of the estimated values of parameters toward their targeted values, and (3) a preliminary inversion using an all-solid model can be employed to identify if a fluid layer exists in the model by showing one layer with a very large value of Young’s modulus (with a value similar to that of the bulk modulus of water) and a mass density similar to that of water. Then the primary GA inversion method, based on a fluid–solid model, can be utilized to adjust the soil characteristics and fine-tune the locations of the fluid layer. If this work is extended to a 3D setting, it can be instrumental in finding unknown locations of fluid–filled voids in geological formations that can lead to ground instability and/or collapse (e.g., natural or anthropogenic sinkholes, urban cave-in subsidence, and so forth).IntroductionDetecting underground cavities (e.g., karstic cavities, caves, tunnels, and so forth) is a challenging task for geotechnical engineering projects due to the geological and hydrogeological complexity of the subsurface environment. Geological hazards, such as collapsing soil or urban ground collapse due to subsurface voids, could induce significant damage to infrastructure. Such a hazard is one of the major issues for land planning, infrastructure operation and maintenance, and disaster management. In addition, cavity evolution that occurs due to hydrogeological process may cause sinkhole collapse. For example, Florida’s karst environment involves active groundwater recharge to the Floridan Aquifer that erodes overburden soils; thus, a cavity could grow gradually, leading to a sudden collapse (cover-collapse type) or gradual subsidence (cover-subsidence type) (Fig. 1) (Beck 1986; Tihansky 1999; Xiao et al. 2016; Kim et al. 2020; Nam et al. 2020). These large underground water-filled cavities hidden below building and transportation infrastructure should be predetected so that prevention and mitigation measures are applied before catastrophic collapse or excessive ground settlement takes place.Cone penetration tests (CPTs) and standard penetration tests (SPTs) have been used to detect and characterize subsurface cavities in a cover soil layer, referred to as a raveled zone in karst areas (Nam et al. 2018; Nam and Shamet 2020). Although SPTs and CPTs provide a continuous subsurface profile—revealing soil type, resistance, and stratigraphy—these invasive methods collect only pointwise profiles of the subsurface environment. It is time-consuming and labor-intensive to apply CPTs and SPTs to wide areas multiple times.Numerous geophysical studies over the last decades have introduced nondestructive evaluation (NDE) methods for detecting and estimating the geometry of subsurface cavities. For example, ground-penetrating radar (GPR), microgravimetry, and resistivity imaging and their combinations can characterize an underground cavity based on its size and depth position (Fehdi et al. 2014; Kiflu et al. 2016). GPR is inefficient in highly heterogeneous media—such as backfills, moist clays with high soil conductivity, and saturated soils below the groundwater table. The microgravimetric method is effective in detecting shallow cavities (Butler 1984; Bishop et al. 1997), but its density contrast cannot be obtained clearly if the size of a void is relatively small compared with its depth. An electrical resistivity method has been applied to detect air-filled or water-filled cavities (Van Schoor 2002; Coşkun 2012), but it requires a large areal space for surveying and cannot detect small-scale irregularities in the geologic interfaces because it measures average resistivity values and is easily ruined by noise sources, such as piping, power lines, and house structures.On the other hand, elastodynamic imaging has been used widely in site characterization (Brown et al. 2002; Kallivokas et al. 2013). It uses elastodynamic waves that are initiated by wave sources on the ground surface and/or borehole sources, and then are reflected or refracted due to material heterogeneity within a medium. The resulting wave motions can be measured on the ground surface or at borehole seismic arrays. Inverse modeling, associated with elastic wave propagation analyses, provides promising results for identifying the material properties of the media. Inverse modeling—employing a partial differential equation (PDE)-constrained optimization method and the finite-element method (FEM) or spectral-element method (SEM)—has been utilized to infer the spatial distribution of material properties of host media at fine resolution (Kang and Kallivokas 2011; Fathi et al. 2015a). For example the geotechnical site characterization methodologies have been investigated in a soil domain that is truncated by perfectly matched layers (PML), in which waves decay and are not reflected off surrounding boundaries (Fathi et al. 2015b), using state-adjoint-control equation-based full-waveform inversion (FWI) approaches (Kang and Kallivokas 2010, 2011; Kallivokas et al. 2013; Fathi et al. 2015a, 2016; Kucukcoban et al. 2019). In particular, there has been progressive development of the material inversion method for site characterization from a simplified one-dimensional (1D) setting to a full three-dimensional (3D) setting followed by field experimental validation. Kang and Kallivokas (2010) investigated a new material inversion method for identifying the vertical distribution of the shear wave speed within a PML-truncated 1D soil column with a hypothesis of horizontal layering without considering geometrical damping. Kang and Kallivokas (2011) further developed their material inversion method for estimating the spatial distribution of the shear wave speed in a two-dimensional (2D) scalar wave setting at fine resolution without considering the realistic behaviors of elastic waves (i.e., vector waves). Kallivokas et al. (2013) continued the investigation of the new material inversion method to estimate the spatial distributions of both P- and S-wave speeds within a 2D linear, elastic, undamped PML-truncated solid setting in consideration of vector wave behaviors. They used an assumption that the material properties of the soils are uniform in the antiplane direction of the considered 2D plane. They validated their numerical method using field experimental data that used a line loading–like application of wave sources, which replicated the 2D plane-strain setting of the computational study. Fathi et al. (2015a) developed the material inversion method for inverting the spatial distribution of both P- and S-wave speeds in a 3D linear, elastic, undamped PML-truncated solid. Fathi et al. (2016) validated the result of the inverse modeling of Fathi et al. (2015a) using field experimental data in a 3D setting. In the experiment, a point wave source was employed on the ground surface, and its corresponding wave responses were measured on the surface and used as measurement data in the inverse simulations. Fathi et al. (2016) compared the inverted spatial distribution of the wave speeds with that of a reference solution from a conventional method, such as spectral analysis of surface waves (SASW).In addition, the Gauss–Newton-based FWI method, which is based on a gradient vector and a Hessian matrix, has been studied for characterizing the material profiles in a truncated two-dimensional solid domain (Tran and McVay 2012; Pakravan et al. 2016) for the application of site characterization. The PDE-constrained optimization also has been used in consideration of the boundary-element method (BEM), the computational efficiency of which is much greater than that of the FEM due to the reduction of the dimensionality, to detect the geometries of wave-scattering objects in host media. There have been studies of inverse scattering algorithms using BEM wave solvers, hinging on the moving boundary concept and the total derivative (Petryk and Mroz 1986), which allows computing the derivative of an objective functional with respect to the geometry variables of scattering objects (Guzina et al. 2003; Jeong et al. 2009). The extended finite-element method (XFEM) has been used as a wave solver in studies to identify cracks or air-filled voids in solid media because it can avoid expensive remeshing processes in inverse iterations while using the flexibility of the FEM for heterogeneous materials (Jung et al. 2013). In addition, a frequency-continuation scheme was devised to help the inversion solver address the multiplicity of solutions of subsurface imaging problems (Bunks et al. 1995; Kang and Kallivokas 2010; Fathi et al. 2015a). That is, when the convergence rate of the inversion solver is decreased due to a local minimum of an objective functional that is composed of measurement data corresponding to a given dominant frequency of excitation, another set of measurement data that are induced by excitation of a different dominant frequency is employed. Such a frequency shift can change the curvature of the objective functional so that a minimizer can escape a local minimum. In general, one may consider increasing the frequency from one inversion set to another because it has been reported that using a lower frequency leads to the inversion result of a lower resolution (due to a larger wavelength); then one can fine-tune it using a higher frequency (due to a smaller wavelength) (Bunks et al. 1995; Kang and Kallivokas 2010; Fathi et al. 2015a).Despite such recent developments in elastodynamic wave imaging techniques, few papers have demonstrated the feasibility of identifying the properties of a solid system that includes fluid-filled voids. Finite-difference method (FDM)-based FWI approaches have been presented, which employ all solid elements in an entire computational domain and detect the areas with smaller values of shear wave speeds Vs (e.g., about 100 m/s) than those of typical soils and rocks to indicate air-filled voids (Tran and McVay 2012; Tran et al. 2013; Mirzanejad et al. 2020). However, those works have not explicitly modeled the interfaces between air voids and the solid domain. Thus, the authors are concerned that, first, modeling a fluid domain as solid elements with a small value of Vs could introduce numerical error in forward wave solutions. Second, using a small value of Vs requires a small solid element for wave simulations. However, the aforementioned FDM works did not use such a small element to model a void, so they could suffer from the additional error of wave solutions. To address this issue, the authors suggest explicitly using fluid elements to model voids filled with water (or air) so that accurate fluid–solid coupling (Everstine 1997; Lloyd et al. 2016) can be incorporated into the wave solver. Then, using such a wave solver, inverse modeling could identify accurately the interfaces between fluid and solids as well as the material profiles of solids. As a prototype work of the suggested method, this paper investigates the feasibility of estimating the interfaces between fluid and solid layers and the material properties of solid layers in a multilayered system using acoustic and elastodynamic waves generated from the ground surface. This investigation proved the feasibility of detecting a water-filled void and identifying its location and geometry in soils, due to which overlaying ground surface potentially could sink.This work used the FEM to study the wave responses of the coupled fluid–solid media in a prototype one-dimensional setting. The solid formation of a semi-infinite extent was truncated by using an absorbing boundary, and a proper fluid–solid interface condition was considered. The inverse problem was cast into a minimization problem, in which the value of an error between the synthetic measured wave response due to target profiles and the computed counterpart due to guessed profiles was minimized. The solution multiplicity of the minimization process was addressed using a multilevel genetic algorithm (GA) in this work. Under this new scheme, the authors suggest the following steps: (1) the optimizer identifies an inversion solution around a local minimum after a number of GA iterations in the first-level GA; (2) a misfit function then can be r-calculated for a different measurement signal induced by a pulse signal of a different central frequency in the next-level GA; (3) by doing so, the GA optimizer is able to escape the local minimum of the previous-level GA; and (4) the same technique is employed to address the solution multiplicity of the next-level GA. The numerical results show that the numerical simulations using the multiple-frequency-level GA accurately detects the interfaces between fluid and solid layers and the material properties of solid layers. The authors also demonstrate that a preliminary GA inversion can determine whether the subsurface media include a fluid layer. Specifically, the preliminary inversion can indicate the existence of a fluid layer by detecting a very large value of Young’s modulus (with a value similar to that of the bulk modulus of water) and a mass density with a value similar to that of water. When the existence of a fluid layer is detected, a primary GA inversion using the fluid–solid model can estimate further the locations of a fluid layer and the material properties of surrounding solids.Problem DefinitionThis study identified the depths of the interfaces between a fluid layer and its surrounding solid layers and the stiffness of the solid layers using a dynamic test, which generates elastodynamic waves into the subsurface media from the top surface. A 1D fluid–solid model, which consists of three different layers, i.e., two solid layers and one fluid layer (Fig. 2), was considered. The governing differential equations for the compressional waves in the three layers are (1) ∂∂x(Es1∂us1∂x)=ρs1∂2us1∂t2,0≤x≤Ls1(2) ∂2P∂x2=1cf2∂2P∂t2,Ls1≤x≤Lf(3) ∂∂x(Es2∂us2∂x)=ρs2∂2us2∂t2,Lf≤x≤Ls2where Ls1, Lf, and Ls2 = depth of interface between top solid layer and fluid layer, depth of interface between fluid and bottom solid layer, and depth of truncation wave absorbing boundary, respectively; us1(x,t) and us2(x,t) = displacement fields of wave motions of solid particles in top and bottom solid layers; and P = fluid pressure field of wave motions in fluid. The wave speeds of the three layers are given by cs1=Es1/ρs1, cf=κf/ρf, and cs2=Es2/ρs2, respectively, where E, ρ, and κ are the modulus of elasticity, the density, and the bulk modulus. The problem of interest also includes the following boundary conditions: (4) Es1∂us1∂x=−f(t),x=0(5) ∂us2∂x=−1cs2∂us2∂t,x=Ls2A pulse signal f(t) is applied at the top of the solid layer at x=0 [Eq. (4)]. The depth of the bottom solid layer is truncated by using the absorbing boundary condition [Eq. (5)].The solid and fluid layers are coupled via the following interface conditions: (6) (7) ∂P∂x=−ρf∂2us1∂t2,x=Ls1(8) (9) ∂P∂x=−ρf∂2us2∂t2,x=LfThe preceding interface conditions originate from the dynamic interface condition between a solid and fluid in a multidimensional setting, presented by Everstine (1997), which is composed of the following. First, the effect of fluid pressure on the structure is imposed as a load. This corresponds to Eqs. (6) and (8). Second, the effect of structural motion on the fluid should be modeled as ∇P·n=−ρf[(∂2us)/∂t2]·n, where n is the normal vector at the interface. This corresponds to Eqs. (7) and (9).The system initially is at rest (10) (11) ∂us1∂t=∂us2∂t=∂P∂t=0,t=0Finite-Element Modeling to Compute the Wave ResponsesTo derive the weak form, the governing wave equations in the strong form Eqs. (1)–(3) are multiplied by the test functions v1(x), w(x), and v2(x). By integrating them by parts and using the boundary and interface conditions of the strong form, the following weak forms are obtained: (12) ∫0Ls1[Es1∂v1∂x∂us1∂x+ρs1v∂2us1∂t2]dx=v(0)f(t)−v(Ls1)P(Ls1)(13) ∫Ls1Lf[κf∂w∂x∂P∂x+ρfw∂2P∂t2]dx=ρfκf[w(Ls1)∂2us1∂t2(Ls1)−w(Lf)∂2us2∂t2(Lf)](14) ∫LfLs2[Es2∂v2∂x∂us2∂x+ρs2v2∂2us2∂t2]dx=v2(Lf)P(Lf)−Es2ρs2Es2v2(Ls2)∂us2∂t(Ls2)Then the functions are approximated as (15) v1(x)=v1Tϕ(x),us1(x,t)=ϕ(x)Tus1(t)(16) w(x)=wTΨ(x),P(x,t)=Ψ(x)TP(t)(17) v2(x)=v2TΩ(x),us2(x,t)=Ω(x)Tus2(t)where ϕ(x), Ψ(x), and Ω(x) denote vectors of global basis functions constructed by local shape functions; and us1, P, and us2 are vectors of nodal solutions. Introducing the finite-element approximations Eqs. (15)–(17) into the weak forms Eqs. (12)–(14) provides the following discrete form: (18) Ks1us1+Ms1∂2us1∂t2=f−L1P(Ls1)(19) KfP+Mf∂2P∂t2=L2ρfκf∂2us1∂t2−L3ρfκf∂2us2∂t2(20) Ks2us2+Ms2∂2us2∂t2=−ρs2Es2L5∂us2∂t+L4Pwhere the matrixes are defined as (21) Ks1=∫0Ls1Es1∂ϕ(x)∂x∂ϕ(x)T∂xdx,Ms1=∫0Ls1ρs1ϕ(x)ϕ(x)T(22) Kf=∫Ls1Lfκf∂Ψ(x)∂x∂Ψ(x)T∂xdx,Mf=∫Ls1LfρfΨ(x)Ψ(x)T(23) Ks2=∫LfLs2Es2∂Ω(x)∂x∂Ω(x)T∂xdx,Ms2=∫LfLs2ρs2Ω(x)Ω(x)Tdx(24) L1=ϕ(Ls1)ΨT(Ls1)=[0⋯0⋮⋱⋮1⋯0],L2=Ψ(Ls1)ϕT(Ls1)=[0⋯1⋮⋱⋮0⋯0](25) L3=Ψ(Lf)ΩT(Lf)=[0⋯0⋮⋱⋮1⋯0],L4=Ω(Lf)ΨT(Lf)=[0⋯1⋮⋱⋮0⋯0](26) L5=Ω(Ls2)ΩT(Ls2)=[0⋯0⋮⋱⋮0⋯1]and the force vector is f=ϕ(0)f(t). The discrete Eqs. (18)–(20) lead to the following coupled discrete form: (27) [Ks1L100Kf00−L4Ks2]⏟K[us1Pus2]+[00000000Es2ρs2Es2L5]⏟C[∂us1∂t∂P∂t∂us2∂t]+[Ms100−ρfκfL2MfρfκfL300Ms2]⏟M[∂2us1∂t2∂2P∂t2∂2us2∂t2]=[f00]which, at a discrete time tn, can be written (28) where qn=[us1T,PT,us2T]T is the vector containing the nodal solutions of solid displacements and acoustic pressures at the nth time step; and q˙n and q¨n denote the first- and second-order derivatives of qn with respect to time. The second-order ordinary differential [Eq. (28)] is solved using Newmark’s unconditionally stable time-integration scheme (Newmark 1959).Inverse ModelingThis study identified the four control parameters (Es1, Es2, Ls1, and Lf) by using measured wave responses initiated by a known excitation signal f(t). The values of the other variables (Ls2, ρs1, ρs2, κf, and ρf) were set to be known during the inversion process, and the excitation signal f(t) also was known. This problem was formulated into a minimization problem to identify estimated control parameters that could lead to the minimum of an objective functional (29) L=∫0T(um(0,t)−u(0,t))2dtwhere um(0,t) = wave signal recorded by sensor at x=0 due to target profile of control parameters; and u(0,t) = computed counterpart at the same measurement location due to guessed profile. In this computational study, um(0,t) was obtained synthetically by running the presented FEM solver using a target profile as an input. To prevent an inverse crime from taking place, smaller values of element sizes were used when um(0,t) was computed than when u(0,t) was computed.This work used a GA, which is a heuristic algorithm, to reconstruct the control parameters that are the fittest to the objective of the problem. In each generation of GA there are individuals, each of which has a set of different values of control parameters. In this work, the GA runs the FEM wave solver for the estimated values of control parameters of each individual to compute u(0,t) and evaluate the misfit [Eq. (29)]. Then the GA evaluates the fitness of each individual using the value of its misfit. The GA then weeds out less-fitting individuals when it creates the next generation of individuals, each of which is characterized by its own set of estimated parameters. In this process, by virtue of the mutation process, in which the GA learns the better-fitting characteristics of control parameters, the GA gives rise to new individuals in the next generation. The GA calculates the sensitivity of the fitness function with respect to the changes of the values of the control parameters and determines how the next generation of individuals is created. Namely, the GA repeats the tasks of weeding out, selecting the best individual, and mutating at each generation. Therefore, the suggested inverse modeling solver could identify targeted values of the parameters using the iterative process of the GA, in which the values of the estimated parameters of the best-fit individual evolve over the progress of generations.In the overall GA process in this problem, constraints were imposed on the allowable ranges of each control parameter. For example, Ls1 (the depth of the upper face of a fluid layer) was set to be always smaller than Lf (the depth of the lower face) so that the thickness of the fluid layer always was positive. In addition, the values of elastic moduli Es1 and Es2 were set to be always positive.It is known that in the minimization problem of a small number (e.g., less than 20) of control parameters, a GA is likely to find a set of control parameters that are close to the global minimum (Jeong et al. 2017; Guidio and Jeong 2021). Because there are only four control parameters, the authors originally hypothesized that the GA is a suitable solution approach to finding the global minimum solution of this inverse problem. However, the numerical experiments showed that the GA suffers from the solution multiplicity, and, thus, the authors tested a new multiple frequency–level GA approach. Namely, after the first level of the GA is finished, the presented inversion solver uses the reconstructed values from the first level to define the upper and lower limits of the control parameters in the next GA level, which uses a pulse signal of a different central frequency from its predecessor, creating a new synthetic um. In the final-level GA, the final best-fit estimated control parameters are obtained as the inversion solution.Numerical ExperimentsThis section presents a set of numerical experiments studying the performance of the presented multilevel GA–based parameter-estimation method. In the first two examples, three targeted fluid–solid models, which differed from each other in terms of Ls1, Lf, and Ls2, were considered as follows. Model 1 had a total length Ls2 of 60 m and its fluid was located between x of Ls1=20 m and x of Lf=25 m. The corresponding geometry parameters of Model 2 were Ls2 of 60 m, Ls1 of 40 m, and Lf of 45 m; and those of Model 3 were Ls2 of 120 m, Ls1 of 50 m, and Lf of 80 m. Parameters Ls1 and Lf are unknown parameters in the inversion. In all the models, the mass density (ρ) of the solid layers was 2,000 kg/m3; the Young’s moduli were Es1=2×108 Pa and Es2=5×108 Pa; and the bulk modulus (kf) and mass density (ρf) of the fluid layer were 2.34×109 N/m2 and 1,021 kg/m3, respectively. Whereas Es1 and Es2 were set to be unknown, ρ, kf, and ρf are set to be known in the inversion.In the presented numerical experiments, a three-level GA–based inversion method was used with a frequency-continuation scheme. Namely, in each GA set, a dynamic force with a different frequency content was used as follows: •The upper and lower bounds of each control parameter were set for the first-level GA, and its last-generation led to the best-fit individual.•In the second- and third-level GAs, the upper and lower limits of the four control parameters were updated with respect to the final-reconstructed values of the parameters obtained in the previous-level GA. Namely, in the second-level GA, the values of the estimated Ls1 and Lf were bounded using the same upper and lower limits as in the first GA level, whereas the values of the estimated Es1 and Es2 were bounded using ±50% derivations of their final reconstructed values obtained during the first GA level.•In the third GA level, the values of estimated Ls1 and Lf were bounded using ±5% derivations of their reconstructed values from the second GA level, whereas ±10% derivations of the values of estimated Es1 and Es2 that were reconstructed from the second-level GA level were used as their bounds in the last GA level.•The number of generations (GN) and population size (PS) both were 50 in all three GA levels.The computational procedure of the presented multilevel GA–based inversion method is summarized in Algorithm 1.Algorithm 1. Multilevel GA–based parameter-estimation algorithmAlgorithm 1. Multilevel GA–based parameter-estimation algorithm1: Set GN and PS=50.2: Set values for known parameters (Ls2, ρs1, ρs2, κf, and ρf).3: for the first-level GA do4: Set the central frequency of an excitation signal f(t).5: Set the upper and lower limits for each control parameter.6: Run the first GA level.7: return the reconstructed value of each control parameter.8: end for9: for the second-level GA do10: Set the central frequency of an excitation signal f(t).11: Set the upper and lower limits of Ls1 and Lf the same as those of the first GA level.12: Set upper and lower limits of Es1 and Es2 using ±50% derivations of their reconstructed values obtained from the first GA level.13: Run the second GA level.14: return the reconstructed value of each control parameter.15: end for16: for the third-level GA do17: Set the central frequency of an excitation signal f(t).18: Set the upper and lower limits of Ls1 and Lf using ±5% derivations of their reconstructed values obtained from the second GA level.19: Set the upper and lower limits of Es1 and Es2 using ±10% derivations of their reconstructed values obtained from the second GA level.20: Run the third GA level.21: return the reconstructed value of each control parameter.22: end forTo avoid an inverse crime, to compute um induced by targeted control parameters, the fluid–solid domain was discretized using an element size of 0.1 m, whereas an element size of 0.2 m was used to compute u due to estimated control parameters. In the forward and inverse modeling, the time step was 0.001 s, and the total observation duration T was set to 1 s.Three examples of numerical experiments are presented. The first example tested the inversion performance of the three-level GA approach by using a dynamic force of a Ricker pulse signal with a central frequency of 5, 10, and 15 Hz, respectively, in each level. In the second example, the inversion performance was examined by using a decreasing frequency-continuation counterpart of 15, 10, and 5 Hz. The last example investigated the utilization of preliminary inversion using all-solid layers for detecting the presence of a fluid–filled cavity by identifying E and ρ of all the solid layers, among which one layer had a very large value of E (with a value similar to that of the bulk modulus of water) and ρ of the value of the mass density of water. Once the preliminary inversion is completed, the optimizer uses the presented three-level GA approach based on the fluid–solid wave model as a primary inversion.To analyze the inversion results, the error between each target control parameter (Ls1, Lf, Es1, and Es2) and its estimated counterpart of the fittest individual at each generation is calculated as (30) E=|targeted value−estimated value||targeted value|×100[%]The average error of all the control parameters of the best-fit individual at each generation is calculated as (31) where Ek = error [Eq. (30)] of reconstruction of kth control parameter.Exemplary Forward Wave ResponsesPrior to the study of the inversion performance, an exemplary forward wave response in the computational domain induced by a pulse loading at the top surface is shown considering the fluid–solid Model 3 with Ls1=50 m, Lf=80 m, Ls2=120 m, Es1=2×108 N/m2, and Es2=5×108 N/m2. The time-dependent value of the dynamic force applied at the top surface is a Ricker wavelet with a central frequency of 20 Hz and a peak amplitude of 1,000 N/m2.The wave responses over space and time are shown in Fig. 3. This seismogram shows the stress field in the solid layers and the acoustic pressure in the fluid layer. The seismogram shows the following behaviors: 1.The elastic wave propagates throughout the first solid layer from the top surface and is transmitted via the first solid–fluid interface at x=50 m to the fluid layer (wave response around x=50 m and t=0.2 s). As the wave enters into the fluid layer, it reflects off the interface back to the solid layer at the same time.2.The acoustic pressure wave in the fluid layer is transmitted into the second solid layer through the second solid–fluid interface at x=80 m and reflects off the interface back to the fluid layer (response around x=80 m and t=0.2 s).3.The stress wave in the second solid layer is transmitted through the absorbing boundary at x=120 m without any reflection (wave response around x=120 m and t=0.3 s).The system repeats Behaviors 1–3 until the amplitude of the wave fades away.In the presented inverse modeling, the sensor on the top surface records the response signal in the time domain. The amplitudes and timings of the recorded signal may provide the inversion solver with information about the material properties of solid layers and the locations of fluid–solid interfaces. In other words, the timings of particular parts of the signal could indicate, primarily, the locations (Ls1 and Lf) of the fluid–solid interfaces, whereas the amplitudes of particular parts of the signal could indicate, primarily, the material properties (Es1 and Es2) in the solid layers.Example 1: Investigating the Inversion Performance by Increasing the Dominant Frequency of Excitation for Each GA LevelIn this example, the performance of the presented multilevel GA was studied for identifying targeted control parameters by increasing the frequency of an excitational Ricker signal for each GA level. That is, a dynamic pulse with a dominant frequency 5 Hz was used in the first-level GA, and, then 10 and 20 Hz pulses were used in the next levels. The control parameters to be identified were Ls1, Lf, Es1, and Es2. Cases 1–3, which used the fluid–solid Models 1, 2, and 3, respectively, were considered. •In Case 1, the targeted values of Ls1 and Lf were 20 and 25 m, respectively, whereas the values of their estimated counterparts, during the first-level GA level, were bounded at 10 m≤Ls1≤22 m and 23 m≤Lf≤35 m. The upper and lower limits were changed for the second and third GA levels as previously discussed.•In Case 2, the targeted values of Ls1 and Lf were set to 40 and 45 m, respectively, whereas the values of their estimated counterparts during the first GA level were bounded at 30 m≤Ls1≤42 m and 43 m≤Lf≤55 m.•In Case 3, the estimated values of Ls1 and Lf were bounded at 40 m≤Ls1≤65 m and 66 m≤Lf≤90 m in the first-level GA, and their targeted counterparts were 50 and 80 m, respectively.For all cases, the targeted values of Young’s moduli were set to Es1=2×108 Pa and Es2=5×108 Pa, whereas the estimated values for both moduli during the first GA level were bounded at 1×108 Pa≤E≤7.5×108 Pa. Their upper and lower limits were changed for the second and third GA levels as previously discussed.Table 1 presents the reconstructed control parameters of the fittest individual that was obtained at the last generation of each GA level per case. The table also includes the average error E¯ [Eq. (31)] of all the parameters reconstructed at the end of each GA level and the error E [Eq. (30)] of each control parameter of the fittest individual at the end of the last-level GA. In all Cases 1–3, the value of E¯ at the end of the last-level GA was smaller than that at the end of the first-level GA, and the final values of E¯ in all Cases 1–3 were smaller than 1%.Table 1. Example 1: reconstructed values of control parameters for each GA level obtained by increasing frequency of dynamic force for each GA levelTable 1. Example 1: reconstructed values of control parameters for each GA level obtained by increasing frequency of dynamic force for each GA levelCaseGA levelLs1 (m)Lf (m)Es1 (Pa)Es2 (Pa)Average error, E¯ (%)Case 1Target20252×1085×108—First20.424.92.13×1085.66×1085.53Second19.824.31.97×1084.94×1081.63Third20.024.82.00×1085.08×1080.60Individual error (%)0.000.800.001.60—Case 2Target40452×1085×108—First41.746.22.20×1085.86×1088.53Second40.646.02.06×1085.16×1082.48Third39.944.71.98×1084.97×1080.63Individual error (%)0.250.671.000.60—Case 3Target50802×1085×108—First53.190.02.28×1085.49×10810.63Second48.979.21.91×1085.33×1083.58Third49.979.91.98×1085.06×1080.63Individual error (%)0.200.121.001.20—The average error for the best-fit individual at each generation tended to decrease during the presented inversion process (Fig. 4). After about 20 GA iterations in the first-level GA, the value of E¯ converged but still had room for improvement (i.e., the optimizer found an inversion solution around a local minimum). To address such solution multiplicity in the first-level GA, a misfit function was used for um induced by a pulse signal of a new central frequency in the next-level GA. By doing so, the optimizer in the second-level GA can escape the local minimum of the first-level GA. By applying the same technique in the third-level GA, the solution multiplicity of the second-level GA can be resolved. Therefore, the authors suggest that it is feasible to reconstruct control parameters effectively using a multilevel GA process combined with a frequency-continuation scheme, by which the frequency is increased progressively for each GA level.Fig. 5 presents the details of the inversion performance of Case 3. Namely, Fig. 5 shows the histograms of the estimated control parameters of the entire population of the individuals during all the GA levels in Case 3. During the initial few generations of each GA level, the inversion solver explored a broader range of estimated values of control parameters. After these initial generations, the values of estimated parameters of individuals converged.Wave responses um induced by the target parameters matched u due to the final reconstructed control parameters, obtained at the end of each GA set, at the sensor on the top surface when a Ricker signal with a central frequency of 5, 10, or 20 Hz was employed in Case 3 (Fig. 6).Example 2: Investigating the Inversion Performance by Decreasing the Dominant Frequency of Excitation for Each GA LevelThis example studied the performance of reconstructing the control parameters using the presented multilevel GA combined with a frequency-continuation scheme that decreases the dominant frequency of the dynamic pulse over the multiple GA levels. The first-, second-, and third-level GAs used input force signals of central frequencies of 20, 10, and 5 Hz, respectively. Cases 4–6, which used Models 1, 2, and 3, respectively, were examined. The upper and lower limits of estimated control parameters for each GA level were set the same as those in Example 1. Table 2 presents the values of final reconstructed control parameters; the average errors for the best-fit individual of all Cases 4–6 decreased over the generations and the GA levels (Fig. 7). The final values of E¯ in all Cases 4–6 were smaller than 2% in Example 2. Thus, the authors suggest that the presented multilevel GA–based optimizer combined with a frequency-continuation scheme, which decreases the frequency for each GA level, is able to identify the values of the targeted control parameters as effectively as the increasing-frequency counterpart. However, this is the case only in the presented example because the presented example contains a small number of control parameters. In a more complex case (e.g., 2D or 3D settings), in which there are a large number of control parameters, the decreasing frequency-continuation scheme may not be as effective as its counterpart of increasing frequency.Table 2. Example 2: reconstructed values of control parameters for each GA level obtained by decreasing frequency of dynamic force for each GA levelTable 2. Example 2: reconstructed values of control parameters for each GA level obtained by decreasing frequency of dynamic force for each GA levelCaseGA levelLs1 (m)Lf (m)Es1 (Pa)Es2 (Pa)Average error, E¯ (%)Case 4Target20252×1085×108—First21.427.32.29×1085.55×10810.43Second20.124.92.04×1085.18×1081.63Third20.225.72.02×1085.04×1081.40Individual error (%)1.002.801.000.80—Case 5Target40452×1085×108—First41.045.72.11×1086.01×1087.44Second40.646.32.05×1085.09×1082.17Third40.044.72.00×1085.01×1080.22Individual error (%)0.000.670.000.20—Case 6Target50802×1085×108—First57.582.32.63×1084.51×10814.79Second49.777.91.98×1084.88×1081.66Third49.779.91.99×1085.18×1081.21Individual error (%)0.600.120.503.60—Example 3: Utilization of a Preliminary Inversion Using Only Solid Layers for Detecting the Presence of a Fluid Layer, Followed by a Primary Inversion Using Fluid–Solid LayersThis example presents a preliminary inversion using only solid layers for identifying if a target model contains a fluid–filled cavity. Here, three target models were used: one which did not contain a fluid layer, and two in which included a fluid-filled cavity. For all models, the preliminary GA parameter-estimation method assumed that the target contained only solid layers. After the preliminary inversion was completed, the authors continued a primary inversion using a fluid–solid model in a manner similar to that described in Example 1. Such a two-step (preliminary to primary inversion) approach was applied to Cases 7–10 •Cases 7 and 8 considered a target Model 4, which did not include a fluid layer, and contained only three solid layers of 5-m length each. The targeted Young’s modulus of each layer was Es1=2×108 Pa, Es2=3×108 Pa, and Es3=5×108 Pa. For Case 7, the targeted mass densities of all the layers were the same ρs1=ρs2=ρs3=2,000 kg/m3, whereas, for Case 8 they all were different from each other: ρs1=1,500 kg/m3, ρs2=2,000 kg/m3, and ρs3=2,200 kg/m3.•Case 9 considered Model 5, which consisted of three layers, i.e., two solid layers and one fluid layer. Model 5 had a total length Ls2 of 15 m, and its fluid was located between x of Ls1=5 m and x of Lf=10 m. The Young’s moduli of the solid layers were Es1=2×108 Pa and Es2=5×108 Pa, and their mass densities were ρs1=ρs2=2,000 kg/m3. The bulk modulus (kf) and mass density (ρf) of the fluid layer were 2.34×109 N/m2 and 1,021 kg/m3, respectively.•Case 10 considered Model 6, which consisted of a fluid layer surrounded by solid layers. The fluid layer was located between x of Ls1=4 m and x of Lf=9.3 m, and its Ls2 was 15 m. The properties of the fluid layer (kf and ρf) and the Young’s moduli of the solid layers (Es1 and Es2) were the same as those in Case 9, except that the mass densities of solid layers were ρs1=1,900 kg/m3 and ρs2=2,100 kg/m3.Cases 7 and 8 studied the capability of the preliminary inversion to detect the absence of a fluid layer by estimating E and ρ of all solid layers. Cases 9 and 10, by using only solid layers during the preliminary inversion, investigated the values of final-converged E and ρ from the preliminary inversion that indicate the presence of a fluid–filled cavity in the model.Preliminary Inversion Using All-Solid-Layer ModelFor all Cases 7–10, the preliminary inversion method utilized only one frequency level of GA, in which a Ricker pulse signal with a central frequency of 50 Hz was used as the dynamic force, and the number of generations and population size were 50 and 100, respectively. Moreover, the estimated values for Young’s moduli in all solid layers were bounded at 1×108 Pa≤E≤50×108 Pa, and the estimated values of their ρ were bounded at 1,000 kg/m3≤ρ≤3,000 kg/m3. Furthermore, the preliminary inversion assumed that the target contained three solid layers of 5-m length each.Fig. 8 shows the reduction of average errors for the best-fit individual for Cases 7–10 as the number of GA generations increased in the preliminary inversion. Table 3 presents the targeted and reconstructed control parameters from the preliminary inversion for Cases 7–10. For Cases 7 and 8, the results show that the preliminary GA-based inversion algorithm determined that all the layers in the target model were solid layers because the estimated Young’s moduli all were 2– 5×108 Pa in all solid layers, and the estimated mass densities all were 1,500– 2,200 kg/m3. The average error in Cases 7 and 8 was 3.82% and 4.82%, respectively. In contrast, in Cases 9 and 10, the preliminary inversion led to estimated values of Es2 and ρs2 of the second solid layer such that it could be interpreted as a fluid layer. Namely, the estimated values of Es2, 22.95×108 Pa (Case 9) and 24.32×108 Pa (Case 10), were similar to the bulk modulus of water (23.4×108 Pa), and the estimated values of ρs2, 1,022 kg/m3 (Case 9) and 1,070 kg/m3 (Case 10), were similar to the mass density of water (1021 kg/m3). The average errors in Cases 9 and 10 were 2.16% and 17.13%, respectively. Thus, after detecting the fluid layers in Cases 9 and 10, the authors continued the primary inversion in Cases 9 and 10 as follows.Table 3. Example 3: reconstructed values of control parameters obtained by preliminary inversionTable 3. Example 3: reconstructed values of control parameters obtained by preliminary inversionCasePreliminary inversionEs1 (Pa)Es2 (Pa)Es3 (Pa)ρs1 (kg/m3)ρs2 (kg/m3)ρs3 (kg/m3)Average error, E¯ (%)Case 7Target2×1083×1085×1082,0002,0002,000—Preliminary values2.01×1082.88×1085.26×1081,9432,0171,8103.82Individual error (%)0.504.005.202.850.859.50—Case 8Target2×1083×1085×1081,5002,0002,200—Preliminary values2.10×1083.08×1084.93×1081,5602,0922,4484.82Individual error (%)5.002.671.404.004.6011.27—Case 9Target2×10823.40×108a5×1082,0001,0212,000—Preliminary values2.01×10822.95×1085.07×1082,0241,0221,8432.16Individual error (%)0.501.921.401.200.107.85—Case 10Target2×10823.40×108a5×1081,9001,0212,100—Preliminary values2.52×10824.32×1086.53×1081,5491,0701,70117.13Individual error (%)26.003.9330.6018.474.8019.00—Primary Inversion Using a Fluid–Solid ModelAfter the preliminary inversion using the solid–only wave solver, the three frequency-level GA-based optimizer that uses a fluid–solid wave model estimated the fluid layer’s depth information and the material properties (E and ρ) of the soil layers. Namely, the preliminary inversion results from Cases 9 and 10 were used to determine the bounds of the guessed parameters in the first-frequency set of the primary GA inversion using the fluid–solid model. The increasing-frequency continuation approach was utilized with a dynamic pulse with a dominant frequency of 5 Hz in the first frequency-level GA and then 10 and 20 Hz were used in the next frequency levels; the values of the estimated parameters were bounded in each frequency set in the following manner: •During the first frequency-level GA using the fluid–solid wave solver, the estimated Ls1 and Lf values were bounded at 2 m≤Ls1≤7 m and 8 m≤Lf≤13 m, the estimated values of Es1 and Es2 were bounded at 1×108 Pa≤E≤7.5×108 Pa, and the estimated values of ρs1 and ρs2 (i.e., the mass densities of surrounding soil layers) were set to 1,500 kg/m3≤ρ≤2,500 kg/m3.•In the second frequency-level GA, the values of estimated Ls1 and Lf are bounded using the same upper and lower limits as the first frequency-level GA. The values of estimated Es1 and Es2 were bounded using ±50% derivations of their final-reconstructed values obtained during the first frequency–level GA, whereas the values of estimated ρs1 and ρs2 were bounded using ±20% derivations.•Finally, in the last frequency-level GA, the estimated Ls1 and Lf values were bounded using ±5% derivations of their final-reconstructed values obtained during the second frequency-level GA, ±10% derivations of the values of estimated Es1 and Es2 were employed to bound the limits of Young’s moduli, and ±10% derivations of the final-reconstructed values of ρs1 and ρs2 were used as theirs bounds in the last frequency-level GA.Table 4 presents the values of final reconstructed control parameters in Cases 9 and 10. For Case 9, the value of E¯ at the end of the last-level GA, 1.74%, was smaller than the error obtained in the preliminary result, 2.16% (Table 3). The same situation was found for Case 10—E¯ at the end of the last-level GA, 2.79%, was much smaller than E¯ from the preliminary inversion, 17.13%. Fig. 9 shows the average error for the best-fit individual over GA iterations of the three-level primary GA inversion using fluid–solid models for both Cases 9 and 10. During the initial several iterations of each frequency level, the algorithm focuses on exploring widely varying values of the estimated control parameters within the given limits, and this behavior increases error. Subsequently, the error of the algorithm is converged by finding better-fit individuals. Fig. 10 shows histograms of the estimated control parameters of the entire population of the individuals during all GA levels in Case 10 during the primary inversion. It demonstrates that, at the beginning of each frequency level, the GA explored wide ranges of control parameters, and eventually converged to the targeted parameters at the end of the last frequency set.Table 4. Example 3 results of primary inversion in Cases 9 and 10: estimated values of control parameters in each GA level obtained by increasing frequency of dynamic force for each GA levelTable 4. Example 3 results of primary inversion in Cases 9 and 10: estimated values of control parameters in each GA level obtained by increasing frequency of dynamic force for each GA levelCaseGA levelLs1 (m)Lf (m)Es1 (Pa)Es2 (Pa)ρs1 (kg/m3)ρs2 (kg/m3)Average error, E¯ (%)Case 9Target5102×1085×1082,0002,000—First5.110.02.01×1085.41×1081,9651,8433.38Second5.19.82.01×1085.45×1081,9691,8433.82Third126.96.36.199×1085.12×1082,0111,9201.74Individual error (%)2.001.000.502.400.554.00—Case 10Target49.32×1085×1081,9002,100—First188.8.131.52×1084.68×1082,4032,23915.36Second4.19.32.02×1085.42×1081,9221,9813.12Third4.09.32.01×1085.37×1081,9221,9392.79Individual error (%)0.000.000.507.401.167.67—Example 3 shows a proof of concept, potentially promising that (1) fluid–filled cavities in a 2D or 3D domain can be detected using the all-solid-element-based preliminary inversion of E, G, and ρ and finding particular solid elements (with a high value of E, a low value of G, and a value of ρ close to 1,021 kg/m3); and (2) a fluid–solid–element-based wave solver can be used to identify the geometry of a fluid–solid interface and E, G, and ρ of solid elements, surrounding fluid, in the 2D or 3D domain. The details of such 2D or 3D extensions are presented subsequently.DiscussionLimitationsThis paper employed a 1D model, which takes into account a horizontally layered geological model and cannot consider realistic three-dimensional random variation of material properties of geological formations. In this study, the authors did not employ field experimental validation because it was infeasible for the 1D wave model, which does not consider the geometrical damping and the vector wave behaviors in a 3D unbounded solid domain, to be validated through the field experiment.RecommendationsThis work can be extended into 2D and 3D settings in a manner such that one can identify the geometry of fluid–filled cavities and the material properties of surrounding soils by using a FWI algorithm (instead of a GA). The GA is a useful method for identifying a small number of control parameters in a feasibility study. However, when the number of control parameters is large (for example, in an inverse problem in which the material properties E, G, and ρ of all the solid elements in a 2D or 3D setting are to be estimated), the GA is less viable than FWI, which is powered by the PDE-constrained optimization method.Under such a robust FWI, the following two-step approach (i.e., preliminary to primary inversion) can be employed to (1) detect the existence of fluid–filled voids, and then (2) accurately identify the geometry of a fluid–solid interface and the material proprieties of soils in 2D or 3D settings. •Preliminary inversion: One can detect approximately the existence of fluid–filled cavities in a 2D or 3D domain by using the all-solid-element-based inversion of E, G, and ρ of solids. In each iteration of this preliminary FWI, the guessed profiles in all the elements (i.e., E, G, and ρ) are updated using a gradient vector. Then it searches an area of interest that has a high value of elasticity but a low value of shear modulus and a mass density that is similar to that of water. In turn, one can consider such an area as a fluid–filled cavity.•Primary inversion: After the preliminary FWI, one uses a wave solver using fluid elements coupled with solid elements. Then one updates the control parameters E, G, and ρ of the solid elements surrounding the potential fluid–filled cavity (modeled as fluid elements) detected from the preliminary inversion, as well as the geometry parameters of the interface between fluid and solid elements. The initial guessed profile of the interface is based on the result from the preliminary all-solid-element FWI. Fig. 11 depicts a possible scenario in which the initially guessed geometry of the fluid–solid interface is approximated based on the result from the preliminary all-solid FWI.The preceding proposed method can be used as both initial screening (i.e., a preliminary inversion using all solid elements) and in-depth cavity characterization (i.e., a primary inversion using both fluid and solid elements). Furthermore, engineers or geologists can somehow know of the existence of water-filled underground cavities or internal soil erosion activities because of surface signs (e.g., depression or subsidence) or geological or hydrogeological phenomena (e.g., sinkhole raveling activities) in most practical geoengineering application sites. Thus, they may have a combination of preliminary data from several different methods, including the proposed preliminary inversion, that will greatly benefit the proposed primary inversion using fluid and solid elements in most cases.ConclusionsThis work presents a new multilevel GA–based inversion method combined with a frequency-continuation scheme. Three levels of GA are used, and in each level, a pulse signal with a different dominant frequency is utilized to tackle the solution multiplicity of the presented inverse problem. Numerical results show the following findings about the proposed inversion method. First, it is feasible to estimate the depths of solid–fluid interfaces and Young’s moduli of the surrounding solid layers by employing an acoustic-elastodynamic wave model and a measured dynamic response at a sensor on the top surface. Second, combining the multilevel GA–based optimizer with a frequency-continuation scheme improves the performance of identifying targeted control parameters. Third, the frequency-continuation scheme in which the frequency is decreased progressively for each GA level is as effective as the conventional frequency-continuation scheme (i.e., increasing the frequency for each level) in the presented 1D setting. Lastly, a preliminary inversion using a solid-only model can address the uncertain presence of a fluid layer in a domain by determining that one of the layers has a very large value of Young’s modulus (with a value similar to that of the bulk modulus of water) and a mass density similar to that of water. After the preliminary inversion, a primary inversion—based on a multiple frequency–level GA using a fluid–solid model—can adjust the fluid layer’s depth information and the material properties of the soil layers. The authors also discussed how to extend the presented method in realistic 2D or 3D settings. Such an extension effectively may detect water-filled underground cavities and estimate the risk of potential urban cave-in and karst sinkholes. In typical engineering projects, engineers do not rely on a single geophysical method no matter how advanced its theories or computation capacities are. Thus, the presented method and its extension can be an important supplement to existing testing methods rather than a standalone solution in practice.Data Availability StatementSome or all data, models, or code generated or used during the study are available from the corresponding author by request, including MATLAB code (.m format) of the presented inverse modeling, and MATLAB datasets (.mat format) of the presented numerical results.AcknowledgmentsThis material is based upon work supported by the National Science Foundation under Awards CMMI-2044887 and CMMI-2053694. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. The authors are grateful for the support by the Faculty Research and Creative Endeavors (FRCE) Research Grant at Central Michigan University. This paper is a contribution of the Central Michigan University Institute for Great Lakes Research. The authors also greatly appreciate the reviewers’ very constructive comments, which substantially helped to improve the paper.References Beck, B. F. 1986. “A generalized genetic framework for the development of sinkholes and karst in Florida, U.S.A.” Environ. Geol. Water Sci. 8 (1): 5. https://doi.org/10.1007/BF02525554. Bishop, I., P. Styles, S. Emsley, and N. Ferguson. 1997. “The detection of cavities using the microgravity technique: Case histories from mining and karstic environments.” Geol. Soc. London, Eng. Geol. Special Publications 12 (1): 153–166. https://doi.org/10.1144/GSL.ENG.1997.012.01.13. Brown, L. T., D. M. Boore, and K. H. Stokoe II. 2002. “Comparison of shear-wave slowness profiles at 10 strong-motion sites from noninvasive SASW measurements and measurements made in boreholes.” Bull. Seismol. Soc. Am. 92 (8): 3116–3133. https://doi.org/10.1785/0120020030. Bunks, C., F. M. Saleck, S. Zaleski, and G. Chavent. 1995. “Multiscale seismic waveform inversion.” Geophysics 60 (5): 1457–1473. https://doi.org/10.1190/1.1443880. Butler, D. K. 1984. “Microgravimetric and gravity gradient techniques for detection of subsurface cavities.” Geophysics 49 (7): 1084–1096. https://doi.org/10.1190/1.1441723. Coşkun, N. 2012. “The effectiveness of electrical resistivity imaging in sinkhole investigations.” Int. J. Phys. Sci. 7 (15): 2398–2405. https://doi.org/10.5897/IJPS11.1063. Fathi, A., L. F. Kallivokas, and B. Poursartip. 2015a. “Full-waveform inversion in three-dimensional PML-truncated elastic media.” Comput. Methods Appl. Mech. Eng. 296 (Nov): 39–72. https://doi.org/10.1016/j.cma.2015.07.008. Fathi, A., B. Poursartip, and L. F. Kallivokas. 2015b. “Time-domain hybrid formulations for wave simulations in three-dimensional PML-truncated heterogeneous media.” Int. J. Numer. Methods Eng. 101 (3): 165–198. https://doi.org/10.1002/nme.4780. Fathi, A., B. Poursartip, K. H. Stokoe II, and L. F. Kallivokas. 2016. “Three-dimensional P- and S-wave velocity profiling of geotechnical sites using full-waveform inversion driven by field data.” Soil Dyn. Earthquake Eng. 87 (Aug): 63–81. https://doi.org/10.1016/j.soildyn.2016.04.010. Fehdi, C., I. Nouioua, D. Belfar, L. Djabri, and E. Salameh. 2014. “Detection of underground cavities by combining electrical resistivity imaging and ground penetrating radar surveys: A case study from Draa Douamis area (North East of Algeria).” In H2Karst research in limestone hydrogeology, 69–82. New York: Springer. Guidio, B., and C. Jeong. 2021. “On the feasibility of simultaneous identification of a material property of a Timoshenko beam and a moving vibration source.” Eng. Struct. 227 (Jan): 111346. https://doi.org/10.1016/j.engstruct.2020.111346. Jeong, C., S.-W. Na, and L. F. Kallivokas. 2009. “Near-surface localization and shape identification of a scatterer embedded in a halfplane using scalar waves.” J. Comput. Acoust. 17 (3): 277–308. https://doi.org/10.1142/S0218396X09003963. Jeong, C., A. C. S. Peixoto, A. Aquino, S. Lloyd, and S. Arhin. 2017. “Genetic algorithm–based acoustic-source inversion approach to detect multiple moving wave sources of an arbitrary number.” J. Comput. Civ. Eng. 31 (5): 04017020. https://doi.org/10.1061/(ASCE)CP.1943-5487.0000664. Jung, J., C. Jeong, and E. Taciroglu. 2013. “Identification of a scatterer embedded in elastic heterogeneous media using dynamic XFEM.” Comput. Methods Appl. Mech. Eng. 259 (Jun): 50–63. https://doi.org/10.1016/j.cma.2013.03.001. Kallivokas, L., A. Fathi, S. Kucukcoban, K. Stokoe II, J. Bielak, and O. Ghattas. 2013. “Site characterization using full waveform inversion.” Soil Dyn. Earthquake Eng. 47 (Apr): 62–82. https://doi.org/10.1016/j.soildyn.2012.12.012. Kang, J. W., and L. F. Kallivokas. 2011. “The inverse medium problem in heterogeneous PML-truncated domains using scalar probing waves.” Comput. Methods Appl. Mech. Eng. 200 (1): 265–283. https://doi.org/10.1016/j.cma.2010.08.010. Kiflu, H., S. Kruse, M. Loke, P. Wilkinson, and D. Harro. 2016. “Improving resistivity survey resolution at sites with limited spatial extent using buried electrode arrays.” J. Appl. Geophys. 135 (Dec): 338–355. https://doi.org/10.1016/j.jappgeo.2016.10.011. Kucukcoban, S., H. Goh, and L. F. Kallivokas. 2019. “On the full-waveform inversion of Lamé parameters in semi-infinite solids in plane strain.” Int. J. Solids Struct. 164 (Jun): 104–119. https://doi.org/10.1016/j.ijsolstr.2019.01.019. Lloyd, S. F., C. Jeong, and M. Okutsu. 2016. “Numerical investigation on the feasibility of estimating the thickness of Europa’s ice shell by a planned impact.” J. Aerosp. Eng. 29 (5): 04016035. https://doi.org/10.1061/(ASCE)AS.1943-5525.0000621. Mirzanejad, M., K. T. Tran, M. McVay, D. Horhota, and S. J. Wasman. 2020. “Sinkhole detection with 3D full seismic waveform tomography.” Geophysics 85 (5): B169–B179. https://doi.org/10.1190/geo2019-0490.1. Nam, B. H., R. Shamet, M. Soliman, D. Wang, and H.-B. Yun. 2018. Development of a sinkhole risk evaluation program. Final Rep. No. BDV24-977-17. Tallahassee, FL: Florida DOT. Pakravan, A., J. W. Kang, and C. M. Newtson. 2016. “A Gauss–Newton full-waveform inversion for material profile reconstruction in viscoelastic semi-infinite solid media.” Inverse Prob. Sci. Eng. 24 (3): 393–421. https://doi.org/10.1080/17415977.2015.1046861. Petryk, H., and Z. Mroz. 1986. “Time derivatives of integrals and functionals defined on varying volume and surface domains.” Arch. Mech. 38 (5–6): 697–724. Tihansky, A. B. 1999. “Sinkholes, west-central Florida.” In Vol. 1182 of Land subsidence in the United States: US geological survey circular, 121–140. Reston, VA: USGS. Tran, K. T., and M. McVay. 2012. “Site characterization using Gauss–Newton inversion of 2-D full seismic waveform in the time domain.” Soil Dyn. Earthquake Eng. 43 (Dec): 16–24. https://doi.org/10.1016/j.soildyn.2012.07.004. Tran, K. T., M. McVay, M. Faraone, and D. Horhota. 2013. “Sinkhole detection using 2D full seismic waveform tomography.” Geophysics 78 (5): R175–R183. https://doi.org/10.1190/geo2013-0063.1. Xiao, H., Y. J. Kim, B. H. Nam, and D. Wang. 2016. “Investigation of the impacts of local-scale hydrogeologic conditions on sinkhole occurrence in East-Central Florida, USA.” Environ. Earth Sci. 75 (18): 1274. https://doi.org/10.1007/s12665-016-6086-3.