Transcranial ultrasound modeling using the spectral-element method Open Access

Accurately modeling the propagation of ultrasound waves is imperative for both ultrasound computed tomography (USCT) and focused ultrasound (FUS) applications. Solving the wave equation within high-contrast media, such as between soft tissue and bone, introduces several considerable complexities that can otherwise be neglected for systems consisting entirely of soft tissue.

First, the significant impedance contrast at the soft tissue-bone interface acts as a strong reflector for the incident ultrasound wavefield, which limits the amount of energy that can be transmitted across the bone layer (Dines , 1981; Fry and Barger, 1978; Ylitalo , 1990). Second, the conversion between compressional (P-) and shear (S-) waves at the fluid-solid interface increases the complexity of the wavefield considerably, especially for waves with steep angles of incidence (Clement , 2004; Marty , 2021; Pichardo and Hynynen, 2007; White , 2006). In addition to the high impedance contrast and fluid-solid coupling, transcranial ultrasound applications introduce the further challenge of requiring the ultrasound waves to propagate through the highly dissipative trabecular (diploë) layer in the skull (Pinton , 2011; Robertson , 2018). This portion of the skull consists of fluid-filled pores, which act to dampen the incoming ultrasound waves, thus further reducing the amount of coherent ultrasound energy that can propagate into the brain tissue.

The combination of these complexities makes it challenging to accurately simulate the propagation of ultrasound waves across the skull. Despite this, recent developments in the area of transcranial ultrasound imaging have shown promise using waveform-based inversion approaches (Guasch , 2020; Marty , 2022b; Mitcham , 2023; Orozco , 2023; Robins , 2023). In contrast to inversion techniques based on travel time measurements, waveform-based inversions make use of the complete time series recorded at each receiver position to better approximate a set of reasonable model parameters that describe the observed data (Fichtner, 2011; Pratt , 2007; Tarantola, 1984). Full-waveform inversion (FWI) is one such technique, which iteratively updates the acoustic and elastic properties of soft tissues and bone by computing sensitivities with respect to a misfit function; a common choice to obtain these sensitivities is to correlate the forward and adjoint wavefields (Fichtner , 2006; Tromp , 2005; Virieux and Operto, 2009). However, for such inversion strategies to be effective, it is paramount that the simulation framework used for computing synthetic wavefields is numerically accurate.

Incorporating these sources of complexity is also of significant relevance when determining the maximum pressure distribution within FUS to ensure that the position, magnitude, and size of the pressure focus are correct. This has significant implications for the safe and effective use of FUS for applications such as neuromodulation (Legon , 2014; Tufail , 2010) or thermal ablation (Elias , 2016; Hynynen and Jolesz, 1998). Several numerical modeling approaches exist for determining the maximum pressure distribution in FUS, many of which use waveform-based modeling techniques (Aubry , 2022).

In addition to the breadth of medical literature on both FWI and FUS, identical topics are also of great relevance within geophysical applications. The process of performing FWI is used for imaging the subsurface on a wide range of spatial scales from local-scale seismic imaging (Pratt, 1999; Shen , 2018) to regional/global-scale seismology (Tape , 2009; Thrastarson , 2024; Thrastarson , 2022). Similarly, quantifying the peak ground velocity (PGV) or peak ground acceleration (PGA) is of significant interest within seismic risk assessments and can be computed using identical techniques to those employed in FUS within medical ultrasound (Baker , 2021; Rodgers , 2019). This abundance of both medical and geophysical literature highlights the importance of these approaches, in addition to the transferability of these types of modeling methods.

There are a plethora of numerical methods that one can use to solve the wave equation. One method which is particularly widely used in regional- and global-scale seismology is a high-order finite-element method known as the spectral-element method (SEM) (Afanasiev , 2019; Faccioli , 1997; Komatitsch and Tromp, 1999; Komatitsch and Vilotte, 1998). The SEM is very well-suited for dealing with geometrically complex domains, while also possessing a number of attractive computational benefits over conventional finite-element methods (Igel, 2017). The SEM has previously been applied in medical contexts in works such as Martiartu (2019); Bachmann and Tromp (2020); Marty (2021), and Ulrich (2022). A key novelty in this work is the consideration of attenuation in both fluid and solid media and the derivation of a spectral-element formulation for the coupled system, which is suitable for matrix-free time domain solvers and explicit time-stepping schemes. Furthermore, leveraging the advantages of boundary-conforming meshes for transcranial applications is explored.

In this work, we present a strategy for precisely capturing the ultrasound wavefield interactions at the soft tissue-bone interfaces through the use of the SEM. The two primary objectives of this work are (1) to provide a formulation of the coupled viscoacoustic-viscoelastic wave equation for the SEM and (2) to illustrate an end-to-end workflow for performing transcranial ultrasound simulations using the SEM. Several numerical examples are performed to show the flexibility of this modeling framework.

This paper is divided into the following sections: First, a brief history of USCT and FUS for high-contrast media is provided along with an overview of the coupled viscoacoustic and viscoelastic wave equations used in the SEM. Next, a summary of the forward modeling workflow is outlined including the hexahedral mesh generation process used for discretizing the domain. A sensitivity analysis using an adaptation of the Shepp–Logan phantom is then conducted to demonstrate the importance for accurately resolving the skull interfaces within the coupled wave equation. A validation of the modeling techniques within the SEM is then provided using the k-space pseudospectral wave solver k-wave (Treeby and Cox, 2010; Treeby , 2020) as a reference. Finally, an in silico demonstration of these coupled viscoacoustic-viscoelastic modeling procedures is given for a realistic skull phantom.

Modeling the propagation of mechanical waves in high-contrast media has been extensively studied in a variety of fields such as medical ultrasound (Wiskin , 2020), seismology (Shen , 2018), and non-destructive testing (Greenhall , 2021). A strong impedance contrast between adjacent materials, such as between soft tissue and bone, leads to a significant reflection of the incoming ultrasound wavefield. The impedance, which is given as the product between the density and the wave speed, is approximately 3.5 times higher in bone when compared to soft tissue (Duck, 1990; White , 2006). This contrast in material properties leads to normal-incidence reflection and transmission coefficients of approximately 0.55 and 0.45, respectively. Thus, in the case of transcranial ultrasound, this severely limits the amount of energy which is transferred into the interior brain tissue.

To motivate the need for full-waveform modeling tools that are both accurate and performant, it is important to discuss the differences between wave-based modeling strategies and ray theory, which is widely used in ultrasound imaging. Ray-based modeling approaches typically approximate the arrival of the ultrasound waves using the travel times of the primary arrivals (Menke, 2012; Ruiter , 2023). This class of techniques has the significant advantage of being computationally inexpensive when compared to solving the wave equation. Historically, this has led to the widespread adoption of ray-based ultrasound modeling strategies due to their fast time to solution (Birk , 2014; Greenleaf and Bahn, 1981; Li , 2009).

While ray-based modeling methods are powerful for a variety of ultrasound applications, the high degree of complexity within the ultrasound wavefield in the presence of high-contrast materials makes the use of wave-based modeling approaches an attractive alternative. In the context of transcranial ultrasound imaging, the first-arriving waves typically stem from either the reflected wave from the surface of the skull or the guided wave that travels through the skull bone (Marty , 2021); both of these first-arrivals are completely insensitive to the interior brain structure.

In addition to ultrasound imaging applications, full-waveform modeling is a useful tool for accurately determining the maximum pressure distribution within transcranial FUS. Waveform modeling tools are particularly well-suited for these applications, given that they can accurately account for the reflections and refractions across the skull boundaries while also being able to capture phenomena such as Scholte waves along the tissue-skull interface or multiple reflections within the skull.

As is evident within Aubry (2022), many numerical codes within the transcranial ultrasound community utilize discretization schemes that are based on regular grids, given that spatially discretizing the domain becomes trivial. While some tools have been developed in recent years based on conventional finite-element (Huthwaite, 2014) or boundary-element (Haqshenas , 2021) methods, the use of the SEM for transcranial applications has thus far been limited to simulations that were either non-attenuative (Marty , 2022a; Marty , 2022b) or purely viscoacoustic (Aubry , 2022).

The propagation of ultrasonic waves in fluid and solid media are fundamentally different due to their resistance to shearing (Aki and Richards, 2002). Fluid media, which have negligible shear moduli and, by extension, possess low S-wave velocities, can be described as acoustic. On the other hand, solid materials, with non-negligible shear moduli, can be described as elastic.

In the context of ultrasound imaging, soft tissues can typically be classified as acoustic and bones can be described as elastic. While elastic waves certainly still propagate through soft tissues, as is evident within techniques such as elastography, the S-wave velocity within soft tissues is on the order of 10 m/s or less (Taljanovic , 2017), while the S-wave velocity for bone is generally within the range of 1350 m/s to 1650 m/s (Pichardo , 2017; White , 2006). This distinction between acoustic and elastic tissues can be made for this particular application, given the temporal scales considered in USCT. That is, the S-wave velocity within the soft tissue is sufficiently low that these shear waves have very little influence in the resulting wavefield when compared to the compressional waves.

The reasoning behind making this clear distinction between acoustic and elastic media is twofold. First, resolving the propagated waves within areas of both low and high velocity equally well would require the same number of degrees of freedom per wavelength for both regions. Due to the inherently shorter wavelengths present within the acoustically slow portions of the domain, an exceptionally fine discretization within these regions would thus be required. This would lead to an immense increase in the overall computational cost while not providing any considerable increase in the signal's information content on the time scales being considered. Second, computing the acoustic (scalar) wave equation is considerably less expensive when compared to the elastic (vector-valued) wave equation (Fichtner, 2011). Given that it is known a priori that the proportion of bone to soft tissue is low, avoiding these unnecessary computations of the elastic wave equations can lead to a significant reduction in the overall computational cost.

Acoustic-elastic coupling has been studied extensively in the context of transcranial FUS (Clement , 2004; Pichardo , 2017; Pulkkinen , 2014; Treeby and Cox, 2014), transcranial imaging (Marty , 2022a; Mitsuhashi , 2017), and seismology (Afanasiev , 2018; Chaljub and Valette, 2004; Komatitsch , 2000; Nissen-Meyer , 2008). The conversion of P- to S-waves and vice versa becomes more pronounced as the angle of incidence between the wavefield and the interface becomes steeper. This phenomenon can be computed analytically using the Zoeppritz equations, where it can be shown that any inbound waves that are not of normal incidence to the fluid-solid interface will have some conversion between P- and S-waves (Yilmaz, 2001). It has been shown in practice that incidence angles of <20° typically lead to minimal acoustic-elastic conversions within transcranial FUS applications (Clement , 2004). However, there is still no clear consensus within the literature over what constitutes acceptable incidence angles for transcranial USCT to allow for these coupling effects to be ignored.

Consider a domain, $x∈Ω⊂Rn$ (where $n=2,3$⁠) with outer boundary $∂Ω≡Γabs$⁠, over a given time interval, $t∈T⊂R$⁠. This domain consists of disjoint regions of acoustic media $Ωa$ and elastic media $Ωe$⁠, which are separated by interfaces $Γi$⁠, as illustrated in Fig. 1. Note that the precise meaning of the outer boundaries $Γabs$ and $Γsponge$ will be discussed in more detail in this section once the absorbing boundaries are introduced.

Both the acoustic and elastic wave equations are parameterized in terms of the space-dependent material properties density $ρ(x)$ and P-wave velocity $vp(x)$⁠, while the elastic wave equation additionally considers the S-wave velocity $vs(x)$⁠. All parameters can vary in space. An acoustic medium is considered to be a fluid where the shear modulus (and, by extension, the shear wave velocity) is negligible, that is, $vs≪vp$⁠.

Dissipation is commonly parameterized using the quality factor $Q(x,ω)$⁠, where the frequency-dependence of $Q(x,ω)$ can be modeled using a system of standard linear solids (Emmerich and Korn, 1987). The resulting attenuation enters the elastic wave equation by modifying the stress $σ(x,t)$ within Eq. (3). This feature, whereby the attenuation enters the spatial derivative term $∇·σ$ instead of the temporal derivative term $ρ∂t2u$⁠, will be of considerable importance concerning the implementation discussed in Sec. II E. To emphasize this, let $σ˘(x,t)$ denote the stress that includes attenuative effects from the viscoelastic medium.

In general, the use of five linear solids has been shown in the literature to yield accurate results for broadband signals (van Driel and Nissen-Meyer, 2014). However, it is possible to achieve a satisfactory level of accuracy using fewer linear solids, provided that the frequency content of the signal is sufficiently narrow. From a practical standpoint, increasing the number of linear solids tends to lead to a minimal increase in compute cost within the time loop, although there is a considerable increase in memory overhead due to the added dynamic fields, which need to be temporarily stored.

In seismic wave propagation, attenuation in the acoustic domain is typically neglected in the coupled system (Komatitsch , 2000; Nissen-Meyer , 2008). However, including such effects for medical applications is of considerable importance given that the majority of the domain predominantly consists of acoustic regions.

A considerable difference between the viscoelastic and viscoacoustic formulations is that the attenuation acts on either the space or time derivative terms, respectively. This has significant implications within the time-stepping scheme when updating the state variables. Details regarding implementation are discussed further within Sec. II E.

In addition to these interface conditions, absorbing boundary conditions are assigned to the outermost extent of the domain on $Γabs$ and in the region of $Ωsponge⊂Ωa$ (see Fig. 1). The absorbing boundary conditions consist of a combination of the Sommerfeld radiation (Sommerfeld, 1912) condition and sponge layers (Kosloff and Kosloff, 1986).

The SEM is a matrix-free high-order finite-element method that uses the tensorized Gauss–Lobatto–Legendre (GLL) basis (Afanasiev , 2019; Faccioli , 1997; Komatitsch , 2000; Komatitsch and Tromp, 1999; Komatitsch and Vilotte, 1998). The primary numerical benefit which the SEM offers is that the global mass matrix is diagonal by construction. The consequence of this is that inverting the mass matrix becomes trivial, meaning that one can use explicit time-stepping schemes without having to solve a linear system.

While this diagonal structure of the mass matrix makes the SEM attractive from a computational standpoint, being able to exploit this unique property requires one to discretize the spatial domain using hexahedral elements (Komatitsch and Tromp, 1999). This requirement may be problematic in some instances as the process of constructing such conforming hexahedral meshes in a (semi-)automatic way is non-trivial and remains an unsolved problem in computational geometry (Pietroni , 2023). The use of tetrahedra does not preclude the use of the SEM, although employing tetrahedra tends to lead to a loss in the diagonality of the mass matrix (Igel, 2017). While mass lumping techniques exist that allow for this diagonality to be preserved using tetrahedral elements (Liu , 2017), computing spatial gradients require one to perform dense operations on these tetrahedral elements (Ferroni , 2017). This is in contrast to the inherently tensorized structure of hexahedral elements where these same operations become sparse.

The acoustic and elastic wave equations must be cast into their corresponding weak formulation to be solved using the SEM. This weak formulation is achieved by multiplying the strong forms of the wave equation [i.e., Eqs. (3) and (17)] by a time-independent test function and integrating each expression over Ω.

Given the weak formulations of the wave equation described in Eqs. (19)–(21) and Eq. (24), one would cast these equations from their integral form to a sum of elemental integrals over a reference element. This step is omitted in the interest of brevity, although the reader is referred to works such as Komatitsch (2000) or Nissen-Meyer (2008) for further details.

The differences in where the viscoacoustic and viscoelastic attenuation enter into Eqs. (26) and (27) has important implications from an implementation standpoint. The memory variables in the viscoacoustic formulation are incorporated within $∂t2ϕ˘$⁠, given the relationship between the acoustic pressure and the second time derivative of the scalar potential, as outlined in Eq. (11). Thus, the acoustic memory variables can be interpreted as modifying the acoustic pressure directly. This differs from the viscoelastic formulation wherein the memory variables enter within the stresses, which corresponds to the elliptic term [see Eq. (3)]. This distinction makes it such that the memory variables enter the stiffness matrix $K˘e$ for the viscoelastic formulation while the memory variables do not enter the mass matrix $Ma$ for the viscoacoustic formulation.

The memory variables are updated using the approach outlined in Sec. 2.4 of van Driel and Nissen-Meyer (2014). These terms are updated for the acoustic and elastic regions during the same step in which the acceleration terms [i.e., Eqs. (30) and (31)] are computed.

The order in which the state variables $φ(x,t)$ and $u(x,t)$ are updated is important, given that $∂t2φ˘(x,t)$ must be known in order to compute $∂t2u(x,t)$⁠. The sequence in which the various fields are updated within the Newmark time-stepping scheme is outlined in Fig. 2.

In order to link the attenuative formulations introduced in Sec. II C with Eqs. (28)–(33), consider the pseudocode outlined in Algorithm 1. The collocation frequencies ω_j, frequency-weighted coefficients a_j, and unrelaxed moduli $CU(x)$ and $κU(x)$ are determined for the N linear solids prior to commencing with the time loop. While updating the displacement and velocity terms is straightforward, updating the acceleration terms involves several smaller intermediate steps.

Updating the acoustic acceleration potential $∂t2φ˘(x,t)$ is achieved by first computing Eq. (30). This acceleration potential is then modified using Eqs. (10) and (11) in order to introduce the attenuation into the acoustic wavefield. The acoustic memory variables $ξ(x,t)$ are then updated by solving a scalar equivalent of the first-order differential equation detailed in Eq. (8).

Computing the elastic acceleration $∂t2u(x,t)$ is quite different from the acoustic case, given that the additional derived quantities of strain $ε(x,t)$ and stress $σ˘(x,t)$ need to first be determined from the displacement $u(x,t)$⁠. The strain $ε(x,t)$ is computed using Eq. (4), from which the attenuated stress $σ˘(x,t)$ can be determined given Eq. (7). After updating the elastic memory variables by solving the first-order differential equation in Eq. (8), the attenuated stress is used to obtain the stiffness term $K˘eu$⁠. The attenuation has already been incorporated within the stress-strain relation, which was used to compute $σ˘(x,t)$⁠, meaning that the stiffness matrix at this stage is already encoding the attenuation (hence the accent over the stiffness matrix). The updated acceleration term $∂t2u(x,t)$ can then be computed without further modification using Eq. (31).

An important consideration within grid-based discretizations is to construct the domain such that so-called staircasing artifacts are mitigated (Afanasiev , 2018; Marty , 2022b; Pelties , 2010). Staircasing artifacts occur when a curved interface is represented on a rectilinear grid, as seen in Fig. 3(a), thus resulting in the appearance of non-physical reflectors. These staircasing artifacts become more pronounced the greater the material contrast is across the boundary. While these staircasing artifacts can be mitigated by reducing the grid size such that the interface appears smooth relative to the wavelengths considered in the subsequent wave simulations, this typically requires the grid spacing to be heavily over-resolved, which can lead to a potentially significant increase in the overall computational cost.

It is important to consider the differences between conforming and non-conforming meshes. Contrary to a conforming mesh, a non-conforming mesh is a discretization within which hanging nodes are permitted. A hanging node is defined as a node which falls on the edge of an element, while not belonging explicitly to the connectivity of said element. This meshing approach allows for a high degree of geometric accuracy to be achieved with minimal additional manual intervention during the mesh generation phase. However, using such an approach leads to the loss of many of the computational benefits that make the SEM such an attractive computational approach. Thus, all meshes considered in this study are exclusively conforming in nature. An example of a non-conforming hexahedral mesh is shown in Fig. 3(b).

A way to achieve a similar degree of geometric accuracy as within a non-conforming hexahedral mesh is to refine the material interfaces such that the conforming nature of the mesh is preserved. Employing such local refinements is generally more complex than applying refinements using non-conforming meshes, although this can still be achieved in a somewhat systematic fashion. An example of a locally refined conforming hexahedral mesh is shown in Fig. 3(c).

Using conforming locally refined hexahedral meshes allows for the mass matrix to remain diagonal, thus preserving one of the most computationally favorable aspects of the SEM. However, applying local refinements along specific material interfaces leads to locally constricted elements, which significantly reduce the global time step due to the Courant–Friedrichs–Lewy (CFL) criterion. The maximum global time step $Δt$ and the minimum element size $Δx$ are related by $Δt∝Δx$ within the CFL criterion, meaning that these refined elements can lead to a significant penalty in the global time step, thus making it prohibitively expensive to calculate a solution to the wave equation.

An alternative to using such a refined hexahedral mesh is to instead construct a boundary-conforming hexahedral mesh. This type of mesh is constructed such that the element boundaries align with the interface which is to be preserved. This allows for elements of relatively uniform size to be maintained throughout the mesh, which means that a penalization in the global time step is avoided. While using such boundary-conforming meshes is ideal from a computational standpoint, the complexity associated with the mesh generation process increases substantially compared to the aforementioned approaches. An example of a boundary-conforming mesh is shown in Fig. 3(d).

It is important to note that there are certainly other discretization strategies in the literature that have also been shown to be effective for dealing with these staircasing artifacts at strong material boundaries, particularly on regular rectilinear grids (Moczo , 2002). However, many of these strategies increase the implementation complexity of the numerical scheme considerably. Thus, including the distinct interfaces within the meshes explicitly serves as a convenient way of depicting the wavefield at these interfaces.

The software tools blender (Blender Foundation, 2023) and coreform cubit (Coreform LLC, 2023) are utilized throughout this study to generate the conforming hexahedral meshes required for use within the spectral-element solver. A brief overview of this meshing workflow is provided here, although the reader is referred to Marty (2022b) for further details.

This discretization strategy involves first constructing a high-quality surface representation of the soft tissue-bone interfaces; these are subsequently used to generate the hexahedral mesh. An implicit assumption for this type of meshing approach is that an approximate surface representation of the bones is available a priori. Such a preliminary surface mesh could be acquired by using a contouring algorithm such as marching cubes (Lorensen and Cline, 1987), which had been applied to a segmented voxel dataset such as a computed tomography (CT), magnetic resonance imaging (MRI), or reflection-based ultrasound reconstruction.

blender is an open-source computer aided design software that is commonly used in animation, motion graphics, and visual effects. blender is used as a preprocessing tool for (1) cleaning up the surface geometries and (2) smoothing the surfaces using its broad range of solid modeling and sculpting tools. The resulting geometry must be manifold (i.e., a closed surface with no self-intersecting faces) to allow for the subsequent meshing process to be applied. The primary use for applying smoothing is to mitigate remnant staircasing artifacts, which may be apparent within the surface mesh as a result of the marching cubes' extraction.

Once a high-quality surface representation has been created, the surface mesh can be used to construct a conforming hexahedral mesh within coreform cubit. The meshing from volume fractions algorithm (Owen , 2014) is used to construct a boundary-conforming hexahedral mesh such that the interfaces between elements precisely follow the surfaces constructed in blender. Meshing from volume fractions generally results in a mesh that consists predominantly of cube-like hexahedral elements with a relatively uniform size throughout the mesh. However, the element boundaries at the material interfaces conform with the location of the initial surface meshes.

This meshing approach has the distinct advantage of using a rectilinear-like mesh for most of the domain, while then locally modifying the connectivity at material interfaces to ensure that the mesh boundaries conform with the boundaries in question. This property is well-suited for domains that contain several regions of similar material parameters, such as is the case with soft tissue and bone. Ideally, the elements in the skull would be larger than the elements in the soft tissue in order to secure a favorable time step given the CFL criterion. While it is, in principle, possible to construct such a mesh, the element sizes in the skull would be of a similar size to the thickness of the skull in some parts of the bone, given the ultrasound source frequencies considered. Thus, this study generally uses a uniform element size for the entire domain to ensure that the geometric complexity of the domain is sufficiently captured.

It is important to mention the difference between this approach of meshing from volume fractions versus conventional hexahedral mesh generation, which is typically performed using domain decomposition. While techniques based on domain decomposition typically split the medium into subdomains such that each subdomain can be meshed independently, meshing from volume fractions acts as a global operation to generate a hexahedral mesh in a semi-automatic fashion. This approach is often less well-suited for meshing entities with sharp corners when compared to decomposition-based techniques, but it serves as a very effective meshing framework for smooth irregular geometries.

The meshes generated in coreform cubit are constructed with first-order accuracy in space; this corresponds to eight nodes per element. As an intermediate processing step, additional nodes that are consistent with the placement of the GLL points within the spectral-element basis are added to each element. All meshes considered for this study are converted to be fourth-order accurate in space such that each element within the computational mesh possesses 125 nodes.

The SEM solver salvus (Afanasiev , 2019) is used throughout this study for solving the wave equation. Simulations are performed for bowl-type transducers using source characteristics that are identical to those utilized in Aubry (2022). These transducers emit a continuous sine source time function with a reference frequency of 500 kHz. The transducers are modeled as a distribution of monopole point sources with a source pressure of 60 kPa.

Unless otherwise stated, discretizations that result in a mesh resolution of approximately three elements per wavelength (EPW) at a reference frequency of $ωref=500 kHz$ are utilized. All simulations are performed with fourth-order accuracy in space and second-order accuracy in time.

The source time function is emitted evenly over the surface of the bowl transducers. These transducers are modeled using a collection of point sources arranged in a phyllotaxis pattern to ensure that the distribution of points over the surface of the transducer is uniform (Vogel, 1979). Note that using the SEM does not require these individual point sources to coincide with the nodes in the computational mesh, meaning that the shape of the transducer does not need to be explicitly meshed.

The forward modeling principles necessary for waveform-based USCT and FUS are identical, as both techniques depend on the realistic modeling of ultrasound waves through heterogeneous and geometrically complex media. However, given that the pressure distributions that are used in FUS consist of a single time-invariant pressure field, comparing these FUS pressure distributions is slightly more convenient and intuitive than interpreting the time-dependent waveforms which are of importance in USCT. As such, the forward modeling examples in this study primarily consist of demonstrations in the context of FUS.

The focal volume is computed for the pressure distribution $pmax(x,t)$ by thresholding the pressures above a certain decibel level. This thresholding tends to isolate several discontinuous regions in the model, typically in close proximity to the focus. The focal volume is computed for the single largest continuous region delineating during this thresholding step.

The full width at half maximum is computed along each of the principal axes (x-, y-, and z-directions) inline with the ultrasound transducer. For domains where the ultrasound transducer is not aligned with the principal axes, a coordinate transform is first applied to the mesh such that the center of curvature of the bowl transducer is located along the z-axis.

In an ideal case, the reference pressure distribution would be computed using an analytical solution, similar to those performed using the focus code (Kelly and McGough, 2009; McGough , 2004) for the homogeneous case within Aubry (2022). However, since all of the cases considered in this work can only be solved numerically, the reference solution is considered to be the numerical solution with the finest level of discretization.

First, a sensitivity analysis was performed using an adaptation of the Shepp–Logan brain phantom (Shepp and Logan, 1974) to demonstrate the importance of resolving the soft-tissue interfaces within the mesh correctly. Two meshing strategies are considered here, namely, rectilinear meshes and cubed sphere meshes (Ronchi , 1996). The interface between the soft tissue and the skull is precisely resolved within the cubed sphere mesh, whereas this interface is superimposed onto the rectilinear mesh. In an effort to emphasize the differences between these various choices of spatial discretizations, the non-attenuative coupled wave equation is considered in this section; the attenuative cases are subsequently explored in more detail within Secs. IV B and IV C. Representative material properties from the literature, as outlined in Table I, were used in order to simulate both soft tissues and bone (Duck, 1990; Pichardo , 2017).

All computations were performed on the Alps supercomputing system from the Swiss National Supercomputing Centre (CSCS). This system utilizes nodes that house four Nvidia Grace-Hopper modules (Nvidia, Santa Clara, CA), where each module contains a 72-core Arm Neoverse V2 central processing unit (CPU) (Arm, Cambridge, UK) with 128 GB of LPDDR5X random access memory (RAM), and an H100 graphics processing unit (GPU).

In order to provide an impression for the compute time required to run each of the simulations in this section, the absolute run times are provided. The number of GPUs used for each simulation varies with the number of elements in the mesh, with approximately 1 GPU being used per 0.5 × 10⁶ elements. Using strong scaling tests, it was found that this ratio yields an excellent balance between time-to-solution and inter-GPU communication. In order to more easily compare the compute costs between simulations, the computational costs in terms of node hours are also provided; this effectively serves as a normalized measure of the compute cost. Note that these computation times are for the time loop of the simulation only and exclude any additional pre- or post-processing of the simulations.

The rectilinear mesh constructed at 2.0 EPW with local refinements has been generated using the unrefined rectilinear mesh at 2.0 EPW as a basis before applying tripling refinements along the expected skull boundaries. This refinement scheme [similar to Fig. 3(c)] involves splitting the sides of each element into three, thus resulting in each element being split into nine in two-dimensional (2D) format or 27 in three-dimensional (3D) format. These tripling refinements preserve the conforming nature of the finite-element mesh while improving the resolvability of the interface.

A bowl transducer that emits a continuous 500 kHz sine wave is placed at a radial distance of 12.5 cm from the origin of the phantom with a declination of 45° as seen in Fig. 4.

The pressure distribution is measured along a plane with normal $n̂=[1.0,−1.0,0.0]$ and which passes through the origin. This corresponds to a sectioning plane, which bisects the center of the bowl transducer and allows for the effective visualization of the maximum pressure distribution emitted by the transducer. All relevant statistics are computed on this slice, which bisects the focus location. A summary of the forward modeling results for the various meshes can be seen in Fig. 5 and Table II.

As can be seen in Fig. 5, the pressure distribution extends behind the transducers. This occurs due to the omnidirectional nature of the individual point sources, which are distributed across the surface of the transducer. A non-physical radiation pattern appears behind the transducer since these point sources radiate energy equally in all directions. However, given that the outer interface immediately behind the transducer corresponds to an absorbing boundary, these non-physical effects do not influence the resulting maximum pressure distribution. These transducers are also modeled as non-rigid sources in the sense that they do not act as scatterers. The consequence of this is that the reflected arrivals from the soft tissues and the skull do not interact with the sources.

It is clear upon inspection of Fig. 5 that the maximum focal pressure in the interior brain tissue is heavily impacted for the cases where unrefined rectilinear meshes are used [Figs. 5(a) and 5(b)]. This is also reflected in the $ℓ2$ and $ℓ∞$ differences illustrated in Table II relative to the cubed sphere reference solution (4.0 EPW), wherein the respective difference metrics are considerably higher than those obtained for their cubed sphere counterparts at the same EPW.

The $ℓ∞$ differences between the rectilinear meshes and the reference cubed sphere solution highlights the challenge with modeling high-contrast curved boundaries on a rectilinear mesh. The point-wise difference is dominated by the discrepancies at the skull boundaries, as is illustrated at the boundary near coordinates $x≈(−0.05,0.05)$ in Figs. 5(a) and 5(b). However, the $ℓ2$ differences across the entire domain do, indeed, decrease with increasing EPW for the rectilinear meshes, thus implying that the agreement across the entirety of the medium improves as the resolution of the mesh increases.

The cubed sphere meshes exhibit a considerably greater degree of convergence with increasing EPW when compared to the rectilinear meshes, both in terms of $ℓ2$ and $ℓ∞$ differences. The point-wise error encoded in the $ℓ∞$-norm shows excellent agreement between the various cubed sphere solutions for $EPW≥2.0$⁠. Similarly, the $ℓ2$ difference for the cubed sphere meshes with $EPW≥2.0$ exhibits very close agreement with the reference solution. This close agreement can be attributed to the strong imprint of the skull being captured with similar levels of accuracy between the various cubed sphere solutions. The deviation for the cubed sphere mesh at 1.0 EPW primarily stems from the accumulation of dispersion errors within the wavefield. These dispersion errors result from the mesh being too coarse relative to the source frequencies which are to be modeled.

Finally, the rectilinear mesh at 2.0 EPW with local refinements along the skull boundaries yields considerably better results than any of the other rectilinear meshes tested, albeit at a considerably higher computational cost. This increase in computational cost is primarily due to a significant penalty in the global time step incurred from the locally restricted elements within the refined mesh. As is evident in Fig. 3(c), the reduction in element size due to the use of these local refinements leads to a more accurate geometrical representation of the domain's boundaries at the cost of introducing locally constricted elements into the mesh.

The penalty on the time step introduced by the use of locally refined elements is further compounded when considering that the size of the elements is proportional to the speed of sound within the CFL criterion. That is, increasing the speed of sound should ideally be compensated for by increasing the element size in order to preserve the time step. Thus, by introducing such a refinement region in an area with locally increased velocities leads to the compounding of these penalties given that the elements are reduced in size instead of being increased.

It is also worth noting that while the compute cost in terms of node hours for the refined rectilinear mesh at 2.0 EPW is high compared to those of most of the other meshes tested, the absolute compute time is considerably longer than those of all other simulations. For example, the simulation that had the second longest compute time was the cubed sphere mesh at 4.0 EPW (101.41 × 10⁶ elements and 46 400 time steps), which ran in 8 min and 45 s on 204 GPUs. By comparison the refined rectilinear mesh at 2.0 EPW (16.84 × 10⁶ elements and 197 540 time steps) took 31 min and 2 s to run on 32 GPUs.

While part of this difference in compute time is obviously due to the dramatically different number of GPUs used, the difference in compute resources stems from the disparity in temporal versus spatial parallelism that can be achieved for this type of time domain solver. The time stepping scheme used in this implementation inherently requires each time step to be computed sequentially, meaning that it is extremely difficult to achieve temporal parallelism without using a different temporal discretization approach. By comparison, spatial parallelism can be achieved much more effectively by partitioning the finite-element mesh and distributing this across multiple GPUs. If one were to increase the number of GPUs used for the refined rectilinear mesh, it would have significantly diminishing computational returns given that each GPU would be under-utilized while also increasing the amount of inter-node communication. Thus, the compute time may be reduced somewhat, but the number of node hours used would increase significantly.

While these computational penalties of using local refinements in this context are indeed severe, much of the staircasing artifacts which were apparent by visual inspection for the other rectilinear meshes have been largely mitigated. This is particularly evident within the focus of the wavefield given that the shape and amplitude between the focal points for the refined rectilinear and the cubed sphere meshes are very similar. However, some differences in terms of multiple scattering remain apparent and are, in turn, also reflected in the increased $ℓ2$ norm relative to the cubed sphere meshes.

In order to validate the modelling techniques performed with the SEM, the results obtained from the software k-wave (Treeby and Cox, 2010; Treeby , 2020) were compared to an identical simulation setup performed within salvus. k-wave is based on the k-space pseudospectral method and was selected for this comparison given its extensive use as a full-waveform modeling tool within the medical ultrasound community in addition to its previous validations to laboratory measurements (Martin , 2020). The axisymmetric implementation of k-wave is considered, given the exceptionally fine spatial discretization (i.e., high points per wavelength), which can be feasibly achieved using this formulation when compared to a full 3D simulation (Treeby , 2020).

An axisymmetric domain identical to that used for PH1-BM6 in the intercomparison study by Aubry (2022) is considered. A domain setup consisting of a layered spherical cap comprised of skin, trabecular bone, and cancellous bone is used, as seen in Fig. 6. The material properties outlined in Table III are used for each of the distinct media within Fig. 6.

Similar to the simulations performed in Sec. IV A, a continuous sine wave was emitted over the course of the 120 μs simulation to reach steady state. A bowl-type transducer with identical dimensions to that considered in Sec. IV A and a source frequency of 500 kHz was used for generating the signal.

A comparison of the results obtained with both the axisymmetric k-wave implementation and salvus can be seen in Fig. 7. The resulting pressure distributions between the two simulation tools are in very good agreement, with the relative $ℓ2$ and $ℓ∞$ norms being 4.9% and 9.0%, respectively.

In addition to the $ℓ2$ and $ℓ∞$ differences, which are computed for the entire pressure field, the characteristics at the foci are also in excellent agreement with each other. In particular, the maximum amplitude of the focus deviates by approximately 3.4%, while the location of the focus is identical between both numerical schemes (accurate to within the spatial discretization of the regular grid). Finally, the difference in the size of the focus (full width at half maximum) is minimal, with differences in the x- and y-directions being 0.0021 mm and 0.0143 mm, respectively.

Despite solving the wave equations using significantly different numerical approaches, the pressure distributions depicted in Figs. 7(a) and 7(b) demonstrate that the numerical results converge to a common solution. This gives confidence that the results obtained using this SEM approach yield comparable results to other already well-established tools used in transcranial ultrasound.

An adaptation of the multimodal imaging-based detailed anatomical (MIDA) model (Iacono , 2015) is considered for demonstrating the use of the previously outlined modeling procedures on a domain with high geometric complexity. The MIDA numerical phantom is a segmented model of the human head which was constructed using several different MRI modalities. Representative material properties from the literature were assigned to each of the discrete anatomical features within the phantom (Duck, 1990; Pichardo , 2017; Webb , 2021). A summary of the material properties used within the adapted MIDA phantom can be seen in Table IV.

A bowl-type transducer identical to that from the Shepp–Logan sensitivity analysis in Sec. IV A is used to model the focusing within this tissue model. All source characteristics such as the source time function and reference frequency remain unchanged from the aforementioned setup.

The conforming hexahedral mesh was constructed using the procedure outlined in Sec. III A. The inner and outer tables of the skull, which are distinctly resolved within the MIDA model, are merged to form a single manifold surface. Using this closed surface representation as a basis, a conforming hexahedral mesh is generated using the sculpting algorithm such that only the boundaries between the fluid regions and the skull are explicitly meshed. This hexahedral mesh can be seen in Fig. 8.

Similar to the elliptical inclusions that are incorporated into the modified Shepp–Logan phantom from Sec. IV A, the soft tissues and the porous diploë layer of the skull are interpolated onto the hexahedral mesh as a post-processing step. The hexahedral mesh contains approximately 8.54 × 10⁶ elements, which, given the fourth-order accurate representation in space, corresponds to approximately 1.07 × 10⁹ nodal values for each unique material property of $vp(x), vs(x)$⁠, and $ρ(x)$⁠, in addition to the attenuation parameters of $αp(x)$ and $αs(x)$⁠.

It is important to reiterate the ease with which these meshes can be constructed despite the high level of geometric complexity, which is made evident by Fig. 8. The level of complexity of constructing these types of fully unstructured hexahedral meshes certainly persists, although the advent of semi-automated meshing strategies such as the parallel hexahedral meshing from volume fractions algorithm introduced by Owen (2014) has made constructing these types of computational domains much more feasible than in the past.

Obtaining such spectral-element meshes for applications of comparable complexity is quoted in the literature as taking on the order of week to months to construct (Igel, 2017). In contrast, the mesh illustrated in Fig. 8 was constructed on the order of minutes on a single workstation with no manual user intervention required. This, of course, implicitly assumes that a surface representation of the features which are to be meshed is available, the construction of which can be non-trivial. Conducting a comprehensive quantitative study of the compute times associated with constructing these types of hexahedral meshes is outside of the scope of this work. However, this provides a qualitative impression as to the power and flexibility of these types of semi-automatic meshing strategies.

Two separate simulations were conducted: (1) treating all tissues as lossless media and (2) modeling both the viscoacoustic and viscoacoustic losses within the domain. The resulting pressure distributions within the domain using the aforementioned bowl-type transducer can be seen in Fig. 9. Furthermore, simulation statistics consisting of the maximum pressure, full width at half maximum, and focal volume computed at −6 dB are shown in Table V.

The differences in terms of the maximum pressure distribution are immediately evident upon inspection of Fig. 9. In particular, the incorporation of attenuation leads to (1) the dissipation of many low-amplitude multiply scattered waves and (2) the maximum focal pressure being reduced.

The maximum pressures at the focus point for the cases with and without attenuation are 640.5 kPa and 516.5 kPa, respectively; this corresponds to a difference of approximately 20%. This difference in peak pressure is comparable to that reported by Treeby and Cox (2014), where similar comparisons between elastic and viscoelastic formulation were performed.

The focal volume computed at −6 dB is similar in shape between the non-attenuating versus attenuating cases. However, the focus becomes slightly larger with the inclusion of attenuation, as indicated by the approximately 5% increase in focal volume; the change in the focal volume is notably lower than the change in the maximum pressure. Similar to the relatively small change in focal volume, the full width at half maximum is quite similar between the attenuative and non-attenuative cases. This suggests that while the maximum pressure changes significantly between the attenuative and non-attenuative cases, the overall shape of the focus is preserved. The position of the focal volume at −6 dB can be seen in Fig. 10.

The use of the SEM has been explored for solving the wave equation related to ultrasound applications. The coupled viscoacoustic-viscoelastic formulation of the wave equation for spectral-elements is presented, and a number of numerical examples are provided in order to demonstrate the efficacy and flexibility of this technique.

A sensitivity analysis using an adaptation of the Shepp–Logan phantom highlights the importance of correctly meshing the interfaces between the soft tissue and bone within the SEM when performing these types of ultrasound simulations. Staircasing artifacts observed at material interfaces, which are inherent to nearly any numerical wave solver, are particularly pronounced at interfaces where strong material contrasts are present. While these effects can be partially mitigated by applying element refinements in parts of the mesh close to these distinct material interfaces, the significant penalty in the global time step does not lend to the scalability of such refinement techniques.

These modeling approaches have been validated against the simulation outputs from the axisymmetric k-wave code, a widely used tool within the medical ultrasound community. The simulation results through the use of the SEM are in excellent agreement with the pressure distributions obtained using k-wave, with $ℓ2$ and $ℓ∞$ differences of 4.9% and 9.0%, respectively. Furthermore, the focus' maximum amplitude, position, and dimensions are in excellent agreement between the two solvers. This strong agreement between these two solutions gives confidence in the quality of the full-waveform modeling results obtained using the SEM, given how fundamentally different the numerical methods of the SEM and the pseudospectral approach are from one another.

Finally, a numerical example using an adaptation of the MIDA is presented in order to demonstrate the SEM's ability to model the propagation of ultrasound waves in domains with a high degree of geometric complexity. The forward modeling results for cases with and without attenuation are compared to highlight the differences in the pressure distributions between these two cases.

While the numerical examples in this study involve transcranial FUS, the same modeling principles are also relevant for transcranial USCT. Being able to accurately solve the wave equation is imperative within frameworks such as FWI, and thus, these forward modeling techniques are directly transferrable to such waveform-based inversion approaches. Solving such an inverse problem is outside of the scope of this study, although this may be the subject of future work.

The techniques demonstrated throughout this numerical study can also be applied to modeling the propagation of mechanical waves in both other medical contexts as well as within other disciplines. The same modeling strategies can be applied to other parts of the body where soft tissue-bone interactions are relevant, such as for joints. Furthermore, these techniques can be applied across other spatial scales within non-destructive testing or in exploration geophysics for use-cases where modeling fluid-solid coupling is also of importance. Examples of such applications could include performing corrosion assessments of fluid-filled pipelines or modeling the attenuative effects of the ocean within offshore exploration.

While the modeling approaches demonstrated in this work possess many attractive properties, several limitations exist within the current modeling framework. First, the compute cost required to solve the wave equation may be significant, depending on the source frequencies which are to be resolved within the simulation. This high compute cost may be exacerbated in cases where meshes of poor quality are used given the potential penalty in the global time step for meshes with locally constricted elements.

In addition to the potentially significant computational cost associated with solving the wave equation, constructing these unstructured hexahedral meshes is significantly more involved when compared to a regularly gridded discretization, which is trivial to construct by comparison. Generating these fully unstructured hexahedral meshes remains an open problem in computational geometry and generally requires some degree of manual user intervention in order to ensure that the mesh is of adequate quality. Furthermore, one needs to have sufficient a priori knowledge of the shape of the skull in order to be able to construct such a spectral-element mesh.

In this work, the use of the SEM is explored for use within transcranial ultrasound modeling. The two primary objectives of this work are (1) to present a coupled viscoacoustic-viscoelastic formulation of the wave equation within the SEM for modeling dissipative effects in coupled fluid-solid media and (2) to demonstrate an end-to-end workflow for performing spectral-element simulations in transcranial ultrasound. Several numerical examples are presented in order to demonstrate the flexibility of this modeling approach within transcranial ultrasound. The importance of precisely meshing the soft tissue-skull interfaces in order to eliminate staircasing effects is shown using an adaptation of the Shepp–Logan phantom. These modeling results obtained using the SEM are then validated against k-wave, given the exceptionally high points per wavelength that can be feasibly achieved using the axisymmetric implementation within k-wave. As k-wave is among the most widely used numerical tools in the medical ultrasound community, the excellent agreement between the axisymmetric k-wave and SEM solutions gives confidence in the results obtained using the SEM. Finally, an adaptation of the MIDA phantom is considered in order to demonstrate the ability for these approaches to effectively deal with domains with a high degree of geometric complexity.

This research was supported by the Centro Svizzero di Calcolo Scientifica [Swiss National Supercomputing Centre (CSCS)] under Project ID Nos. s1238 and c20 on the Piz Daint and Alps supercomputers, respectively. The authors would also like to thank Scott Keating for the many constructive discussions throughout the manuscript writing process.

The authors have no conflicts to disclose.

Raw data were generated at the Swiss National Supercomputing Centre (CSCS) large scale facility. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.