Numerical simulation of impulse-induced surface acoustic waves for elastography purposes using k-Wave simulation toolbox Open Access

Elastography is an emerging medical imaging technique that quantitatively assesses and visualizes biological tissues' elasticity and stiffness.^1–8 Because pathological alterations in tissues often lead to changes in tissue mechanical properties, including elasticity or stiffness, elastography can provide more diagnostic information in the form of tissue mechanical properties than traditional medical imaging methods such as ultrasound B-mode imaging, x-ray imaging, and MRI. In addition to disease diagnosis, elastography can also be applied to treatment monitoring^9–11 and image-guided therapy.^12,13

According to the tissue mechanical information provided, elastography techniques can be categorized into two major types: qualitative and quantitative.^14–17 Qualitative elastography measures indirect elastic parameters, such as strain, for diagnosis purposes, while quantitative elastography measures actual elasticity, such as shear modulus and Young's modulus, to interrogate a tissue's physiological and pathological conditions. Although both qualitative and quantitative elastography methods have been introduced to the medical workflow, quantitative elastography has more practical applications in the long run because it provides actual properties of the tissue under test. In most technical development efforts, quantitative elastography relies on generating and tracking shear waves, such as in shear wave elastography^2,18–20 and magnetic resonance elastography,^4,6,21,22 to quantify tissue elasticity for medical diagnosis. Because shear wave is a type of bulk wave propagating inside a tissue at some distance away from the tissue interface, shear wave based elastography can only image deep tissue and cannot image surficial tissue. To interrogate the elastic properties of the superficial tissue, surface acoustic waves have been utilized. The initial technical development is mainly in the form of optical coherent elastography^23–26 and recently expanded to ultrasound surface acoustic wave elastography.^27–30 

Typically, numerical simulation plays an important role in the development of medical imaging techniques because it can provide valuable information for understanding the biophysical mechanisms of imaging processes and optimizing the technical parameters for specific applications. For elastography method development, numerical simulation can help avoid potential complications due to the hardware and software variations of different platforms or manufacturers, which potentially hinder its broader clinical adoption.^31–33 So far, experimental studies are the primary approach in developing surface acoustic wave based elastography techniques, and no numerical simulation method has been developed. To fill this gap, this study aims to provide a numerical simulation platform using k-wave, a MATLAB-based open-source acoustic simulation toolbox, for simulating surface acoustic waves for elastography purposes.

k-wave was initially developed as a numerical tool to perform a more realistic wave simulation and image reconstruction in biomedical photoacoustic tomography.^34,35 Later, its simulation capacity was enhanced by adding models for simulating nonlinear ultrasound propagation³⁶ and elastic wave propagation³⁷ in absorbing media. So far, k-wave can simulate wave propagation in both homogeneous and heterogeneous³⁸ media, in one, two, and three dimensions, and in time forward and reverse manners. Due to such versatile simulation capabilities, k-wave has been selected as a computational investigation tool for various applications, including elastography. To date, k-wave has been used to simulate shear wave generation and propagation in shear wave elastography^39–41 and to simulate soft tissue in strain imaging,⁴² but has not been reported to be applied to simulate surface acoustic waves. This work explores the possibility of using k-wave to simulate surface acoustic waves.

In the rest of this paper, we will first introduce the simulation procedures for generating acoustic radiation force impulse and the corresponding elastic waves, including surface acoustic waves, in a two-medium model. We will then briefly describe the methods used to theoretically and experimentally verify the results. Then, we will provide a detailed description of the results with a brief discussion of the major aspects of this study. We will conclude the paper with a brief summary of the major findings.

This section introduces the methodology of simulating acoustic radiation force (ARF) impulse and its resultant elastic waves, including surface acoustic waves, using the k-wave acoustic wave simulation toolbox in the MATLAB (R2023a, MathWorks, Massachusetts, USA) environment.⁴³ The introduced methodology is applicable for simulating any ARF-induced surface acoustic waves. Out of our research interest, this study will use the Scholte wave as an example to introduce the method procedures. In the meantime, a medical application of surface acoustic waves, i.e., Scholte-wave-based elastography,²⁸ is considered in this study when selecting specific simulation details such as medium materials and transducer parameters. The general procedure of the simulation process is first introduced, and then, the critical steps are explained in detail.

The physical process of acoustic radiation force impulse-induced surface acoustic wave generation can be treated as a two-step process: ARF generation and surface acoustic wave generation. In the first step of the process, an ultrasound transducer excited by an electric pulse generates an ARF impulse in the region of a fluid-tissue interface. In the second step, the medium displacement caused by the ARF impulse induces surface acoustic wave motions propagating along the fluid–tissue interface. According to this physical process, the simulation process is divided into two stages, compressional wave simulation and elastic wave simulation, as shown in Fig. 1. In this simulation, there are three key steps: (i) defining the two-phase media with a flat interface where the surface wave will be generated, (ii) simulating the compressional wave field from which the ARF impulse can be derived, and (iii) simulating the elastic wave field at the fluid–tissue interface region. The second step can be achieved using the k-wave function kspaceFirstOrder3D, and the third step using the k-wave function pstdElastic2D or pstdElastic3D. The details of these three key steps and other associated steps are described in the following subsections.

Unlike the conventional shear wave simulation, which is performed in a single medium model, the surface acoustic wave simulation must be performed in a two-phase medium model. This medium model consists of liquid (water was selected for this study) and tissue for Scholte-wave-based surface acoustic wave elastography, as demonstrated in Fig. 2.

In the k-wave simulation, the following material properties are needed for each medium: mass density, compressional wave speed, shear wave speed (only in tissue), the power law absorption coefficient and absorption exponent. In this study, the material properties for Scholte wave simulation are listed in Table I. Please note that the artificially chosen shear wave speed at a tissue does not accurately represent the actual tissue. This intentional choice facilitated reasonable computing time and memory management throughout the simulation. Please also note that the power law absorption exponent for water is not 1.5, but 1.5 was taken in the simulation because the current version of the k-wave function can only take a scalar number for the whole medium and does not allow for different values for different media.

Compressional wave simulation is the first stage of surface acoustic wave simulation, and its purpose is to simulate the ARF impulse. The predefined ultrasound transducer driven by an electric impulse can generate compressional wave propagation in the two-phase media. The resultant acoustic field of this propagating compressional wave is the source of the ARF impulse, which will further excite the tissue, generating surface acoustic waves at the liquid–tissue interface. The detailed simulation process is introduced below.

The initial step in a k-wave simulation is specifying the computational grid properties required to define other simulation details. In k-wave, a function makeGrid is predefined for this purpose. This function takes two input parameters: the number of grid points $( N x, N y, N z)$ and grid point spacing $( d x , d y , d z)$ in each Cartesian direction. Please note that throughout the simulation process, the $x$-, $y$-, and $z$-direction refers to the lateral, elevation, and axial directions with respect to the ultrasound transducer (Fig. 2), respectively. The grid points denote specific spatial locations within the computational space where the governing equations are solved. In the makeGrid function, the grid points and spacings are used to generate matrices containing wavenumbers and Cartesian grid coordinates. Subsequently, an object of the kWaveGrid class is generated, which contains parameters such as directional wavenumbers, total number of grid points, grid spacings, grid coordinates, number of spatial dimensions, and spatial frequency utilized by the utility functions within k-waves. k-wave also utilizes a specific type of anisotropic absorbing boundary layer known as the perfectly matched layer (PML). PML is important for ensuring sufficient absorption to significantly attenuate outgoing waves and prevent wave reflection from the boundary to the medium.

For our compressional wave simulation, we used a three-dimensional grid of 492 points in the axial/depth direction, 176 points in the lateral direction, and 176 points in the elevation direction without PML. The spatial step size in all three directions was selected as 0.03 mm (approximately 10 grid points per wavelength) to ensure a high spatial resolution for demonstrating the generated waves. Hence, the overall grid dimensions were 14.76 mm in depth, 5.28 mm in lateral, and 5.28 mm in elevation. To model the two-phase medium comprising water and tissue, as demonstrated in Fig. 2, we designated the upper 200 points (6 mm) in the z-direction with water properties, as described in Table I, while the remaining bottom 492 points (8.76 mm) were assigned with tissue properties. A 10-point PML in all directions was used around the computational domain.

Once the kWaveGrid object is created, the time array for the simulation duration can be defined using the built-in function makeTime. This function considers the total simulation time, longitudinal wave speed in the medium, and a non-dimensional number named the Courant–Friedrichs–Lewy (CFL) number. A CFL number is a non-dimensional ratio quantifying the distance a wave propagates within a single time step relative to the spatial resolution (grid spacing). The time step for the simulation is automatically computed using the input of the makeTime function. In our simulation, a CFL number of 0.3 and a total simulation time of 20 μs were chosen. As a result, the automatically calculated time step was approximately 6 ns, and this value was employed to simulate the compressional waves.

Another essential structure for k-wave simulation is the source, which illustrates the properties and location of the acoustic source within the computational space. In k-waves, a source can generally be specified as initial pressure, time-varying pressure, or time-varying particle velocity. Regardless of the source type, two parameters must be defined: a source mask outlining grid points associated with the source and the properties of the excitation signal. An ultrasound transducer can also serve as a source by specifying appropriate physical transducer parameters. However, assigning grid points to individual physical transducer elements and correctly applying delayed input signals to each element's points can lead to complex indexing challenges.⁴³ To address this, k-wave offers a dedicated transducer class that efficiently manages the creation of masks and the assignment of input signals.

A transducer can be modeled using the kWaveTransducer function, which takes inputs of the computational grid properties and the physical transducer properties. The detailed transducer properties used in this simulation are shown in Table II. We used an L11-4v linear phased array transducer (sold by Verasonics) positioned at $z=0 mm$⁠, centered within the x–y plane. The transducer's physical characteristics, including the number and size of its elements, were specified within the simulation. Notably, element sizes were defined in units of grid points, implying a direct relationship between the transducer's physical dimensions and the grid point spacing $( d x , d y , d z)$ in each Cartesian direction.

To define the transducer as an ultrasound source, an input signal must be defined to drive each transducer element. The beamforming delays are automatically computed according to axial focus, elevation focus, and steering angle. In our simulation, the input signal is a tone burst signal with a frequency of 5 MHz and a duration of 1000 cycles. As detailed in the k-wave manual,⁴³ the input signal, represented as a 1D vector, must be a time-varying velocity (or force) source along the axial direction.

As a final structure of the k-wave simulation, sensor points need to be defined to record the acoustic field at each time step in the simulation. The placement of these sensor points within the computational domain is determined by utilizing a sensor mask. Typically, a binary mask with the same size as the computational grid is assigned as a sensor mask. This mask can occupy any location within the computational grid to capture the acoustic field of interest, but in this study, we put the mask in the zero-elevation plane, capturing the dominant sound field generated by the transducer. After defining the sensor mask, k-wave enables recording acoustic parameters at the mask during each time step of the simulation process. Time-varying pressure, particle velocity, and intensity were recorded at the specified mask for our simulation.

Once the four input structures—grid, media, source, and sensor—have been specified, the simulation functions, kspaceFirstOrder3D, can be executed to simulate compressional wave propagation in the defined three-dimensional space. This function is developed based on a first-order k-space model, which allows for the simulation of wave propagation and interactions with structures or boundaries considering those four input structures. Our simulation takes approximately 1 h and 50 min on an Intel Xeon 3 GHz PC with 64 GB RAM (Intel, Santa Clara, CA). After the simulation, the values of the specified acoustic field parameters at the grid points defined by the sensor mask are saved at each time step.

Calculating ARF is the final step and the goal of compressional wave simulation. The calculated ARF will be the input of the elastic wave simulation because it is the source of surface acoustic waves in the physical process. This subsection will illustrate a comprehensive procedure for computing the ARF, utilizing the time-varying pressure and component-wise particle velocity recorded at the sensor mask during simulating compressional wave propagation. As Prieur and Sapozhnikov⁴¹ summarized, three methods exist for modeling ARF in an inviscid fluid: the second-order (SO) approximation, the quasi-plane wave (QPW) approximation, and the attenuated plane-wave approximation. The formulations of these three methods are shown in Eqs. (1)–(3), respectively.

In this study, the ARF components were acquired using all three methods for comparison. The compressional wave simulation can obtain all the required variables, i.e., particle velocity components, pressure, and acoustic intensity, at each grid point (Sec. II C 4). For the elastic wave simulation, we specifically chose the force from the zero-elevation plane as input. This plane contains the peak ARF and is the monitoring plane for mapping the elasticity field in physical experiments.⁴⁷ 

Elastic wave simulation is the second stage of the surface acoustic wave simulation, and its purpose is to simulate the elastic waves, surface acoustic wave included, propagating in the two-phase media. The input of this stage of simulation is an ARF impulse, and the outcomes are propagating elastic waves. However, because the k-wave elastic wave function (pstdElastic2D or pstdElastic3D) does not take force as a valid input, the ARF impulse must first be converted to its equivalent version of particle velocity.

With a proper representation of ARF allowed by k-waves, the next step of the simulation is to redefine the simulation inputs involved in the physical process of elastic wave generation and propagation. Like the compressional wave simulation process, the simulation for the generation and propagation of elastic waves requires four main structures—grid, medium, source, and sensor. Because our focus is on the elastic waves in tissues instead of water, we adjusted the spatial domain of the simulation to cover mostly the tissue and only a small portion of water close to the interface. In simulation, we reduce the axial space $( N z)$ to 236 grid points, in which 20 points (0.6 mm) were selected from the water layer and 216 points (6.48 mm) from the tissue layer. To capture the propagation of the generated elastic waves, we extended the lateral space $( N x)$ to 320 grid points (9.6 mm). The simulation time was extended to 5500 μs, providing enough time for observing surface wave propagation. The elastic wave simulation utilized the same media properties outlined in Table I. Instead of using a transducer, the source structure for elastic simulation was formed by the time-varying velocity components derived from the ARF components at each grid point on the zero-elevation plane, as shown in Eq. (4). Since the space of interest in elastic simulation is different from that in compressional wave simulation, the source data (the time-varying velocity components) are only populated in the overlapped space of these two simulations, and the source values at the rest of the space were set to zero. To mimic a force impulse, the source data were defined to last for 20 μs.

As the last step, the simulation process can be completed by the k-wave built-in functions, pstdElastic2D and pstdElastic3D, which conduct the time-domain simulation of elastic wave propagation within a two-dimensional and three-dimensional medium, respectively. In such a simulation, the medium can be either homogeneous or heterogeneous. As we are only interested in observing the surface wave propagation within the zero-elevation plane, we used pstdElastic2D instead of pstdElastic3D. With the redefined inputs of computational grid, medium, source, and sensor, this function returns the time-varying pressure, component-wise particle velocity, and intensity at the sensor grid points.

To assess the reliability of the introduced method for simulating elastic waves, the simulation results were verified by comparing both theoretical expectations and experimental results. First, the wave speeds of both Scholte wave and shear wave were evaluated from the generated waves and compared with their theoretical expectations. Then, the generated waves and their propagation were compared with the experimental results acquired under conditions similar to those in the simulation.

Once the elastic simulation was completed, we evaluated the most important properties of the generated waves—the speeds of both Scholte wave and shear wave. The evaluated shear wave speeds will be compared with their theoretical expectation, i.e., the shear wave speed adopted in the simuation, as defined in Table I. The evaluated Schotle-shear wave speed ratio and the penetration depths will be compared with their theoretical expectation.²⁷ 

Another method to verify the feasibility of k-Wwves for simulating surface acoustic waves is to compare the simulated waves with the experimentally acquired waves under similar conditions. In this study, the simulated waves are compared with the experimental observations of Scholte wave generated at the water-gelatin interface reported in the literature.²⁸ As reported in this literature, a tissue-mimicking phantom is made of 5 w/v % gelatin from porcine skin (Sigma-Aldrich, Co., St. Louis, MO, USA) and 0.2 w/v % silica gel (Product No. 288500, Sigma-Aldrich, Co., St. Louis, MO, USA). The experiments were made on a programmable ultrasound research system (Vantage 128, Verasonics Inc., Redmond, WA, USA) using a Verasonics linear array transducer L11-4v. The phantom was excited by an acoustic radiation force impulse (ARFI) at 6 MHz, and the generated waves were recorded at 9 MHz. The water layer between the transducer and the phantom is 5 mm, and the focal length of the transducer is set at 7.5 mm away from the transducer surface.

Please note that not all the experimental parameters are the same as those in the simulation. The major differences are the material properties of the tissue-mimicking phantom used in the experiment and the tissue model used in the simulation. Because this comparison is not meant to be an evaluation of the accuracy of the numerically simulated waves but the feasibility of the simulation method (k-waves) in simulating surface acoustic waves, the comparison will be focused on whether surface acoustic waves can be generated using the proposed numerical method or not.

In this section, we will first report the acoustic field generated in the compressional wave simulation, and then compare the ARF fields calculated using three different ARF formulations. Followed is the report of the Scholte wave generated at a water–tissue interface as well as the shear wave generated inside the tissue. By analyzing the wave propagation process, two critical properties of the generated Scholte wave, i.e., wave speed and penetration depth, are reported demonstrating the feasibility of the Scholte wave in imaging superficial tissue elasticity. Finally, the experimentally measured Scholte waves and shear waves are reported and compared qualitatively with the simulated waves.

The acoustic field of the transducer defined in Sec. II C 2 in the water–tissue medium is crucial for the surface acoustic wave generation, and it was first calculated in the compressional wave simulation. Theoretically, the acoustic field can be represented by various physical properties of the medium, such as pressure, particle velocities, and acoustic intensity. In this study, the acoustic field was expressed by the normalized root-mean-square acoustic pressure and its spatial distribution in the zero-elevation plane $(y=0)$ is shown in Fig. 4(a). The axial (z axis) pressure profile located at $x=0$ in Fig. 4(a) is shown in Fig. 4(b).

According to Table I, the acoustic impedances (the product of sound speed and mass density) of water and tissue are 1.480 $× 10 6rayl$ and $1.617× 10 6rayl$⁠, respectively. Because these two acoustic impedances are very close (about 8% difference), the strength of the acoustic field in water and tissue does not show high contrast. However, this acoustic impedance difference leads to noticeable reflections from the interface (z = 0 mm) existing in water (z < 0 mm), and these reflections are visible as lines in the pressure map [Fig. 4(a)] and amplitude fluctuations in the pressure profile [Fig. 4(b)].

The focal region of the acoustic field shown in Fig. 4 is located in the region from −0.81 to 2.91 mm in the z axis, with the peak at z = 0.75 mm. This generally matches the designed focal region, which is 1.5 mm in the z axis (7.5 mm away from the transducer surface). The main reason that the peak amplitude is not exactly at z = 1.5 mm is that k-wave uses the compressional wave speed of water to calculate the phase shifts for the transducer elements involved, neglecting the compressional wave speed difference between water and soft tissue.

The spatial distributions of ARF calculated using the second-order approximation, the quasi-plane wave approximation, and the attenuated plane wave approximation are demonstrated in Figs. 5 and 6. Figure 5 shows the distribution of component-wise forces $f z$⁠, $f x$⁠, and $f y$ at a zero-elevation $(y=0)$ plane in a water–tissue medium. For each method, the three component forces were normalized by the maximum magnitude of the three forces. It is evident in Fig. 5 that axial force $f z$ contains most of the wave energy and, therefore, contributes much more to the resultant force than the other two components. Although small, $f x$ has more energy than $f y$⁠. Like pressure distribution, the ARF field attains its highest amplitude at the focal point.

Among the three methods, the force fields obtained from the attenuated plane wave approximation contrast with those obtained using the other two methods due to the extremely low amplitude of force in water. Disregarding this difference, the axial force $f z$ obtained from all three methods is nearly identical, and the lateral force $f x$ and the elevational force $f y$ obtained with the quasi-plane wave approximation and the attenuated plane wave approximation share more similarity than those obtained with the second-order approximation.

To provide a quantitative comparison of the ARF generated using the three methods, the actual ARF fields (the combination of $f z$⁠, $f x$⁠, and $f y$⁠) with a unit excitation (at each transducer element) at the zero-elevation plane are shown in Fig. 6. In Fig. 6, a vector plot was overlapped with each force field demonstrating the force direction (arrow direction) and amplitude (arrow length) in the axial–lateral plane.

By examining the vector plots, we can see that $f z$ is the main force for all three methods because the major directions of all arrows are downward (along the z axis), especially in the region around $x=0$⁠. A closer examination of the arrow directions can also show that the quasi-plane wave approximation and the attenuated plane wave approximation produce similar results, which are less divergent than the force direction distribution of the second-order approximation. The force direction in the x–z plane, as in Fig. 6, is determined by the relative size of the force components $f z$ and $f x$⁠; these observations are consistent with the results shown in the first two columns of Fig. 5, which clearly show (i) $f z$ is almost identical among the three methods, and (ii) $f x$ provided by the quasi-plane wave approximation and the attenuated plane wave approximation are very similar and smaller than that provided by the second-order approximation.

For the force amplitude, the color intensity of force distribution in Fig. 6 shows that the second-order approximation and the quasi-plane wave approximation methods produce forces with the maximum amplitude $(∼4× 10 4)$⁠: two orders of magnitude higher than the attenuated plane wave approximation (the maximum amplitude being $∼4× 10 2$⁠). Although the force amplitudes of the second-order approximation and the quasi-plane wave approximation methods are in the same order of magnitude, the maximum force amplitude of the quasi-plane wave approximation $(∼5× 10 4)$ is nearly 70% higher than that $(∼3× 10 4)$ of the second-order approximation.

To further demonstrate the difference in the amplitude of the ARF calculated using the three simulation methods, their lateral amplitude profiles at z = 1 mm (going through the focal region) are plotted in Fig. 7. It is clear in Fig. 7 that all three methods provide similar lateral force distribution profiles but dramatically different force amplitudes: the second-order approximation and the quasi-plane wave approximation methods produce forces with similar amplitudes, which are much higher than those produced by the attenuated plane wave approximation.

So far, Figs. 5 and 6 have shown the spatial distribution of ARF accumulated over time, but the real-time ARF field is not like those in these figures. To demonstrate the evolution of ARF in time, the instantaneous ARF distributions in space at three different time points were obtained using the second-order approximation method, as shown in Fig. 8. It is evident in Fig. 8 that as time goes on, the ARF expands in space and grows in magnitude, accumulating energy in the focal region. It should be noted that the instantaneous ARF spatial distribution is strongly related to the length of the source signal exciting the transducer elements. In this study, the source signal is a 5 MHz tone burst signal of 1000 cycles, which is a $200μ s$ long signal and can occupy 300 mm in space (much larger than the z axis dimension in Fig. 8) with a sound speed of 1500 m/s. This is why, in Fig. 8, the time evolution of ARF behaves like a continuous wave instead of a pulsed wave.

In response to impulsive ARF excitation, elastic waves can be generated in the tissue and propagate away from the excitation site over time. Given a water–tissue interface in this study, both Scholte wave (surface acoustic wave at a liquid–solid interface) and shear wave can be generated, as shown in Fig. 8. In Fig. 8, the ARF impulse was applied vertically (along the z axis) at $x=0 mm$⁠, and the generated waves propagate laterally along the x axis. In this study, the elastic waves were presented by $v z$⁠, the local particle velocity in the axial direction. The wave amplitudes were normalized by the maximum wave amplitude generated in their respective simulation methods for comparison purposes.

Figure 9 shows that, for all three methods, elastic waves can be successfully generated. Early after the acoustic radiation force impulse (the left column of Fig. 9), the generated Scholte wave and shear wave propagate together without distinct separation. During propagation, the two elastic waves start to separate due to their different speeds, as in the middle column of Fig. 9. Eventually, as in the right column of Fig. 9, these two waves completely separate from each other, each propagating at different depths of the tissue: Schotle waves propagating in the near-surface region and shear wave propagating in the deeper region of the tissue. To compare the three methods in elastic wave generation and propagation, the elastic waves in Fig. 9 were normalized by the maximum amplitude of the generated waves of their respective methods. Thus, the amplitude difference of the elastic waves generated using different methods can be neglected, allowing us to focus on comparing the spatial distribution of the waves generated using these methods. By comparing the generated waves side by side, we can see that these three methods have no difference in simulating elastic wave generation and propagation. This observation is the same as the shear wave simulation alone.⁴¹ The generation and separation of surface and shear waves after ARF excitation were also observed in a phantom study.²⁸ 

This study used the time-of-flight method to determine the Scholte and the shear wave speeds. Considering the Scholte wave's dominance in the close vicinity of the water–tissue surface, a horizontal line in Fig. 3(a) was established at 0.1 mm below the water–tissue interface to sample Scholte waves for speed evaluation. Similarly, a horizontal line 1.5 mm below the interface was chosen to evaluate the shear wave speed. The Scholte wave speeds and shear wave speeds acquired from the simulations using three ARF modeling methods, i.e., the second-order approximation, the quasi-plane wave approximation, and the attenuated plane wave approximation, were listed in Table III. In the meantime, the average and standard deviation and their ratio were also evaluated. Moreover, the speed ratio between the Scholte and shear waves was evaluated for verification purposes.

Table III shows that all data, including speed measurements and statistical parameters, are consistent among the three methods. The low ratio between the standard derivation and the average speed (STD/AVE), which is around 0.03% for all cases, further suggests the high consistency of the speed measurements in each case. For all the methods, the Scholte and shear wave speeds were estimated to be approximately 8.5 and 10.1 m/s, respectively. Compared with the shear wave speed adopted in the simulation, 10 m/s, as shown in Table I, the measured shear wave speed is slightly higher (1%) than the expectation but in an acceptable range. Furthermore, the Scholte-to-shear wave speed ratio for all three methods is approximately 0.84, closely matching the theoretical value of 0.839.^48,49

To assess the feasibility of k-waves in simulating surface acoustic waves, the simulated waves were compared with the experimentally acquired waves under similar conditions. Because surface acoustic waves at a water–tissue interface have been experimentally observed in our pioneering work of establishing Scholte waves as a surface acoustic wave elastography approach,²⁸ we will compare the simulated waves with the experimental observations. Figure 10 shows the experimentally generated elastic waves under similar simulation conditions. Again, for feasibility assessment purposes, instead of accuracy evaluation of the simulation against experiments, we only selected three typical snapshots of the elastic waves observed as in Fig. 9: (i) Scholte waves not separated from the shear wave, (ii) Scholte waves nearly separated from the shear wave, and (iii) Scholte waves fully separated from the shear wave.

This study demonstrated three different approaches to estimating acoustic radiation force (ARF). The ARF fields generated by the three methods are presented in Figs. 5–7. A close comparison of these three figures shows that, in the tissue (not water), these ARF fields share similar energy distribution (or field shape) but distinctly different amplitudes. It also shows that the ARF field amplitudes calculated by the second-order approximation and the quasi-plane wave approximation are at the same order of magnitude, which is two orders of magnitude than that calculated by the attenuated wave approximation.

Although the purpose of this study is not to determine the accuracy of the ARF estimations using different methods, it is still meaningful to ask the question: Why are the ARF amplitudes so much different among these three methods? As discussed by Prieur and Sapozhnikov,⁴¹ the major difference among these three methods is the applicable conditions of approximation. The second-order approximation has no assumption on the wave geometry for simplifying the mathematical expressions of the ARF source wave, but the other methods have explicit assumptions for the wave geometry. The quasi-plane wave approximation assumes that the wave is a bounded but wide quasi-plane wave, and the attenuated plane wave approximation assumes that the wave is a purely plane wave. Considering that the compressional wave field for ARF generation in this study is highly focused, it is reasonable to conclude that the second-order approximation is a more suitable ARF model, and therefore, the ARF field calculated by this method is more reliable.

For acoustic radiation force (ARF) simulation, the primary interests are the spatial distribution and amplitude of the generated forces. As discussed in the section above, the three methods produce similar spatial force distributions but highly different force amplitudes. Because elastic waves are induced by ARF, the ARF properties will affect the elastic wave properties. First, the wave geometries generated by different methods are nearly identical, and this can be clearly observed in Fig. 9. Second, the elastic wave amplitudes are different among the three methods, although this information is eliminated from Fig. 9 by normalizing the wave amplitudes of each method with their respective maximum amplitude values. For elastography purposes, the simulated wave amplitude is not of interest because its accuracy has no effect on the actual wave amplitude in experiments and so will not affect the reliability of numerical simulation as a tool for the effective study of surface acoustic waves in elasticity imaging applications. As an important wave property for elastography, wave speeds provided by the three different methods are nearly identical. In summary, the three methods are equally eligible for simulating elastic waves for elastography purposes.

The critical question to answer in this study is the feasibility of k-waves as a numerical tool for simulating surface acoustic waves. The results shown in Fig. 9 first confirm that surface acoustic waves, Scholte waves in this study, can be generated in the simulation using all three methods. The comparison between the wave speeds evaluated from the simulated waves and their theoretical expectations (Sec. IV D) suggests that numerical simulations are reliable in modeling the wave speed, the main feature of wave propagation. A qualitative comparison of Fig. 9 (numerical results) and Fig. 10 (experimental results acquired under similar conditions) demonstrates the similarity in the process of elastic wave generation and propagation, further confirming the capability of k-waves in simulating surface acoustic waves. In summary, the successful simulation of surface wave generation and propagation validates the k-wave simulation tool as a viable approach for modeling surface acoustic waves and surface acoustic wave elastography.

Although the validity and effectiveness of the proposed surface acoustic wave simulation method are demonstrated in a water–tissue model (Fig. 2) with a flat interface and homogeneous tissue in this study, the simulation methodology can effectively be extended to other and more complicated scenarios because k-wave allows for an arbitrary distribution of the heterogeneous material properties parameters defined by a user.³⁴ Therefore, the proposed method is not limited to a liquid–solid interface (Scholte wave generation, as in this study) but can be applied to both air–solid (Rayleigh wave) and solid–solid (Stoneley wave) cases by choosing appropriate material properties for the corresponding medium. It can also simulate surface acoustic wave propagations in heterogeneous media and at interfaces of any geometry other than a flat interface.

As mentioned in Sec. II C 4, our compressional wave simulation took about 110 min on an Intel Xeon 3 GHz PC with 64 GB RAM (Intel, Santa Clara, CA). However, the same computational task can take different times to complete under different computational resources. The computational system and platform and the type of codes used determine the variations in computational efficiency. First, the k-wave MATLAB toolbox can be implemented on all platforms where MATLAB is available, and the hardware performance of the platform primarily determines its efficiency. To improve the computational efficiency for large-scale (large grid sizes) simulation in 2D and 3D, optimized versions of the functions kspaceFirstOrder2D and kspaceFirstOrder3D written in C++ are also available for 64-bit Linux and 64-bit Windows systems. Depending on the exact properties of the system and the simulation domain sizes, these optimized codes typically outperform the MATLAB code several times.⁴³ 

The method proposed in this study simulates the physical process of surface acoustic wave generation and propagation in two steps: first, the acoustic radiation force is simulated, and then the resulting elastic waves are. This two-step approach simplified the numerical simulation process but inaccurately represented the physical process. In reality, elastic waves are generated immediately when the acoustic radiation force impulse is applied to the tissue. As the acoustic radiation force impulse interacts with the tissue from the interface to deep inside, the resulting elastic waves, especially the generated shear waves, will have a tilted wavefront propagating obliquely in the x–z plane because these waves are generated at different depths (the z axis) of the tissue at different times. In practice, considering the two orders of magnitude difference between the speeds of the excitation wave (longitudinal wave with a speed of ∼1500 m/s) and the generated wave (shear wave with a speed of ∼5 m/s), the oblique wave can be approximately treated as a plane wave horizontally propagating in the imaging plane (x–z plane in this study). However, in the application where the speeds of the excitation wave and the generated elastic wave are close, or the actual geometry of the wavefront is critical, this two-step method should be modified.

Meanwhile, the proposed method also inherits the inherent limitations of k-waves. The core of k-waves is the pseudo-spectral method for acoustic and elastic wave simulations. While this method provides high accuracy in capturing wave dynamics, it presents technical challenges. The major challenge of k-wave is its significant computational cost, particularly for large-scale simulations, due to the inherent global operations, such as FFTs, in the method.^50,51 Additionally, the global nature of the basis functions results in considerable memory requirements, especially in three-dimensional simulations or when achieving very high resolution, is necessary.^50,51

Although the similarity between Fig. 9 (numerical results) and Fig. 10 (experimental results acquired under similar conditions) confirmed the feasibility of the proposed method, their difference still leaves the question about the accuracy of the proposed method in simulating surface acoustic waves unanswered. For some research projects where the accuracy of the surface acoustic wave simulation is vital, a comparison study should be conducted with well-controlled experimental and simulation conditions.

Another effort to assess the proposed method's capability in simulating surface acoustic waves is to conduct a sensitivity analysis. It will be important to identify the key factors that affect the properties of the generated surface acoustic waves and quantify their effects. Such analysis will provide practical guidelines for future numerical investigation of surface acoustic wave applications, including surface acoustic wave elastography as a new medical imaging strategy.

This study introduces the k-wave acoustic simulation toolbox as a numerical simulation platform for simulating surface acoustic waves induced by acoustic radiation force (ARF) impulse. Based on the physical process of the ARF-induced surface acoustic wave generation, we divided the simulation process into two stages: compressional wave simulation and elastic wave simulation. The goal of the first simulation stage is to find the spatial distribution of the ARF impulse, and the second stage of simulation is to generate elastic waves and record their propagation in time and space. Through the simulation on a water–tissue medium model, we first demonstrate the feasibility of k-waves as a numerical tool for the surface acoustic wave simulation. In the meantime, we compared the performance of three ARF modeling methods, i.e., the second-order approximation, the quasi-plane wave approximation, and the attenuated plane wave approximation, in ARF generation and elastic wave generation. Our results show that these three methods provide similar spatial distribution of ARF fields but different ARF amplitudes. Despite the difference in simulating AFR fields, these three methods provide nearly identical results regarding elastic wave generation and propagation in space and time domains. Considering elastography's requirement for the reliability of wave propagation in time and space, the methodology developed in this work can be an effective and low-cost platform for the further development of surface acoustic wave-based elastography methods.

The authors are grateful to Texas Tech University Graduate School for providing the Doctoral Dissertation Completion Fellowship to partially support this work.

The authors have no conflicts to disclose.

Abdullah A. Masud: Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal). Jingfei Liu: Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Methodology (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.