Superposition method for modelling boundaries between media in viscoelastic finite difference time domain simulations

Finite-difference time domain (FDTD) methods have been used extensively to model the propagation of ultrasound in heterogeneous media.¹ The use of FDTD simulations has been validated for planning transcranial ultrasound treatments,² which is an emerging modality for treatment of brain disorders including thalamotomy for treatment of essential tremor³ and chronic neuropathic pain,⁴ and thermal ablation of brain tumours.⁵ Applying transcranial focused ultrasound presents several challenges due to the presence of the skull, which is a problematic biological tissue to focus through due to its strong acoustic mismatch with water, conversion between compressive and transverse waves, scattering, and attenuative properties, which gives rise to phase aberrations, standing waves, and undesired tissue heating.^6–9 

The Virieux-staggered grid FDTD scheme¹⁰ is widely used to model the viscoelastic wave propagation in heterogeneous media^11,12 due to its computational efficiency, smaller memory requirement, and easy implementation,¹³ and has been validated for the case of shear mode conversion.¹⁴ Staggered grid methods provide greater accuracy and stability than traditional FDTD methods, and form the basis of reverse-time migration and full-waveform inversion.^12,15

A source of error in FDTD solutions of the viscoelastic wave equation arises when the liquid-solid interface does not align with the Cartesian grid.¹⁶ This staircase effect was first described by Fornberg.^11,17 Staircasing refers to the spatial approximation required to define continuous geometries onto a discrete Cartesian grid in two and three dimensions (2D and 3D), resulting in many vertices or edge positions that do not correspond to points on the grid.¹ This error is exacerbated by end-point misregistration, which occurs when the angle of the interface does not form a Pythagorean triangle on the grid, leading to line-source endpoints that are not coincident with specific grid point positions. Both the staircase representation of a continuous line and the endpoint misregistration are considered aspects of the staircase effect.¹ The discrete discontinuities in the representation of the interface produces circular wavefronts emanating from points where the discrete staircase approximation causes the interface to appear “jagged” or discontinuous on the rectangular grid, which act as a source of unphysical noise and artifactual waves.¹⁷ 

Previous approaches to satisfy the boundary conditions at a material interface include the use of a hybrid grid-/finite-difference- (FD-) method, where the grid method¹⁸ is applied in the fluid for its flexibility in handling complex fluid-solid interfaces, while the FD method is applied in the solid due to its lower numerical dispersion.⁹ Other approaches may account for surface topography, such as the technique developed by Tessmer et al.,¹⁹ which is based on a pseudo-spectral method that maps a rectangular grid onto a curved surface, but requires significantly more time and memory and is unstable for large surface curvature. The FD method on a rotated staggered grid²⁰ uses rotated finite-difference operators to align the grid to the material interface. Finite differences on deformed grids^21–24 are used for their relative simplicity and preservation of stencil computations,²⁵ but are incompatible with the use of the fully staggered Virieux scheme, and so require more computationally expensive approaches. Lisitsa et al.²⁵ developed a combined Galerkin/FD method that exploits the polyhedral mesh of the discontinuous Galerkin method in the regions around the interface in order to more accurately handle complex topography, while employing the more efficient FD method in the rest of the computational domain. The use of finite element methods provides flexibility for their ability to accurately represent complex geometries, but construction of irregular polyhedral meshes is a problem of high complexity, and the use of global finite element methods is more computationally expensive than FDTD methods on deformed grids.²⁵ 

Previous examinations have demonstrated that this staircasing effect decreases with finer spatial sampling, and that reasonably accurate results may be obtained from the standard FDTD method by using grid point spacing of 30 points per minimum wavelength (PPW).²⁶ For 2D simulations, the memory requirements increase proportional to Δx^–2, where Δx is the spatial step in the FDTD implementation, while computational time requirements increase proportional to Δx^–3, which imposes limits on the use of high-resolution FDTD algorithms in practical applications.

Due to these limitations of these numerical schemes, there is a demand for solutions that can reduce the staircasing effect with manageable increases to the memory and time requirements, while maintaining the simplicity, accuracy, and efficiency of FDTD methods. To address this problem, we present a modification to the standard Virieux staggered-grid FDTD scheme that increases the accuracy of simulations without the constraining resources required for increased spatial resolution. Conceptually, we applied the principle of superposition for linear systems²⁷ to represent the liquid-solid interface as a superposition of different configurations. Our results demonstrate that this approach is especially advantageous to address the problem of staircase artifacts in FDTD modelling.

The proposed method is inspired by the principle of superposition, which states that for linear systems the overall state may be broken into a sum of separate pieces, where the equation can be solved in each separate piece, the pieces may be added back together to obtain the solution for the whole.²⁷ The goal of this method is to address the ambiguity of the liquid-solid interface by having the locations of the step discontinuities distributed at different locations along the interface in different simulations by using different approximations of the interface based on the location of the interface relative to the grid positions, and then combine the results of those different simulations to obtain a solution that contains contributions each point along the interface.

The superposition method aims to reduce the artifactual wavefronts produced by the discrete staircase approximation of angled or curved interfaces by recombining the waves from multiple simulations, where each simulation exhibits differing approximations of the interface based on the location of the interface relative to the grid positions. The application of the superposition method will be considered analytically in 1D, and numerically in 2D.

The staircasing effect in FDTD simulations arises when the location of the material interface does not align with the discrete grid used in numerical simulations. In order to examine the application of the superposition method analytically in 1D, the transmission and reflection of waves is considered from an interface that can only be defined at discrete points that are offset from the location of the physical interface. To begin, supposed that a boundary between two media exists at x = 0. Incident wave ψ_i at frequency f₀ is propagating toward increasing x, giving rise to the reflected wave ψ_r and transmitted wave ψ_t. This scenario is shown in Fig. 1. The wave speed c is given by

The incident, reflected, and transmitted waves will then be given, respectively, by

where the time-dependence has been suppressed. The wave numbers are given by

and the reflection and transmission coefficients satisfy

The staircasing effect in FDTD simulations arises when the location of the material interface does not align with the discrete grid used in simulations. To approximate the staircase problem analytically in 1D, the location of the interface is restricted to discrete points that are offset from x = 0. When these points do not align with the location of the interface, the interface can either be located at the point that immediately precedes the location of the interface, or the point immediately after. For the first case, the location of the material interface is considered to be at position $x=−(1−β)Δx$⁠, while in the second case the location of the material interface is located at $x=βΔx$⁠. The quantity β can range from 0 to 1, and represents the fraction of the distance between the two grid points that is occupied by the solid. The length $Δx$ is analogous to the spatial resolution in numerical simulations. The transmitted wave in scenarios 1 and 2 are denoted $ψt1$ and $ψt2$⁠, respectively, and are shown in Fig. 2.

The goal of the superposition method is to use a linear combination of the basis functions $ψt1$ and $ψt2$ to create a transmitted wave $ψt′$ that is mathematically more similar to the exact transmitted wave ψ_t than either of the basis functions. For 1D, the resultant transmitted wave will be given by

Compared to waves that are transmitted through the interface at the actual location ψ_t, the transmitted waves in the situations shown above experience a phase shift, such that

and

where

To quantify the correlation between the reconstructed transmitted wave $ψt′$ and the analytical transmitted wave ψ_t, consider an error function e_t that is given by the difference between the reconstructed wave and the analytical transmitted wave, such that

The relative error for the transmitted wave, ε_t, can be given as the magnitude of the error function, relative to the magnitude of the transmitted wave, which can be found using

which for the transmitted wave reconstructed using the superposition method yields

This can be compared with the error $ϵtst$ calculated from the stepwise approximation to the transmitted wave, where the location of the interface is assigned to whichever grid point is closer to the location of the actual interface, and the transmitted wave is given by

Solving for the magnitude $ϵtst$ yields

Because the stepwise approximation preserves the amplitude of the reconstructed waves and differs from the exact solution only by a phase shift, all error for the waves in the stepwise approximation can be attributed purely to phase error and not amplitude error. The same analysis was applied to quantify the reconstructed reflected wave error $ϵrSPP$ which yields

This can be compared against the error for the reconstructed reflected wave in the stepwise approximation $ϵrst$⁠, which yields

While the error for the transmitted wave reconstructed using the superposition method does not depend on the direction of propagation, the reflected wave error depends on the wave number in the incident/reflected medium, which for a given interface will depend on the direction in which the wave is incident on the boundary. This same analysis may be applied for a wave that is reflected from the opposite side of the boundary due to an incident wave propagating towards –x to obtain the error, $ϵ−r$⁠. For the reflected wave reconstructed using the superposition method, this yields

while for the wave reflected from the opposite side of the interface calculated using the stepwise approximation, this yields

To evaluate the magnitude of the error, wave speeds were selected based on the properties of water and acrylic, such that $c1=1482 m/s$ and $c2=2848 m/s$⁠. To be consistent with numerical comparisons, the magnitude of $Δx$ is given by 1/6th of the smallest wavelength, and a frequency $f0=1 MHz$ is selected, giving a $Δx$ value of 0.247 mm. The magnitude of the transmitted and reflected wave error for the superposition method and the stepwise approximation are shown in Fig. 3. The error is lower for the superposition method than for the stepwise approximation for all wave components and all values of β, and converge with the exact solutions as the value of β approaches 0 or 1.

To evaluate the effects of amplitude and phase error separately, the waves reconstructed using the superposition method given in Eq. (7) allows for the reflected and transmitted waves to be written in terms of the corresponding analytical wave as

where $Aij$ is a complex amplitude (with $i=t,r,−r,$ and $j=SPP,st$⁠). The amplitude and phase of A for the superposition and piecewise approximations are shown in Fig. 4. The piecewise approximations to the various wave components only differ from the analytical solutions by a phase difference, so that the absolute value of the amplitude coefficient is always equal to 1, and so are not included in the figures. While the superposition approximation shows reduction in the overall amplitude of the various waves, the overall reduction in total error by the superposition method compared to the stepwise approximation shown in Fig. 3 can be attributed to the significant reduction in phase error achieved by the superposition method.

The Virieux scheme is a finite difference method of modelling viscoelastic wave propagation on a fully staggered grid with the P-SV formulation that defines the viscoelastic equations given as a first-order hyperbolic system according to^10,14,28

and

where ρ is the density, and λ and μ are the Lamé coefficients. The r_ij values are memory variables used to include the standard linear solid model of viscoelasticity, which are given by

and

where τ_l and τ_s are the stress relaxation times for the longitudinal and shear stress, respectively, and $τσ=1/2πf0$ with f₀ given by the central frequency. The full procedure for optimizing material parameters based on speed of sound and attenuation coefficient is given in Ref. 28.

The viscoelastic wave equation was solved using fourth order in space and second order in time^29–31 finite-difference operators. Harmonic averaging of the bulk and shear moduli and arithmetic averaging of the density was applied in all simulations,³² which provides greater accuracy across fluid-elastic boundaries even when the interface is aligned with the grid.³³ The full procedure for the implementation of the FDTD algorithm is given in Ref. 31.

The superposition method is a modification to the Virieux FDTD method. It consists of running multiple simulations of the physical space, which will be referred to as “configurations,” where each configuration contains a representation of the physical space constructed with different staircase approximations of the material interfaces. Measurements of the stress/velocity are based on the average values across all configurations.

For the 1D solution, the error in the location of the interface is related to the quantity β, and there are only two possible waves that can be recombined, allowing for the precise definition of the superposition weighting values used in Eq. (7). In contrast, the error in the location of the numerical interface in 2D numerical simulations can vary continuously for each cell that is intersected by the interface. The superposition weighting is approximated by averaging multiple simulations of the same domain with a different staircase approximation for each configuration. For a cell that is intersected by the interface, the fraction of configurations in which that cell is designated a solid is proportional to the amount of space in the cell occupied by the solid.

For FDTD simulations, each grid point must be defined as either a solid or a fluid, with the corresponding material properties. In 2D this is done by creating a map of material properties, designated $χ(x,z)$⁠, which is assigned a value of 0 to denote the fluid or 1 to denote the solid. The physical properties at each point in the grid are determined by the value of the material map at that point, such that

with a similar relation for the Lamé parameters and relaxation coefficients, which determines the physical quantities (ρ,λ,μ,τ) at each point in space.

For a given cell position, the fraction of configurations in which a cell is designated as the material that lies beyond the interface is proportional to how much of the space within that cell is beyond the line defining the physical interface. Implementation is done using a weighting map $β(x,z)$ defined at each point in space, where the value of β is given by the fraction of space within each cell that is beyond the line of the interface, and building material maps with a configuration-specific threshold that must be met for the cell to be designated as the material that lies beyond the interface. The grid coordinates are considered to be the location of the longitudinal stresses in the cell, such that the cell at coordinates (x₀, z₀) is considered to range from $x0−Δx/2≤x<x0+Δx/2$ and $z0≤z<z0+Δx$⁠. For regions where only a single material is present, the values are 0 or 1. At a material boundary that does not align with the Cartesian grid used for the simulation, the weighting values fall between 0 and 1. This is shown in Fig. 5(a).

For FDTD simulations, only a single χ map is required. For the superposition method, a unique map is created for each separate configuration. For a solid/fluid interface where the solid lies beyond the line of the interface, each element in the $nth (n=1,2,…,N)$ material map χ_n will be designated as the solid if the weighting value (β) satisfies

where N is the total number of configurations. This threshold ensures that the number of configurations where a given cell is designated as a specific material is approximately equal to the relative area that that material occupies in that cell. For the example shown in Fig. 5(a) with N = 4 configurations, the corresponding χ_n maps are shown in Fig. 5(b). For N = 1, the superposition method reduces to the unaltered Virieux FDTD scheme, with a material map constructed using Eq. (30) with a threshold of 0.5.

The FDTD method is applied to each configuration independently. The resulting values of pressure/stress are given by the average value across all configurations, for example the stress components σ_ij it yields

and root-mean-square (RMS) stress amplitude maps used for error quantification are calculated according to according to

To determine the accuracy of the superposition method, the model was tested under three different conditions, as shown in Fig. 6. The first scenario uses a highly simplified arrangement with a half space occupied by water and a lossless solid. The second scenario uses the same fluid/solid half space arrangement but introduces attenuation in the solid half-space. For the third and final test, the scenario consists of transmission through layered attenuating media, which provides a more practical application of the method.

The simulations for various values of incident angles (θ), spatial resolutions (⁠$Δx$⁠), and number of configurations (N) were compared to high spatial-resolution simulations. Low spatial resolutions (6, 8, 10, 12, and 14 PPW) were compared against high spatial resolution simulations (50 PPW), which provide sufficient spatial resolution to produce reasonably accurate results for all incident angles.²⁶ 

The physical scenario is shown in Fig. 6. A solid/fluid interface is located 5 cm from the source plane. The solid has the longitudinal/shear speed of sound, longitudinal/shear attenuation (⁠$αl/αs$⁠) and density of acrylic as determined experimentally in Ref. 14. For the lossless half space simulation, the wave speeds and density of the solid are unchanged, and attenuation is set to zero. The surrounding medium is a lossless liquid with the density and longitudinal speed of sound of water. These quantities, as well as the corresponding relaxation values, are given in Table I. For trials with the layered attenuating media, the thickness of the solid layer was selected to be 8 mm, which is within the measured range of thickness of human skulls.^34,35

The size of the domain was 10 cm along x and 11 cm along z. Spatial step $Δx$ was set to $minλs/PPW$ based on a frequency of 1 MHz and shear speed of sound given in Table I, while the time step $Δt$ was calculated to incorporate the effects of attenuation on the stability criteria³⁶ and provided a Courant-Friedrichs-Lewy number of 0.7. The domain is surrounded by a perfectly matched layer³⁷ (PML) to absorb outgoing waves. Detectors were placed at the edges of the domain (just before the start of the PML) in order to detect the outgoing wave. To vary the incident angle θ, the material interface was tilted at various angle values with respect to the z-axis.

The source was located at z = 0 m and generated waves propagating towards the solid material. The source plane was aligned with the Cartesian grid to ensure that no artifacts arise from staircase approximation of the source plane, which is not addressed by this method. The impact of staircasing on acoustic line-sources has previously been examined in Refs. 43 and 44. An unfocused source with a Gaussian amplitude distribution along x was used. To minimize the spectral width, the source signal S(x, t) in time was a Gaussian wave packet modulated at a central frequency f₀ of 1 MHz, given by

where $τ=4/f0$ and $a0=1 cm.$

To quantify the staircasing effect, the 2D RMS stress amplitude maps were compared against the RMS stress amplitude maps from high-resolution simulations. High resolution data were downsampled using a nearest-neighbour method to match the same spatial points as those from low-resolution simulations. Each RMS amplitude map $σRMS(x,z)$ was compared with the downsampled high-resolution simulations $σ̂RMS(x,z)$ using the structural similarity index (SSIM), which has been validated for spatial model assessment due to its inclusion of mean, variance, and correlation between two continuously valued maps.³⁸ The structural similarity is calculated over the simulation domain using³⁹ 

where μ_x and μ_z are the average pixel values over x and z, η_x and η_z are pixel value standard deviations along x and z, and η_xz is the pixel value covariance along x and z. C₁, C₂, and C₃ are constants introduced to stabilize instabilities when any denominator is close to zero. The mean SSIM index is calculated over all points in space to evaluate the similarity over the entire map. The error was then quantified using the structural dissimilarity index (DSSIM), which is related to the SSIM as⁴⁰ 

which has a range in the interval of $[0,2]$⁠, where two identical maps will have a DSSIM value of zero, while maps that are negatives of each other will have a DSSIM value of 2. The mean square error (MSE) was also calculated between the RMS amplitude maps and the downsampled high resolution maps using the formula

where K is the total number of spatial points. The MSE error was calculated for different values of incident angle, spatial resolution, and number of configurations. While the DSSIM provides a relative metric for correlation between the two maps, the RMS MSE provides a direct quantitative metric of the absolute difference between the two RMS amplitude distributions.

The mean square error Ξ_S was also calculated for the detectors along the boundaries of the domain using

where M is the total number of samples and x is the spatial coordinates. The high-resolution detector signals were downsampled to match the space and time coordinates of the low-resolution detector values.

Simulation code was implemented using the cuda language (V9.0.176, Nvidia, Santa Clara, CA, USA) and executed in a high-performance computing cluster counting with graphic processor units (GPUs) (GP100, Nvidia, Santa Clara, CA, USA) available at Compute Canada facilities. matlab software package (R2019a, MathWorks, Natick, Massachusetts, USA) was used to prepare and analyze data for simulations.

Examples of the reduction in scattering artifacts in the lossless solid/fluid half-space can clearly be seen in the RMS stress amplitude map as shown in Fig. 7. The Virieux method introduces multiple scattered wave artifacts that are significantly reduced when the superposition method is applied.

The reduction of scattering artifacts is reflected in the reduced MSE error values as a function of N shown in Fig. 8, and the reduced DSSIM as a function of N shown in Fig. 9. These results show that the spatial distribution of the reflected and transmitted waves is improved with the application of the superposition method for all incident angles and spatial resolutions. The increase in error with increasing value of N shown for some incident angles may be attributed in part to the reduction of amplitude in the reflected wave with increasing value of N that is caused by destructive interference between waves in different configurations. As the value of N increases, more unique solutions are present which gives rise to more destructive interference between configurations that reduces the overall amplitude, which manifests as a magnitude error in the reflection and transmission coefficients. While this interference introduces error in the magnitude of the reflection and transmission coefficients, the 1D superposition method Sec. II A indicates it is also the source of the significant reduction in imaginary error that is the source of scattering artifacts. This is consistent with the analytical consideration of the error arising from the staircase approximation and effects the superposition are presented in  Appendixes A and  B, respectively.

The MSE calculated for the detectors around the boundaries of the domain show more limited reduction of error with increasing value of N, as shown in Fig. 10. This can be explained by the time delay that is introduced in low-resolution simulations due positional error in the detector positions and to the effects of numerical dispersion, to which this error metric is particularly sensitive,²⁶ and which is not addressed by the superposition method.

The RMS stress amplitude maps are shown in Fig. 11, where the effects of attenuation in the solid can be seen to cause rapid attenuation of the transmitted wave. The MSE as a function of N is shown in Fig. 12, which illustrates that the superposition method is less effective in reducing the error when attenuation is present in the solid. This can be explained with reference to the attenuative effects in the solid, which reduces the overall contribution of transmitted waves to the MSE and causes the MSE to be more heavily weighted towards the error in the reflected waves. When compared with the results of the lossless half space test, this lends further support to the conclusion drawn from the 1D analytical case, where both amplitude and phase errors from the interface are most significant in the waves that are reflected from the interface. While the RMS MSE shows some inconsistency in the ability of the superposition method to reduce scattering artifacts, the DSSIM metric shown in Fig. 13 for the attenuating half space illustrate how the application of the superposition method improves structural characteristics in the spatial distribution of the reflected and transmitted waves.

As was the case with the lossless half space, the MSE for the detectors at the domain boundaries show limited reduction in error with the application of the superposition method, which may be attributed to the same explanation as the detector MSE in the lossless half-space.

Examples of the reduction in scattering artifacts in the lossless solid/fluid half-space can clearly be seen in the RMS of the stress amplitude map as shown in Figs. 14–16. The Virieux method introduces multiple scattered wave artifacts that are significantly reduced when the superposition method is applied, even in cases where multiple interfaces are present.

The MSE was also calculated for the RMS of stress amplitude map exclusively in the regions beyond the solid plate, where only the transmitted wave is present. These results are shown in Fig. 17, which shows a reduction in the MSE with increasing value of N due to the reduction of artifactual transmitted waves.

The application of the superposition method significantly reduces the effect of staircasing in viscoelastic FDTD simulations. The case of the fluid/solid half-space provides a simplified model for verification of the ability of the superposition method to reduce scattering artifacts in 2D that arise from the staircase approximation. The attenuating half space demonstrates that the superposition method is compatible with the inclusion of viscoelastic relaxation for modelling attenuative effects. The third case of transmission through layered media in 2D demonstrates the applicability of more practical and complex scenarios, like transcranial focusing, which may include multiple interfaces.

The RMS of the stress amplitude maps show that the profile of the reflected and transmitted waves are less affected by the diffractive effects from the staircase artifacts. The reduction in artifactual scattered waves is reflected in the reduced DSSIM for the RMS stress amplitude maps, which show a reduction in error with increasing number of configurations N. While the RMS MSE metric also shows a reduction in error with the application of the superposition method, the improvements are generally more limited, especially when the MSE is calculated over the whole simulation domain, which may be a product of the ability of SSIM to incorporate structural information from continuous valued maps, where such additional information is helpful in spatial model assessment.³⁸ 

The differences between the SSIM and the MSE reflect the improved distribution of energy from reflected and transmitted waves, at the expense of precision in the energy carried in those waves. The ability of the superposition method to improve the distribution of energy while reducing wave amplitude can be seen in the difference maps showing the absolute error, along with the corresponding high-resolution RMS map, shown in Fig. 18 for 5° and 40° interface angles.

The 1D case presented in Sec. II A illustrates how the superposition method significantly reduces phase error, while introducing an error in the magnitude of the transmission and reflection coefficients, and that the error is highest in the reflected wave due to the magnitude of the wavenumber in the incident medium that produces a larger phase difference between configurations. As shown in Fig. 12, the superposition method produces RMS MSE error greater than uncorrected case for some interface angles when the transmitted wave is rapidly attenuated in the solid. It was explained in Sec. IV B that the attenuative effects in the solid reduces the overall contribution of transmitted waves to the MSE, causing the error to be more heavily weighted towards the error in reflected waves. The MSE metric is considered to be sensitive to changes in amplitude, compared to structural changes (such as the elimination of the diffraction artifacts), which was among the reasons why the DSSIM was included. The increase in error with increasing N shown in Fig. 12 can be attributed to the reduction of the reflection coefficient caused by the interference between different configurations, the attenuation of the transmitted wave that minimizes its contribution to the overall error, and the sensitivity to changes in amplitude that the MSE metric exhibits.

 Appendix A provides an analytical examination of how phase errors arise from staircase approximations.  Appendix B uses the same analytical framework to examine how the superposition method reduces the imaginary component of the error, while introducing an error in magnitude. This is consistent with the results shown in Fig. 18 where scattering artifacts caused by phase error are significantly reduced, but the amplitude of the reflected and transmitted waves are also reduced. The results of the analytical approach in  Appendix B provides a basis to conclude that as the value of N increases, the magnitude error caused by the application of the superposition method will converge to a finite value, which decreases with increased spatial resolution.

As mentioned in Sec. IV A, the detector MSE showed limited improvements with the application of the superposition method. This may be attributed to time delays that can arise due the effects of numerical dispersion during propagation as well as to positional errors in the detector positions, to which this error metric is particularly sensitive,²⁶ and which is not addressed by the superposition method. This delay can be seen in the pressure vs time in Fig. 19, where the transmitted waves in the low-resolution simulations experiencing a slight delay. Additional methods are available for further reduction of numerical dispersion, but aside from the higher order FD approximation for differential operators, they have not been included here.⁴¹ While detector MSE metrics are presented here for completeness, dispersion analysis has previously been excluded from examinations of interface errors.⁴¹ Despite this limited improvement of the detector MSE, the RMS map error metrics are more effective at capturing the erroneous spatial distribution of waves that are produced by the staircase approximation, and which the superposition method is intended to address.

The memory required to implement the superposition method scales proportionately to the number of configurations. Results indicate that improvements in accuracy are minimal beyond $N≥10$⁠. An advantage of the superposition method is that each different simulation configurations can be run to completion, and the final signal can be progressively reconstructed as each configuration is simulated. This is in contrast to high-resolution simulations, where all stress, pressure, and memory variables at all points in space must be stored in memory at once. Additionally, while the time required for a standard FDTD 2D simulation scales with $Δx−3$⁠, the total time required to complete a simulation for the superposition method scales with $NΔx−3$⁠, with a larger value for $Δx$⁠. For the 50 PPW simulations, each stress, velocity, and memory variable component must be calculated at $2.12×1012$ different points in space and time, which requires $2.67×107$ values be stored in memory for each variable array. For the lowest resolution (6 PPW) simulations with N = 8, accounting for all the different configurations, the values of the stress, velocity, and memory variable component must be calculated at $3.17×1010$ different points in space and time, requiring only $4.13×105$ values to be stored in memory for each variable array at a time.

Both the Virieux FDTD and superposition methods show a significant reduction in the transmitted amplitude of shear waves at high incident angles when compared to high-resolution simulations. This is evident by the decreased amplitude of the transmitted waves in 40° RMS figures shown in Fig. 14, as well as an increased amplitude in the waves reflected from the first material interface. This error is especially pronounced at incident angles above the longitudinal critical angle, indicating a reduced transmission coefficient for shear waves. Plane wave analysis has previously indicated that the numerical reflection and transmission coefficient error is purely imaginary, implying that the amplitude of the coefficients remain the same but introduce a phase shift.⁴² These results suggest that this analysis may break down for shear wave conversion and that an error in the magnitude of the reflection/transmission coefficients is introduced due to the staircase effect.

Some inconsistencies in the relationship between the error and the spatial resolutions, and specifically the increased error at higher spatial resolutions, as can be seen in Fig. 17, may be explained by diffractive effects that arise from the periodicity of the staircase approximation, the misalignment of the staircase representation of the liquid-solid interface¹ and the compounded effect of modelling multiple interfaces simultaneously, which has been shown to produce higher errors than simulations which include only a single interface.¹ 

Lack of error reduction for angles approaching 45° is an expected characteristic of the superposition method. The number of unique material map representations is limited by the number of possible positions available for the step discontinuity at a given value of x, which is proportional to $1/ tan θ$ for angles between 0° and 45°. As N increases, many of the new material maps are exact duplicates of already existing map in other configurations. This is apparent in the 6 PPW N = 20 simulation, where there are only three unique material maps configurations for the top interface between all 20 configurations, while for simulations employing low θ values, all 20 material map configurations are unique. Cases of identical configurations are only expected to occur when the interface angle is constant and close to 45°. In these cases, the limited ability of the superposition method to improve results may be mitigated by the natural error minima that occurs at 45° due to the interface aligning with grid points.²⁶ In realistic simulations, it is expected that curved interfaces will result in a range of interface angles, making the duplication of exact configurations unlikely.

For cases in which identical material map configurations are produced, an increased value of N will lead to the proportion of configurations in which a given cells is designated a specific material approached the exact fraction of space within that cell that the material occupies. Increasing the value of N is therefore expected to provide improved accuracy even when there is a limited number of unique material map configurations, though the improvements in accuracy will be limited compared to increases in N that include unique material maps. Redundant calculations may be avoided without altering the basis of the superposition method by checking for configurations that are identical, ensuring that these configurations are only used in a single simulation, and using a weighted averaging when recombining the final field.

The superposition method is effective in reducing the staircasing artifacts, and does not require any modifications to the implementation of the Virieux FDTD algorithm in individual simulations. For this reason, the superposition method is compatible with the application of other other numerical modifications, such as the spatial averaging developed by Moczo et al.³² 

While the Virieux method has been verified to provide reasonably accurate results at all incident angles with 30 PPW spatial resolution,²⁶ some artifactual scattered waves are still present at low incident angles, as shown in the maps in Fig. 20, which illustrate the difference between the high resolution RMS maps and the various low-resolution RMS maps. These scattered waves correspond to the first order diffraction maxima that are produced by the diffraction grating-like structure of the interface arising due to the periodicity of the staircased approximation. While these artifacts remain in high-resolution simulations, they are not present in low-resolution simulations where the superposition method is applied. While other artifacts remain in the low-resolution superposition RMS amplitude map, this result may indicate that at small incident angles the superposition method is capable of reducing the diffractive effects that arise from the staircasing error effect more effectively than increasing the spatial resolution, which is consistent with the larger reduction in error that occurred with smaller tilt angles. For this reason, the spatial resolution for high-resolution simulations used for comparison in this paper was increased from 30 PPW to 50 PPW, which was the maximum resolution available with current resources, and where no artifactual waves were observed.

Significant reduction in the staircasing artifacts in viscoelastic FDTD simulation can be achieved by the application of the superposition method. The examples provided demonstrate a significant reduction in the artifactual scattering of waves at the solid/fluid interface, without any alterations to the FDTD method applied in any of the individual simulations.

This research was enabled in part by support provided by Compute Ontario and Compute Canada. We acknowledge the support of the Discovery Program of the Natural Sciences and Engineering Research Council of Canada (NSERC).

In this appendix, the effect of staircasing on the sampling of a propagating plane wave is considered analytically. To begin, a sensor line shown as the dotted line in scenario 1 in Fig. 21.

An incident wave with wave vector $k→$ is propagation along $−ẑ$ towards the sensor line given by $z=−x tan θ$⁠. The field at that line, $S0(x)$⁠, is given by

where $k=2π/λ$⁠.

In scenario 1 shown in Fig. 1, the tilted sensor line is approximated by a discontinuous staircase approximation with step size 2a along $ẑ$ and step width w along $x̂$⁠, shown as a solid line. The effect of staircasing introduces a position-dependent phase error, due to the spatial misalignment of the actual sensor line and the staircase-approximated sensor. For the staircase approximation to the exact sensor line, this phase delay of the sensor staircasing can be described by

where

Based on the diagram in Fig. 21, the periodicity of the staircasing error creates a periodic phase error that from $−w≤x≤w$ can be given by

which satisfies the requirement based on Fig. 21 that $ϕ(−w)=1$ and $ϕ(w)=exp (−2ika)$⁠.

In the rotated coordinate system $(x̂′,ẑ′)$ shown in Fig. 21 with $x̂′$ parallel with the sensor line, such that $x=x′ cos θ$⁠, Eqs. (A2) and (A4) can be written as

and

where $w′$ is the distance along the sensor line from $x=−w$ to x = w, with

The saw-tooth-like phase error characteristic of $ϕ(x′)$ caused by the periodicity of the staircase error will produces multiple discrete k-vector components along the $x′$ axis, which can be found by writing $ϕ(x′)$ in terms of its Fourier series components, which yields

with

Solving for the wave number along the $x̂′$ axis yields

which is equivalent to the grating equation with a periodicity equal to $2w′$⁠.⁴⁵ 

This demonstrates that the staircased approximation of the sensor line gives rise to artifacts that can be expressed as a discrete set of k-vector components along $x′$⁠. In the case of the staircased approximation to the sensor, the exact solution corresponds to the n = 0 order, and it can be shown that the amplitude of the zeroth order is given by $sinc(ka)$⁠, with all other orders correspond to staircasing artifacts. This same behaviour has been observed in the diffraction-grating-like behaviour of reflected waves when the interface between two media is approximated by a staircase approximation. This is illustrated in Fig. 22, which show the RMS wave amplitude map from a 3° interface with the incident wave removed, and is consistent with the behaviour predicted by the analytical diffraction equation.

Consider now the case where the field is sampled at two different staircased approximations of the sensor plane, separated by a distance 2a, as shown in Fig. 21. In scenario 1, the staircased approximation of the sensor is above the sensor line, while in scenario 2 the staircased approximation is below the sensor line. The field at these planes can be given by

and

The case for scenario 1 has already been examined in  Appendix A, which gives

$ϕ1(x′)$ and $ϕ2(x′)$ differ by a positional difference along z that produces a phase shift, which allows $ϕ2(x′)$ to be written as

To apply the superposition method, first a superposition weighting factor $β(x′)$ is introduced, which gives the relative distance from the first staircased sensor to the actual sensor line

This weighting factor is used as coefficients to combine the two fields $S1(x′)$ and $S2(x′)$ into the superposition-weighted field, $S′(x′)$⁠, given by

or

with

In contrast with the individual fields which satisfied $|ϕ1|=|ϕ2|=1$ for all values of x, in this case

Using the values functions for $ϕ1(x′), ϕ2(x′)$⁠, and $β(x′)$⁠, it can be shown that $ϕ′(x′)$ can be represented in a Fourier series as

For the zeroth order term, we can determine that

The two staircase approximations giving in Fig. 21 represent the two basis functions for the superposition method. However, scenarios 1 consistently overestimate the position of the sensor position along z, and scenario 2 consistently underestimates the position along z. A more accurate staircase approximation would intersect the interface at x = 0, and would provide

with the first term amplitude, a₀, given by

The parameters in these equations can be expressed in terms of relevant FDTD variables by substitution $a=Δx/2, Δx=λ/PPW$⁠, and $k=2π/λ$⁠. For both the staircased approximation and the superposition approximation, the amplitude of the zeroth order is reduced. In the staircase approximation, the reduction of intensity in the central order is due purely to the energy scattered into higher diffraction order modes. With the application of the superposition method, the combination of the different basis functions leads to destructive interference that reduces the amplitude of all modes, but disproportionately affects the n > 0 modes. The absolute, real, and imaginary components of the amplitude for $|n|≤10$ modes are shown in Fig. 23(a), 23(b), and 23(c). The effect on the n = 0 mode amplitude as a function of spatial resolution is shown in Fig. 23(d).