## 3 Results  
The accuracy of the proposed computational framework needs to be verified first before performing parametric studies to investigate the effects of the inclusion on the temperature distribution during ultrasound therapies. Theoretical solutions of the pressure fields and the numerically simulated temperature fields in a parallel bilayer system (similar to case b in Table 2), which is a degenerated form of Figure 1 under the plane incident wave from the ultrasound transducer are compared with the experiment data and analytical solutions in the literature [15; 17]. Details are shown in Section S3 of the Supplementary Data. Hereafter, a parametric study of the ultrasound diathermy is performed. \begin{table}
\begin{tabular}{l c c c c} \hline \hline
**Material Property** & **Tissue** & **Stainless Steel** & **HDPE** & **Bone** \\ \hline Density (kg/m\({}^{3}\)) & 1190 & 7710 & 960 & 1912 \\ Specific heat (J/kg \({}^{\circ}\)C) & 3431 & 502 & 2300 & 1313 \\ Speed of sound (m/s) & 1512 & 5790 & 2460 & 4400 \\ Thermal conductivity (W/m \({}^{\circ}\)C) & 0.6 & 16.27 & 0.442 & 0.32 \\ Attenuation coefficient (db/m) & 54 & 124 & 66 & 900 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Material properties of the soft tissue, inclusions, and bone [17]. The considered composite system, as shown in Figure 1, consists of two layers. The minimum thickness of the soft-tissue layer \(h_{1}\) was 15 mm, while the maximum thickness of the bone layer \(h_{2}\) was 15 mm. In addition, both layers have width, \(L\), of 80 mm. It is noted that the tissue-mimicking hydrogel phantoms are adopted to replace the real soft tissues in this study. The initial temperature of the composite system was 26 \({}^{\circ}\)C based on the initial ambient temperature considered for the ultrasound heating experiments. The temperature of the front surface of the composite system (near the ultrasound transducer) was held constant at 26 \({}^{\circ}\)C, while the other three outer surfaces were adiabatic. The area beneath the ultrasound probe was the heating region, with non-uniform heating implemented by the user subroutine DFLUX (see Equation (S21) of the Supplementary Data). A time step of 1.0 s was used in the thermal transient analyses. In addition, to fit the results of the in vitro experiments [15], the amplitude coefficient \(A^{int}\) (see Equation (S1) of the Supplementary Data) was calibrated to values of 1.9 \(\times\) 10\({}^{5}\) and 1.25 \(\times\) 10\({}^{5}\) for 1 MHz and 3 MHz ultrasound heating, respectively. To investigate the effects of different parameters, such as inclusion position \(S\), inclusion size \(D\), inclusion material \(M\), radius of the convex bone \(R\), and ultrasound operation frequency \(f\), on the temperature profiles, a parametric study was performed with ten FE models. The parameters for these ten models are listed in Table 2. The first two cases are the control experiments, in which the composite system is only composed of soft tissue in case a and both the tissue layer and bone layer are flat in case b. It should be noted that an inclusion size of zero denotes a lack of an inclusion and a convex bone with an infinite radius means that the bone is flat. Stainless steel (Steel) and high-density polyethylene (HDPE) were considered as the two inclusion materials. The material properties of the soft tissue, inclusion, and bone used in the FE model are listed in Table 1.  
Figure 3a-j present the pressure distributions for ten cases. In these figures, the purple dashed lines represent the interface between the soft tissue and the bone, while the purple dotted lines represent the circumference of the circular inclusions. Note that only the illuminated regions under the ultrasound probe of the composite systems are shown. The pressure values in each region are initially normalized with respect to the incident amplitude coefficient (i.e., \(A^{int}\) in Equation (S1) of the Supplementary Data).