Inverse Scattering

Inverse Scattering

Fig11-TSWaveform tomography and in particular inverse scattering are essential to a broad spectrum of science and technology disciplines including geophysics, oceanography, medical imaging, energy production, and non-destructive material testing. In general the relationship between the wavefield scattered by an anomaly and its geometry (or physical characteristics) is nonlinear, which invites two overt solution strategies: i) linearization via e.g. Born approximation, or ii) pursuit of the nonlinear minimization approach. Over the past two decades, however, a number of sampling methods have emerged that both consider the nonlinear nature of the inverse problem and dispense with iterations. Commonly, these techniques deploy an indicator functional that varies with coordinates of the trial i.e. sampling point, and projects observations of the scattered field onto a functional space reflecting the baseline wave motion in a reference (“anomaly-free”) medium. Such indicator functional, designed to reach extreme values when the sampling point strikes the anomaly, accordingly provides a tomogram via its (thresholded) spatial distribution. Examples of such imaging paradigm include the Linear Sampling Method and the Method of Topological Sensitivity described below. In situations involving distinct anomalies, a tractable solution to the waveform tomography problem can also be sought via the boundary integral equation (BIE) framework and the concept of shape and material sensitivities. Our work also entails studies of the interior transmission problem which, despite its central role in the theory of inverse scattering, is not covered by the existing theories in partial differential equations.

Method of Topological Sensitivity

Since its birth in the context of shape optimization, the idea of Topological Sensitivity (TS) has been generalized and applied to tackle a variety of inverse scattering problems. Our research in this area is largely focused on the inverse problems in acoustics [27345361] and elastodynamics [171922293335374851565758]. In the imaging approach the TS, which quantifies the perturbation of a given cost functional due to nucleation of an infinitesimal defect in the (reference) background medium, is used as an effective anomaly indicator through an assembly of sampling points where it attains extreme negative values. Typically, the TS formulas are amenable to explicit representation in terms of the wavefields computed for the reference domain, which is the root of computational efficiency of this class of waveform tomography solutions. Our recent work has involved both experimental verification [56] and theoretical justification [61] of the TS approach to inverse scattering.

TS reconstruction of a crescent-shaped anomaly from sparse cross-hole seismic measurements [35]. Arranged clockwise: images obtained using low, intermediate, high, and multi- frequency illumination.

Linear Sampling Method

The Linear Sampling Method is a non-iterative technique for reconstructing obstacles from the scattered field data. This method entails the solution of a linear Fredholm equation of the first kind. Here, the linearity of the equation does not stem from an approximation of the physical conditions (e.g. the weak scatterer hypothesis); rather, it is the result of an equivalence between the nonlinear inverse scattering problem and the linear integral equation. Theory shows that the norm of a solution to the germane (Fredholm) equation becomes unbounded for sampling points outside of the support of a scatterer, which de facto yields a characteristic function of the anomaly. Notable advantages of this approach to inverse scattering include computational efficiency and insensitivity to the boundary conditions present on the surface of a scatterer. Our work in this area [182631324253] is focused on extending the Linear Sampling Method to inverse problems involving elastic waves, near-field sensory data, heterogeneous reference domains, variety of anomalies (inclusions, fractures, interfaces), and more recently unknown background media.

Three-dimensional LSM reconstruction of a tunnel in a layered half-space [31] from surface seismic measurements.

Obstacle sensing by shape-material sensitivity framework

This facet of our research [151636] focuses on the elastic-wave reconstruction and characterization of discrete heterogeneities (cavities, inclusions) in an otherwise homogeneous solid from the limited-aperture waveform observations taken on its surface. On adopting the boundary integral equation (BIE) framework as basis for describing the forward scattering problem, the inverse query is cast as a minimization task involving sensory data and their simulations for a trial anomaly that is defined through its boundary, elastic moduli, and mass density. For an optimal performance of the gradient-based search methods suited to solve the problem, explicit expressions for the shape (i.e. boundary) and material sensitivities of the misfit functional are obtained via an adjoint field approach and direct differentiation of the governing BIEs. Making use of the message-passing interface, the featured sensitivity formulas are implemented in a data-parallel code and integrated into a nonlinear optimization framework based on the direct BIE method and an augmented Lagrangian — whose inequality constraints are employed to avoid solving forward scattering problems for physically inadmissible (or overly distorted) trial inclusion configurations. Numerical results for the reconstruction of an ellipsoidal defect in a semi-infinite solid show the effectiveness of the proposed shape-material sensitivity formulation, which constitutes an essential component of the defect identification algorithm.


Shape sensitivity-driven evolution of a trial subterranean void in the boundary-only imaging process [15]: iterations k=0 and k=26 denote respectively the initial “guess” and final shape, the latter overlapping with the “true” defect. The ground i.e. testing surface is indicated in dark blue.

Interior transmission problem

This facet of our research deals with the theory of the interior transmission problem (ITP) for heterogeneous and anisotropic (visco-) elastic solids [4050]. The subject is central to the so-called qualitative methods (e.g. Linear Sampling Method) for inverse scattering involving penetrable obstacles. Although simply stated as a coupled pair of elastodynamic wave equations, the ITP for elastic bodies is neither self-adjoint nor elliptic. The aim of this work is to provide a systematic treatment of the ITP for elastic bodies, considering a broad range of material-contrast configurations. In particular, we investigate questions of the existence and uniqueness of a solution to the ITP, the discreteness of its eigenvalues, and the existence of such eigenvalue spectrum. Necessitated by the breadth of material configurations studied, the relevant claims are established via a suite of variational formulations, each customized to meet the needs of a particular subclass of eigenvalue problems. For instance the analysis shows that in the case of a viscoelastic scatterer embedded in an elastic solid, the ITP eigenvalue spectrum is empty.

Elastodynamic ITP configuration [40]: scatterer composed of inclusions D∗ and cavities Do (left) and scatterer support, D, occupied by the background material (right)

Elastodynamic ITP configuration [40]: scatterer composed of inclusions D∗ and cavities Do (left) and scatterer support, D, occupied by the background material (right)

Experimental Verification

We are particularly interested in validating our theories by way of carefully designed laboratory experiments. In this vein the method of topological sensitivity (TS) for solving inverse scattering problems, previously supported by a multitude of numerical simulations, is put to test experimentally [56] with the focus on 2D obstacle reconstruction in a thin aluminum plate using elastic waves. To this end, in-plane measurements of transient elastodynamic waveforms along the edges of the plate are captured in a non-contact fashion by a 3D Scanning Laser Doppler Vibrometer (SLDV). Using an elastodynamic (time-domain) finite element model as a computational platform, the TS reconstruction maps are obtained and analyzed under varying experimental conditions. The results show significant agreement between the defect geometry and its reconstruction, thus demonstrating the utility of the TS approach as an efficient and robust solution tool for this class of inverse problems. For completeness, the experimental investigation includes a set of parametric studies geared toward exposing the effect of key problem parameters on the quality of obstacle reconstruction such as the (dominant) excitation frequency, the source aperture, and duration of the temporal records. On the analytical front, it is shown that the use of a suitable temporal windowing function in specifying the L2 cost functional (that underpins the TS formulation) is essential from both theoretical and computational points of view.

Elastodynamic scattering experiment [56] featuring a 3D Scanning Laser Doppler Vibrometer (top), and experimental TS maps of an aluminum plate containing a hole and a slit (bottom).

Nonlinear Waves


Nonlinear Waves in Soft Solids

We aim to understand the fundamentals of weakly nonlinear wave motion in soft, tissue-like solids due to action of a focused, high-intensity ultrasound. We are interested in both (i) developing the asymptotic models of nonlinear ultrasound and (ii) exposing its mean effects, synthesized via the phenomenon known as acoustic wind” or acoustic radiation force. Our research has lead to the development of a new experimental technique for quantifying locally the nonlinear elasticity in tissue-like solids, that permits measurement over a mm-sized support of the acoustic radiation force — supplied by a focused ultrasound beam.

Acoustic Radiation Force (ARF) in soft tissues

The focus of our study is a sustained body force (the so-called acoustic radiation force) in homogeneous tissue-like solids [4349] generated by an elevated-intensity, focused ultrasound field (Mach # = O(10^(-3))) in situations when the latter is modulated by a low-frequency signal. This intermediate-asymptotics problem, which bears relevance to a range of emerging biomedical applications, is characterized by a number of small (but non-vanishing) parameters including the Mach number, the ratio between the modulation and ultrasound frequency, the ratio of the shear to bulk modulus, and the dimensionless ultrasound attenuation. On approximating the response of soft tissues as that of a nonlinear viscoelastic solid with heat conduction, the featured second-order problem is tackled via a scaling paradigm wherein the transverse coordinates are scaled by the width of the focal region, while the axial and temporal coordinate are each split into a “fast” and “slow” component with the twin aim of: (i) canceling the linear terms from nonlinear field equations, and (ii) precisely tracking the effect of ultrasound modulation. In this way the acoustic radiation force (ARF), giving rise to the mean tissue motion, is exacted by computing the “fast” time average of the germane field equations. A comparison with the existing theory reveals a number of key features that are brought to light by the new formulation, including the contributions to the ARF of ultrasound modulation and thermal expansion, as well as the precise role of constitutive nonlinearities in generating the sustained body force in tissue-like solids by a focused ultrasound beam. Motivated by this result, we further investigate the possibility of estimating locally the mixed (deviatoric-volumetric) modulus of nonlinear tissue elasticity [60] from the amplitude of ARF-generated shear waves. Our prototype technique is verified experimentally by comparisons with an established methodology known as acoustoelasticity.

Experiment versus simulation: shear waves in a paraffin-based phantom (G=7kPa) due to ARF generated by a 4 MHz Verasonics linear array [60].

Experiment versus simulation: shear waves in a paraffin-based phantom (G=7kPa) due to ARF generated by a 4 MHz Verasonics linear array [60].

Dual-time simulation of nonlinear ultrasound fields

This research introduces a KZK-type model [52] for describing nonlinear wave propagation due to “slowly” modulated ultrasound excitation, where the ratio between the modulation frequency and carrier frequency is o(1). The problem is relevant to a variety of diagnostic techniques where the acoustic radiation force, generated by a focused ultrasound beam, is used to palpate tissue locally within the focal region. An indispensable prerequisite for computing the ARF however, is a precise knowledge of the nonlinear wave field generated by the ultrasound transducer at hand. In this class of applications, the MHz carrier signal is commonly modulated by frequencies on the order of 10^2÷10^4 Hz which, in the context of time-domain numerical algorithms, dramatically increases the computational effort necessary to simulate the nonlinear wavefield. To alleviate the impediment, our approach makes use of a dual time-scale approach towards the development of an intermediate  (KZK-type) asymptotics solution that fits the gap between the quasi-static approximation and the full solution. For completeness, we also develop a matching computational platform [54] for solving the nonlinear parabolic wave equation. The key advantage of the latter (dual-time) algorithm lies in the fact that the nonlinear solution is approximated, for an arbitrary modulation envelope, using a set of precomputed steady-state solutions. This makes the method computationally inexpensive, and allows for the fast simulation of long excitation signals for which the time-domain simulations are inherently fruitless. Numerical results demonstrate that the accuracy of our approximation is proportional to the ratio between the modulation and carrier (i.e. ultrasound) frequency, as suggested by theoretical predictions.

Schematics of a modulated ultrasound beam [54] generated by phase-focused transducer array (left) and convex-shaped transducer (right).

Waves in Complex Media

Fig31-chessWave Motion in Periodic & Random Media

We are interested in two fundamental questions: (i) how does the presence of a “microstructure”, whose characteristic length scale is decade(s) smaller than the germane wavelength, affect the mean wave motion, and (ii) what knowledge about the “microstructure” of an anomaly can be gained remotely, by considering the inverse scattering problem with long-wavelength illumination. 

Second-order homogenization of waves in periodic solids

The goal of this research is to better understand the mathematical structure and ramifications of the second-order homogenization of wave motion in periodic media. To this end, multiple-scales asymptotic approach is applied to the scalar wave equation in one and two spatial dimensions [59]. In contrast to previous studies where the second-order homogenization has lead to the introduction of a single fourth-order derivative in the field equation, our investigation demonstrates that such (asymptotic) approach results in a family of PDEs uniting spatial, temporal, and mixed fourth-order derivatives – that jointly control the incipient wave dispersion. Given the consequent freedom in weighing the homogenization parameters, the notion of an optimal asymptotic model is next considered in a 1D setting via its ability to capture the salient features of wave propagation within the first Brillouin zone, including the onset and magnitude of the phononic band gap. In the context of 2D wave motion, on the other hand, the asymptotic analysis is first established in a general setting, exposing the constant shear modulus as sufficient condition under which the second-order approximation of a bi-periodic elastic solid is both isotropic and limited to even-order derivatives. On adopting a chessboard-like periodic structure (with contrasts in both modulus and mass density) as a testbed for in-depth analytical treatment, we show that the second-order approximation of the mean wave motion is governed by a family fourth-order PDEs that: (i) entail exclusively even-order derivatives; (ii) describe anisotropic wave dispersion via the “sin4θ + cos4θ” term, and (iii) include the asymptotic model for a square lattice of circular inclusions as degenerate case. In this vein, we further consider the phenomenological models of gradient elasticity (GE) in that we investigate (a) how are the GE coefficients linked to the characteristics of a microstructure, and (b) whether the GE models of can be used to describe the presence of the first band gap [55].


Anisotropic wave dispersion in a chessboard-like solid [59]

Scattering by obstacles with periodic microstructure

Typically, higher-order homogenization analyses of waves in periodic media postulate unbounded domains in that they do not consider the asymptotic contribution of domain boundaries. From the mathematical investigations of Neumann and Dirichlet problems, however, it is known that such contribution is not only relevant but also difficult to quantify. To help bridge the gap, we study the homogenization of a transmission problem [63] arising in the scattering theory for bounded inhomogeneities with periodic coefficients, modeled by the anisotropic Helmholtz equation. The coefficients are assumed to be periodic functions of the fast variable, specified over the unit cell with characteristic size ε. By way of multiple scales expansion, we focus on the O(εk), k=1,2 bulk and boundary corrections of the leading-order (O(1)) homogenized transmission problem. The analysis in particular provides the H1 and L2 estimates of the error committed by the first-order- corrected solution considering i) bulk correction only, and ii) boundary and bulk correction. We treat explicitly the O(ε) boundary correction for the transmission problem when the scatterer is a unit square, and we show that is L2-limit as ε→0 exists provided that the boundary cutoff of cells is fixed. We further establish the O(ε2)-bulk correction of the transmission problem, and we distill its effect on the mean wavefield inside the scatterer. The analysis also highlights a previously established, yet scarcely recognized fact that the O(ε) bulk correction of the mean wavefield vanishes identically.

Scattering by an obstacle with periodic coefficients [63]

Wave dispersion by fractal microscopic structures

Typically, the ability to sense microstructure is determined by the ratio of scatterer size to probing wavelength. Here, we address the question of whether macroscopic waves can report back the presence and distribution of microscopic scatterers despite several decades difference in scale between the wavelength and scatterer size [62]. In our analysis, monosized hard scatterers 5  μm in radius are immersed in lossless gelatin phantoms to investigate the effect of multiple reflections on the propagation of shear waves with millimeter wavelength. Steady-state monochromatic waves are imaged in situ via magnetic resonance imaging, enabling quantification of the phase velocity at a voxel size large enough to contain thousands of individual scatterers, but small enough to resolve the probing wavelength. We show in theory, experiments, and simulations that the resulting coherent superposition of multiple reflections gives rise to power-law dispersion at the macroscopic scale if the scatterer distribution exhibits apparent fractality over an effective length scale that is comparable to the probing wavelength. Since apparent fractality is naturally present in any random medium, microstructure can thereby leave its fingerprint on the macroscopically quantifiable power-law exponent. Our results are generic to wave phenomena and carry great potential for sensing microstructure that exhibits intrinsic fractality, such as, for instance, vasculature.


Shear-wave dispersion experiment. (a) Confocal microscopy image of a gel sample with mono-sized glass spheres. (b) Phase of the shear-wave displacement as measured by MRE within a single slice of the gel sample. The wavelength λ (∼3 mm) is comparable to the effective correlation length whereby the apparent fractality can be probed. The third image shows the derived phase velocity cp from which mean values of cp are calculated. The latter are acquired at several vibration frequencies, yielding the apparent power-law dispersion (4th image).

Applied Inverse Problems


Applied Inverse Problems in NDE, Geophysics & Medicine

We aim to apply our expertise in solid mechanics, wave propagation, and inverse problems to a variety of engineering and technological problems where non-invasive diagnosis of subterranean regions, materials, structures, or living tissues is required.

3D seismic sensing of heterogeneous fractures

The goal of this research is to establish a comprehensive platform for the 3D reconstruction and interfacial characterization of arbitrarily-shaped, partially closed fractures by way of seismic or ultrasonic waves. In the ongoing work [64] that is relevant to both mine safety and hydraulic fracturing , this goal is made possible via an extension to the so-called Generalized Linear Sampling Method (GLSM) to enable accurate geometric reconstruction of subsurface fractures regardless of their interfacial condition. Aided by this result, our inverse solution deploys a 3-step hybrid approach where: 1) the fracture surface is reconstructed without the knowledge of its contact condition; 2) given the geometry, the fracture opening displacement (FOD) profile is recovered from an integral transform relating FOD to the scattered seismic waveforms; and 3) given FOD, the spatial distribution of specific stiffnesses is resolved from the boundary integral equation for a fracture, incorporating its (inhomogeneous) elastic contact condition. The proposed developments are being verified in a laboratory setting, making use of the recently acquired Scanning Laser Doppler Vibrometer (SLDV) as a motion sensing tool. In parallel, the concept of topological sensitivity (TS) is extended [58] to enable simultaneous 3D reconstruction of partially closed fractures and qualitative characterization of their interfacial condition by way of elastic waves. Interactions between the two surfaces of a fracture, due to e.g. presence of asperities, fluid, or proppant, are described via the Schoenberg’s linear slip model. The proposed TS sensing platform is formulated in the frequency domain, and entails point-wise interrogation of the subsurface by infinitesimal fissures endowed with interfacial stiffness. For completeness, the featured elastic polarization tensor – central to the TS formula – is mathematically described in terms of the shear and normal specific stiffnesses sn) of a vanishing fracture. Simulations demonstrate that, irrespective of the contact condition between the faces of a hidden fracture, the TS is capable of reconstructing its geometry and identifying the normal vector to the fracture surface without iterations. On the basis of such geometrical information, it is further shown via asymptotic analysis – assuming “low frequency” illumination, that by certain choices of (κsn) characterizing the trial (infinitesimal) fracture, the ratio between the shear and normal specific stiffness along the surface of nearly-planar (finite) fractures can be qualitatively identified.

Elastodynamic reconstruction of a partially-closed cylindrical fracture by the Linear Sampling Method [64]

Viscoelastic characterization of soft tissues

In this suite of works, we develop customized solutions for the vibration- i.e. wave-based identification of viscoelastic properties of soft tissues [47] with an emphasis on skin [3844].  One example is a theoretical and computational platform behind the operation of a piezoelectric sensor array for the surface monitoring of tissue motion that can be used as a basis for the reconstruction of layered viscoelastic skin properties [44]. This is accomplished by the scale reduction of the Multi-channel Analysis of Surface Waves (MASW), a methodology that is successfully used in engineering geophysics for the seismic-wave reconstruction of vertical geological profiles. The utility of the new sensor, containing an array of hair-thin PVDF sensors that are sensitive to surficial tissue motion, is enhanced through a comprehensive skin-fiber interaction analysis that furnishes integral information, cumulative over the length of a fiber, about the attenuation and dispersion of surface waves. Having such predictive model as a lynchpin of the back-analysis of electric charges furnished by the fibers, the methodology allows for an effective reconstruction and viscoelastic characterization of cutaneous and subcutaneous tissue sublayers on a millimeter scale. Another example of research in this direction revolves around the use of modulated ARF as a vibration source [38], which enables (via the analysis of Scholte waves) viscoelastic characterization of skin layers less than a millimeter in thickness.

Testing configuration [38] for the viscoelastic characterization of mm-thin tissue mimicking phantoms.

Probing nonlinear tissue elasticity by ARF

Prompted by a recent finding that the magnitude of the acoustic radiation force (ARF) in isotropic tissue-like solids depends linearly on a particular third-order elastic modulus — hereon denoted by C,this study [60] investigates the possibility of estimating C from the amplitude of ARF-generated shear waves. The featured coefficient of nonlinear elasticity, which captures the incipient nonlinear interaction between the volumetric and deviatoric modes of deformation, has so far received only a limited attention in the context of soft tissues due to the fact that the latter are often approximated as (i) fluid-like when considering ultrasound waves, and (ii) incompressible under static deformations. On establishing the analytical and computational platform for the proposed sensing methodology, the study proceeds with applying the prototype technique toward estimating via ARF the third-order modulus C in a series of tissue-mimicking phantoms. To help validate the concept and its implementation, the germane third-order modulus is independently estimated in each phantom via an established technique known as acoustoelasticity. The C-estimates obtained respectively via acoustoelasticity and the new theory of ARF show a significant degree of consistency. The key features of the new sensing methodology are that: (a) it requires no external deformation of an organ other than that produced by the ARF, and (b) it estimates the nonlinear C-modulus locally, over the (mm-sized) focal region of an ultrasound beam—where the shear waves are being generated.

Schematics (left) and picture (right) of the experimental setup [60] for measuring the third-order moduli C and D in tissue-mimicking phantoms by way of acoustoelasticity.

MR elastography and Vibro-acoustography

This research deploys the concept of Topological Sensitivity (TS) toward boosting the performance of existing techniques for tissue elasticity imaging such as Magnetic Resonance Elastography (MRE) and Vibro-acoustography [485157]. The work in [57], for instance, investigates the application of TS toward reconstructing and identifying tissue anomalies from MRE measurements. The basic idea behind MRE imaging is to apply time–harmonic vibration to a tissue inside the magnetic resonance (MR) scanner, and to capture thus induced 3D wave motion via suitable sequencing of the phase-coordinated, freeze-frame MR scans. To account for the facts that (i) the displacement of a material point captured by the MRE signifies the volume average over a reference voxel size, and (ii) any given measurement voxel may contain, at least partially, the lesion of interest, the concept of TS is extended to allow for (1) averaged volumetric measurements and (2) overlapping interrogation and measurement domains. To circumvent the difficulties involved in numerical modeling of an entire organ, the latter is subdivided into an array of cubic subdomains, each consisting of N3 voxels and playing the role of a reference body for TS imaging. To accomplish the sought partitioning, triaxial displacements along the boundary of each subdomain (captured by the MR scanner) are deployed as Dirichlet data — thus specifying a separate boundary value problem, whereas the internal measurements within the subdomain are taken as the observations for lesion reconstruction. The performance of our methodology for lesion reconstruction and tissue characterization is examined via both numerical examples and a preliminary application to in vivo MRE data. The results highlight the potential of the TS method toward elevating the resolution of tissue (visco-) elasticity reconstruction by MRE.


MRE experimental setup [57] (adapted from Sinkus, 2005). The breast tissue is vibrated by a shaker acting on its side; both the patient and the shaker are placed inside the MR scanner

Nondestructive evaluation of roads

In this work, we apply elastodynamic and electromagnetic model-based inversion toward an enhanced diagnosis of pavement systems via Falling Weight Deflectometer (FWD), Light Weight Deflectometer (LWD), and Ground Penetrating Radar (GPR) [1320304546]. The FWD test is one of the most commonly used tools for nondestructive evaluation of flexible pavements. Although the test is intrinsically dynamic, the state-of-practice backcalculation techniques used to interpret the FWD data are primarily based on elastostatic solutions due to high computational cost of dynamic multilayered analyses. It has long been known that the foregoing discrepancy may lead to systematic errors in the estimation of pavement moduli, especially in the situations of pronounced inertial phenomena due to shallow bedrock or seasonal stiff layer. In this investigation a simple, yet effective algorithm is proposed that allows the static backcalculation analyses to perform well even when dynamic effects are significant. The technique is based on the use of discrete Fourier transform as a preprocessing tool to filter the dynamic effects and extract the static pavement response from transient FWD records. Using the filtered (i.e. zero-frequency) force and deflection values in lieu of their peak counterparts, the static backcalculation can be further performed in a conventional manner, but free of inconsistencies associated with the neglect of inertial effects. Our results based on synthetic deflection records demonstrate a marked improvement in the elastostatic prediction of pavement moduli when the proposed modification is used. This approach has been implemented, together with a state-of-the-art waveform analysis of GPR data [46], into a computer platform GopherCalc for high-fidelity diagnosis of pavement systems.


GopherCalc graphical user interface [30] for the interpretation of FWD and GPR data.

Dynamic site characterization

The goal of this research [1102339] is to develop advanced inverse solutions for dynamic site characterization via e.g. Spectral Analysis of Surface Waves that are free of the customary modeling simplifications. For example in [23] we investigate shallow seismic profiling via the dispersion and attenuation analysis of Love surface waves. Although the Rayleigh waves are now commonly used as an engineering tool for nonintrusive identification of shallow subsurface stratigraphy, little attention has so far been paid to utilizing their horizontally polarized, Love-wave counterpart. The proposed methodology revolves around a causal viscoelastic model for the wave motion in a layered half-space caused by action of a surficial, torsionally vibrating disc. It is shown that the use of Love waves as a sounding tool reduces the number of material parameters relevant to viscoelastic site characterization by precluding the effect of the Poisson’s ratio and the P-wave quality factor in each geologic layer. For practical purposes, solution to the inverse problem is reduced to the minimization of a spectral misfit between experimental observations and theoretical predictions of the horizontally polarized surface ground motion. To maintain the rigor and robustness of the inverse solution, sensitivities of the forward model with respect to layer parameters are computed analytically, while the germane cost functional is established within the framework of the maximum-likelihood inverse theory. The Love-wave testing methodology for which this article establishes a necessary theoretical framework may be used either in a stand-alone fashion, or as a complement to the spectral analysis of Rayleigh waves. The proposed technique may be especially useful for (i) shallow surveys where material dissipation is of interest, and (ii) accurate characterization of saturated layers in which vertically polarized waves find limited use owing to their sensitivity to the pore fluid.

Experimental setup for the Spectral Analysis of Love Waves [23]

Experimental setup for the Spectral Analysis of Love Waves [23]




We aim to develop a comprehensive analytical and computational platform toward better understanding of the wave motion in heterogeneous semi-infinite solids. Our research applies primarily to engineering-scale geodynamics — relevant e.g. to the seismic response of objects of civil infrastructure

Viscoelastic waves in layered media

In this research, we establish a rigorous treatment of the singular visco-elastodynamic solutions for a semi-infinite multi-layered solid [8]. We show explicitly via an asymptotic analysis of the propagator matrices that the singular components of the dynamic Green’s functions, which are critical to the theoretical foundation of boundary integral equation methods, correspond fully to the static point-load solutions for an appropriate bi-material full-space. With the aid of the analytical expressions for the bi-material response [6], a computational formulation [11] for the multi-layered Green’s functions is also developed where the integral representation of the solution is decomposed into a closed-form singular part and a residual component which is amenable to numerical contour integration. Thanks to the foregoing treatment, the multi-layered fundamental solutions can be accurately and efficiently evaluated for a wide range of material and geometric configurations, including the special cases of elastic strata and the source points at the interface between two layers. As an illustration, we demonstrate the performance of the method in simulating the exact solution for a functionally-graded half-space with a linear wave velocity profile [8].


Singular nature of the elastodynamic Green’s function urr for a multi-layered half-space: comparison with the Mindlin’s solution and elastostatic bi-material Green’s function.

Wave motion in functionally-graded solids

As a fundamental step toward understanding the mechanics of wave motion in functionally-graded solids, we study 3D elastic wave propagation in a vertically-heterogeneous half-space with a linear shear wave velocity profile [3]. With the aid of a customized displacement-potential representation, Hankel transforms, and Fourier decompositions, the dynamic response of the semi-infinite solid to an arbitrarily distributed buried source is shown to admit integral representation in terms of modified Bessel functions whose order is either real (at lower frequencies) or imaginary (at higher frequencies). Specific aspects of the problem such as multiple poles along the inversion path in the complex plane and the characteristics of wave propagation in the vertical direction are elucidated. To cater for the computational treatment of elastodynamic boundary value problems in functionally graded solids, the solution was degenerated to a set of ring-load and point-load Green’s functions which are fundamental to boundary integral equation formulations [4]. A canonical solution for a disk vibrating vertically on the surface of a linear-wave-velocity half-space [5] is derived to highlight salient features of the inhomogeneous half-space problem, including the presence of local maxima in the dynamic foundation compliance function and the existence of cutoff frequencies below which the energy generated by the foundation does not propagate toward infinity.

Wave motion in a functionally-graded half-space whose moduli increase quadratically with depth [4].

Dynamic soil-structure interaction

The general subject of dynamic soil-structure interaction is relevant to a wide range of applications such as earthquake-resistant design of structures, foundations, tunnels and lifelines. One common challenge to these problems is the need to deal with complex stress wave propagation patterns in a supporting geological medium. To help address the problem, we developed a regularized format of the direct boundary integral equation (BIE) method [7] for three-dimensional elastodynamics in semi-infinite solids. Founded on the premise of a complete decomposition of the Greens functions into their singular and regular parts, the regularized BIE formulation is shown to be analytically appealing without demanding mathematical and numerical complexities such as Cauchy principal values. Extended to deal with general seismic soil-structure interaction problems in semi-infinite media, the formulation is implemented computationally together with a robust singularity treatment of multi-layered viscoelastic half-space Greens functions [811] as well as a rigorous account of the radiation condition for semi-infinite solids [24] and singular contact tractions [25] in the related mixed boundary value problems. On the component level, we have also investigated the dynamics of surface foundations resting on a semi-infinite solid [5912] via indirect BIE formulations.

Dynamic contact tractions beneath a vibrating block foundation [7].

Thin Films


Thermal Cracking of Thin Films

There is an abundance of natural and manufactured quasi-brittle systems whose one dimension is decade(s) smaller than the other two — commonly referred to as thin films. We aim to better understand the mechanisms of distributed and diffuse failure in such systems due to straining caused by e.g. cooling, drying, or tectonic deformation.

Crack spacing in strained coatings & strata

When the thermally induced stress in a shrinking pavement layer reaches the tensile strength of asphalt, regularly spaced thermal cracks form across the width of the pavement. To better understand this phenomenon, a one-dimensional analytical solution for the stress distribution in a thermally shrinking elastic layer placed on an elastoplastic, cohesive–frictional base is developed and validated by comparison with a 2D numerical solution [1421]. From the analytical model, a prediction of a length scale that provides bounds on the crack spacing is obtained. Predicted bounds on crack spacing are validated by comparison with field observations in asphalt pavements. It is demonstrated that the proposed formulation can also be applied to estimate the average crack density observed in thin ceramic films subjected to the application of an axial strain; in the latter system, the crack spacing is six decades smaller than that observed in pavement systems. Similarly, when a layered system of rock beds is subjected to a sufficiently large extensional strain, joints form in the competent layers. Previous analyses have shown that the ratio between joint spacing and competent layer thickness decreases as the applied strain increases. Further, if the entire interface between the competent layer and the matrix fails in shear (slip), no new joints can form and a lower bound on the joint spacing is reached. By extending the work in [14], a joint spacing analysis is developed to explicitly account for the effects of overburden depth [28]. The resulting model is a pair of non-linear equations that can be solved for the characteristic joint spacing as a function of layer thickness, applied strain, and overburden depth. The model results show that, for a given applied strain, the joint spacing first decreases and then increases with depth. This behavior is controlled by the opposing effects of depth-increasing shear strength along the competent layer–matrix interface and depth-increasing compressive prestress. The analysis also reveals that the lower bound (saturation) joint spacing is strongly dependent on depth. Within the choice of realistic physical parameters, predicted values of the saturation-spacing-to-layer-thickness ratio span the range of values observed in the field.

Joint spacing in layered rocks [28]: Schematics of the layered bed system (top), representative domain (bottom left), and (c) approximate shear stress distribution at a section x=const. in the competent layer. Dark green is the competent layer, and light blue indicates the matrix.

Smeared crack formulations

The smeared-crack approach of strain-softening is a popular engineering concept for the numerical simulation of tensile failure processes in quasi-brittle materials such as rock and concrete. Despite numerous publications on the equivalent constitutive formulation of cracked solids, little is known about the predictive value of the smeared-crack approach, particularly when non-proportional load histories are considered. To this end, two failure indicators are examined at the constitutive level in order to assess the predicted failure mode of fixed and rotating crack concepts [2]. The directional stiffness parameter is used to examine singularities of smeared-crack formulations, which signal continuous failure due to material branching. In contrast, the localization indicator is used to examine the formation of weak discontinuities, which precede incipient failure process when the continuum degrades into a discontinuum. In this work, both failure indicators are used to assess the failure modes of the fixed and rotating crack formulations for two loading scenarios — involving pure shear and combined tension with shear.

Layout of combined tension-shear load scenario [2]

Click here to add your own text