Quantitative validation of ultrasound thermal strain imaging sensitivity using simulation and oil–gelatin phantom experiment

Ultrasound thermal strain imaging (US-TSI) exploits temperature-dependent changes in speed of sound (SoS) to monitor tissue temperature change and characterize tissue composition. However, quantitative relationships among composition, temperature, and thermal strain remain insufficiently characterized. We developed a quantitative framework using measured SoS–temperature relationships, Wood’s equation-based mixture modeling, Field II simulations, and experimental validation. Thermal strain was estimated using a slope-based displacement method. Thermal strain showed a linear dependence on temperature rise (simulation: R 2 = 0.975–0.996; experiment: R 2 = 0.868–0.968) and on oil mass fraction (simulation: R 2 = 0.998–0.999; experiment: R 2 = 0.925–0.993) within a temperature-rise range of 1 °C–4 °C. These results demonstrate strong agreement among theory, simulation, and experiment and support the use of US-TSI for quantitative tissue characterization.

When tissue temperature rises, local changes in the speed of sound (SoS) alter the time-of-flight of ultrasound waves, while tissue thermal expansion also contributes to temperature-induced echo shifts. If beamforming is performed with an assumed constant SoS, these combined effects produce apparent axial displacements in reconstructed images. The spatial gradient of this temperature-induced axial displacement is defined as thermal strain (Maass-Moreno and Damianou, 1996; Seo et al., 2011). Early studies identified the temperature dependence of ultrasound propagation and explored ultrasound-based thermometry (Gammell et al., 1979; Ueno et al., 1990). Subsequently, Maass-Moreno and Damianou (1996) provided an analytical framework that relates temperature-induced echo shifts to changes in SoS and thermal expansion, establishing a theoretical foundation for thermal strain imaging. Thermal strain carries useful information in biomedical applications in two primary ways so far. First, it provides a measure of local temperature change, enabling noninvasive thermometry for hyperthermia treatment monitoring and radiofrequency (RF) ablation guidance (Seip and Ebbini, 1995; Varghese et al., 2002; Arthur et al., 2005; Souchon et al., 2005), and more recently for real-time temperature feedback during focused ultrasound–induced mild hyperthermia (Lee et al., 2026). Second, because the temperature dependence of SoS differs between water-bearing and lipid-bearing tissues, the sign of thermal strain encodes tissue composition, which can be used for tissue characterization, more specifically, for lipid identification (Seo et al., 2011; Shi et al., 2005; Kim et al., 2008).

Early approaches of thermal strain for tissue characterization employed electromagnetic energy deposition to induce temperature rise. For example, microwave-induced thermal imaging (MITI) demonstrated that temperature-induced changes in SoS could be estimated from speckle-based displacement tracking, enabling differentiation between lipid-bearing and water-bearing tissues (Shi et al., 2003). However, electromagnetic heating cannot be focused at depth. In addition, it often creates non-uniform temperature distributions that complicate quantitative interpretation and limit clinical translation. This limitation motivated the development of ultrasound thermal strain imaging (US-TSI), in which ultrasound is used for both imaging and heating. Early implementations demonstrated that thermal strain could be both induced and measured using ultrasound arrays, including configurations with separate heating and imaging probes as well as systems using a single linear array for both functions (Kim et al., 2008; Huang et al., 2007). These studies established the feasibility of generating sufficient temperature rise (∼1 °C–2 °C) and measurable thermal strain using clinically compatible ultrasound systems. Subsequent developments have focused on improving temporal resolution and motion robustness, for example, through high-frame-rate imaging techniques such as plane-wave acquisition (Sheng et al., 2025; Chen et al., 2025). This approach enables controlled energy delivery within a defined focal region while maintaining co-registered imaging and tracking, which is more promising toward translation in future clinics.

For small temperature rises, thermal strain is approximately proportional to the fractional change in SoS, resulting in near-linear dependence on temperature in the range of approximately 10 °C near physiological temperature (Nasoni and Bowen, 1989; Miller et al., 2004; Seo et al., 2011). Within this temperature increase range, negative strain exhibits in water-bearing tissues and positive strain exhibits in lipid-rich tissues. Kim et al. (2008) demonstrated the technique’s potential using excised rabbit kidneys ex vivo, and subsequent in vivo studies validated it across multiple animal models: lipid detection in atherosclerotic plaques in rabbits and pigs (Mahmoud et al., 2013; Sheng et al., 2025), and hepatic fat assessment in obese mice (Khalid et al., 2021). More recently, Chen et al. (2025) developed a US-TSI probe capable of millisecond-level transient heating, achieving approximately 2 °C in 50 ms in a porcine model.

Despite these advances, most existing TSI studies rely on assumed or literature-based SoS–temperature relationships and simplified material representations. As a result, thermal strain has mostly been interpreted qualitatively for contrast, but the quantitative relationship between mixture of material composition, temperature-dependent SoS, and systematic measurements of thermal strain remains incompletely characterized. In medium such as lipid–water or oil–gelatin mixtures, effective SoS depends on the contribution of compressibility and density of each constituent, which can be described by classical mixture models such as Wood’s equation (Wood, 1930; Wilson and Roy, 2008). However, integrated studies that combine direct SoS and temperature measurements, density characterization, mixture modeling, theoretical strain derivation, and thermal strain imaging validation within a unified framework have not been explored.

In this study, we establish a quantitative framework to investigate the sensitivity of US-TSI using oil–gelatin phantoms. Oil-based materials are commonly used as surrogates for lipid content in tissue-mimicking phantoms (Endsley et al., 2025; Ohmae et al., 2018). Within the temperature range of 24 °C–40 °C, SoS as a function of temperature was measured for both pure oil and gelatin solutions, accompanied with independent density measurements. This specific temperature range was selected to consider both the room-temperature conditions (approximately 25 °C) employed in phantom experiments and the physiological temperatures (approximately 37 °C) relevant to in vivo applications. Higher temperatures were not investigated in the phantom studies, as the gelatin matrix transitions into a liquid state under such high-temperature conditions. Phantom experiments were first conducted under controlled heating conditions to systematically characterize the relationship between oil mass fraction, temperature rise, and measured thermal strain. Based on experimentally measured sound speed–temperature relationships of the pure gelatin and oil materials, a physics-based model was then developed using Wood’s mixture equation to predict effective SoS as a function of oil mass fraction and temperature. Based on the predicted SoS, we performed Field II simulations (Jensen, 1996) to generate ultrasound radiofrequency (RF) signals. The simulated RF signals were then beamformed to obtain B-mode images, followed by displacement tracking and thermal strain estimation as the axial gradient of the displacement. By comparing simulation results with experimental measurements, this study validates the theoretical and simulation model of US-TSI for lipid identification, which further quantitatively provides insight of US-TSI sensitivity for detecting lipid contents.

Quantitative framework and model formulation

US-TSI is fundamentally governed by three key variables: oil mass fraction, temperature rise, and thermal strain. Therefore, we established a quantitative framework to find the quantitative relationships among these variables.

Thermal strain vs. temperature rise at fixed oil mass fraction

Thermal strain as a function of temperature rise was first analyzed at fixed oil mass fraction. According to thermal strain theory, thermal strain ε is approximately given by where β is the thermal expansion coefficient, λ represents the tissue characteristics and Δ T denotes temperature rise (Seo et al., 2011; Nasoni and Bowen, 1989). To assess whether neglecting thermal expansion was reasonable for the present materials, we compared the magnitude of the expansion-related term with the measured SoS-related term.

At 37 °C, the measured SoS–temperature relationships gave λ = 9.82 × 10 −4 /°C for the gelatin solution and λ = −1.42 × 10 −3 /°C for castor oil. Using the specific-volume data of Guignon et al. (2010), the castor-oil volumetric expansion coefficient was estimated as 7.24 × 10 −4 °C −1; assuming isotropic free expansion, β ≈ α/3 = 2.41 × 10 −4 /°C, approximately 17% of |λ|. For the highly hydrated gelatin matrix, a water-dominant linear coefficient of approximately 1.2 × 10 −4 /°C (Medina, 2012) corresponds to approximately 10%–12% of the measured gelatin SoS coefficient. Thus, the SoS-induced term remains dominant, although expansion is not zero. The castor-oil estimate represents a conservative free-liquid upper bound because oil droplets embedded in the solidified gelatin matrix are mechanically constrained. Because a directly matched effective expansion coefficient was unavailable for the solidified composite, the Field II models evaluate the SoS-induced component and are interpreted as an SoS-dominant approximation. Therefore, thermal strain (ε) can be approximated using Equation 2.

Accordingly, a linear model, as given in Equation 3, can be used to describe the relationship between thermal strain ε and temperature rise ΔT.

The slope a (strain/°C) quantifies the thermal strain sensitivity of the medium at a given oil mass fraction, while b represents the intercept of the linear fit, accounting for baseline offsets.

Thermal strain vs. oil mass fraction at fixed temperature rise

To relate thermal strain to SoS, thermal strain can be expressed in terms of the two sound-speed states before and after heating (c 0 and c 1), as shown in Equation 4.

This approximation assumes that heating is spatially uniform within the phantom, leading to a consistent change in SoS (c 0 → c 1) under a fixed temperature rise.

Then the thermal strain could be expressed as a function of: where ϕ denotes oil mass fraction. T 0 denotes the temperature before heating and T 1 denotes the temperature after heating. The effective SoS in the mixture (c mix)was modeled using Wood’s equation under the long-wavelength approximation, assuming the medium behaves as acoustically homogeneous (Wood, 1930; Wilson and Roy, 2008):

Here, ρ o i l and ρ g e l denote density of oil and gelatin, respectively. And c o i l and c g e l denote SoS in oil and gelatin, respectively. V is the oil volume fraction. Both density and SoS depend on temperature.

Since Wood’s equation is defined in terms of volume fraction, v was obtained from the mass fraction ϕ as:

Mixture density ρ mix was then calculated as: combining Equations 6–8, the mixture SoS in terms of oil mass fraction and temperature:

Although the resulting expression ε ϕ by substituting Equation 9 into Equation 5 is not strictly linear, the nonlinearity over the investigated range (0%–40% oil) is negligible. Therefore, at each fixed heating level, a linear regression was performed between strain and oil mass fraction according to Equation 10, yielding a slope k that quantifies strain–oil mass fraction sensitivity.

A larger | k | indicates stronger compositional contrast and improved sensitivity of oil mass fraction, while m represents the intercept of the linear fit.

Parameter measurement

To provide model inputs, density and SoS of pure components were measured independently.

Densities of the gelatin solution and castor oil were measured using a gravimetric method (three repeats per sample). The measured densities were 0.9250 ± 0.002 g / c m 3 for the gelatin solution and 0.9029 ± 0.003 g / c m 3 for castor oil, reported as mean ± standard deviation. These measured values were used to convert oil mass fraction to volume fraction and to calculate the mixture density in Wood’s equation.

SoS in pure oil and gelatin solutions was measured over a temperature range of 24 °C–40 °C, covering both room-temperature phantom conditions (∼25 °C) and physiologically relevant temperatures (∼37 °C). Temperature was monitored using a calibrated Type-T thermocouple (MT29, Physitemp Instruments, LLC, Clifton, NJ, USA), and the mean temperature during acquisition was used for analysis. Two parallel suture lines with known physical separation (d physical) were embedded in the medium. The axial distance between their echo peak (d acoustic was extracted from B-mode images, and SoS was calculated using Equation 11.

Where c beamforming denotes the assumed SoS used during beamforming which is 1,540 m/s. For each temperature point, SoS was averaged across 10 ultrasound imaging frames. Linear regression was used to obtain SoS–temperature relationships for oil and gelatin.

Phantom preparation

Oil–gelatin phantoms (Figure 1) with oil mass fraction ranging from 0% to 40% were fabricated. First, gelatin powder (Strength 300, Type A, Sigma-Aldrich, St. Louis, MO, USA) was dissolved in distilled water at 55 °C under constant stirring to prepare a 5% by weight gelatin solution. A water-soluble surfactant, Tween 80 (Sigma-Aldrich, St. Louis, MO, USA) was then added to this solution. Separately, the oil phase was prepared by heating the castor oil (Thermo Fisher Scientific, Waltham, MA, USA) to 80 °C and mixing it with an oil-soluble surfactant, Span 80 (Sigma-Aldrich, St. Louis, MO, USA). The oil mixture was slowly added dropwise into the gelatin solution while stirring at 70 °C and 1,100 rpm to form a stable emulsion. A guar gum solution was subsequently added for increasing adhesion. After cooling to 38 °C, 1% by weight cellulose (Type 20, Sigma-Aldrich, St. Louis, MO, USA) was added as a scattering agent. The mixture was stirred in a mixer and poured into small containers (<200 mL), which were rapidly cooled in an ice–water bath. After cooling, the phantom was covered with pure gelatin to prevent oil leakage. Optical microscopy of the emulsified mixture was performed after mixing and before gelation using a stage-micrometer calibration of 1 division = 0.01 mm. Representative, clearly distinguishable droplets in the selected field were on the order of tens of micrometers, approximately 30–80 μm. At the 7.8 MHz imaging frequency, the corresponding acoustic wavelength in the phantom is approximately 0.19–0.20 mm, giving representative diameter-to-wavelength ratios of approximately 0.15–0.4. At the 3.47 MHz heating frequency, the corresponding wavelength is approximately 0.43–0.44 mm, giving representative diameter-to-wavelength ratios of approximately 0.07–0.19. Thus, the long-wavelength approximation is better satisfied at the heating frequency than at the imaging frequency. These observations are consistent with the use of Wood’s equation as a first-order bulk effective-medium approximation, but they do not demonstrate strict microscopic homogeneity. The imaged field showed a broad apparent droplet-size distribution and localized clustering, and the long-wavelength condition was not satisfied for every inclusion. Because the microscopy sample was imaged as a liquid mixture under a coverslip before gelation, the measurements were used only to estimate the droplet-size order of magnitude.

Experimental workflow and signal processing

An overview of the experimental workflow is illustrated in Figure 1. Phantom experiments were performed on oil–gelatin phantoms with different oil mass fraction using a Vantage 256 system (Verasonics Inc., WA, USA) with an external power supply (QPX600DP, Aim and Thurlby Thandar Instruments, UK) connected through a UTA-408 adapter. In this study, a single focused heating beam was used. Different heating levels were applied by varying the transmit voltage from 10 V to 40 V. Each heating event used a 120-cycle pulse repeated over 50 ms with a 50% duty cycle. Ultrasound imaging was performed before and after each heating period. For imaging, a 5-angle compounded plane wave sequence (−10°, −5°, 0°, 5°, 10°) was used with a center frequency of 7.8 MHz and a frame rate of 500 Hz. Each imaging period lasted 1.2 s, providing enough frames for stable displacement estimation. The extended acquisition allows selection of frames with consistent signal quality and reduces the influence of transient effects, including shear wave motion.

A thermocouple was inserted into the phantom at approximately 25 mm in the axial direction and 0 mm in the lateral direction, as shown in Figure 4a to measure the temperature. A pre-heating reference frame (frame 400) was selected, and axial displacement was computed between this reference and subsequent post-heating frames using an axial kernel corresponding to 0.5 mm in depth through Loupas tracking (Pinton et al., 2006). Displacement fields were stored as three-dimensional maps (axial × lateral × time). To suppress isolated outliers while preserving displacement structure, a 3D median filter was applied. Displacement profiles were then extracted along several scanlines centered at the thermocouple lateral position, averaged across neighboring five scanlines to suppress speckle noise and tracking variability. A three-scanline average was additionally evaluated as a processing-sensitivity analysis. For previous in vivo studies, axial thermal strain was obtained by spatial differentiation of displacement using a two-dimensional Savitzky–Golay (Luo et al., 2005) derivative method (filter based). So, we also applied this method for visualization as shown in Figure 4. However, this approach can be affected by filter settings and smoothing, which may lead to underestimation in the results. In this study, to enable a more controlled and direct comparison between simulation and experiment, we used a simple linear fitting of displacement within a selected region to compute strain. This slope-based method avoids the influence of filter tuning and provides a more consistent estimation under homogeneous phantom conditions.

A direct numerical comparison between the slope-based and Savitzky–Golay estimators was performed using the 10% oil datasets. The two methods were compared over both the matched 21–26.5 mm ROI and a localized 25–26 mm ROI near the thermocouple. Based on this comparison, the slope-based estimator was used for all quantitative analyses, whereas the Savitzky–Golay method was retained for spatial strain-map visualization.

Independent repeatability was additionally evaluated using newly prepared 0% and 40% oil endpoint phantoms. For each composition, three independent measurements were performed at low and high heating voltages (10 and 40 V). Temperature rise and five-scanline slope-based strain were summarized as mean ± standard deviation. These supplementary measurements were conducted separately from the original dataset used in Figure 6 and were intended to assess endpoint repeatability rather than reproduce the exact temperature–strain points in that figure.

Temperature data processing and extraction of ΔT

Temperature time series were recorded using a thermocouple for each heating voltage. To avoid artifacts immediately following the heating pulse (e.g., electrical noise and rapid transient variations), data within this time window were excluded. Instead, the post-heating cooling curve was fitted using a reduced-form heat diffusion model, providing a smooth temperature estimate and enabling extrapolation back to the peak temperature. The temperature rise (ΔT) was then defined relative to the baseline temperature measured prior to heating. Post-heating cooling behavior was modeled using a reduced-form solution of the two-dimensional heat diffusion equation. For an initial localized Gaussian temperature distribution, the analytical solution at the heating center predicts the decay described by Equation 12. where T ∞ represents the asymptotic equilibrium temperature, A is the initial temperature elevation at t = t 0, and τ = a 2 / 4 k is the effective thermal time constant determined by the initial spatial width a and thermal diffusivity k. The parameters T ∞, A, and τ were estimated via nonlinear least-squares regression.

The fitted peak temperature at t = t 0 was obtained from the model, and the temperature rise was defined according to Equation 13.

To validate the cooling-model-based temperature estimate, repeated thermocouple measurements were acquired at four heating voltages (15, 25, 30, and 35 V; n = 3 per condition). The directly measured temperature at t 0 was calculated by averaging the thermocouple measurements within ±0.1 s of t 0. The fitted temperature was obtained by evaluating the two-dimensional cooling model at the same time point. The sensitivity of the estimated Δ T rise to the fitting duration was additionally assessed using cooling-fit endpoints of 40, 50, and 60 s, with the 60-s fit used as the reference.

Simulation workflow

To investigate the effect of SoS variations predicted by the Wood’s equation, we verified through Field II simulation. An overview of the simulation workflow used in this study is illustrated in Figure 2. The thermal strain is simulated under different selected oil mass fraction (ϕ) matching experiment condition. For each simulation, RF echo signals were generated using the same transducer parameters as the experimental imaging setup and subsequently processed with the same displacement tracking and strain estimation pipeline.

A static scatterer phantom (40 mm × 30 mm, 10,000 scatterers) was generated with a fixed random seed to ensure identical scatter positions for all simulations. Two SoS states were defined for each heating condition: baseline (c 0) and temperature-rise (c 1), computed from measured component c(T)relationships and Wood-based mixture modeling. Scatterer positions were unchanged between states; only SoS was varied. Beamforming was performed assuming a constant reference SoS of 1,540 m/s, introducing time-of-flight mismatch and producing apparent displacement. RF data were simulated using a 192-element linear array (center frequency 7.8 MHz, element pitch 0.2 mm) with five plane-wave steering angles within ±10°. Delay-and-sum beamforming and IQ demodulation were applied. Axial displacement between one pre-heating frame and one post-heating frame was estimated using a phase-based Loupas algorithm (Pinton et al., 2006) with an axial kernel corresponding to 0.5 mm. Median filtering was applied to suppress outliers. To obtain thermal strain in the simulations, we used a slope-based displacement method same as the experimental processing. Simulations were repeated for heating conditions corresponding to experimentally measured temperature rises. For each heating condition, the SoS pair c 0, c 1 was updated accordingly, while all other parameters were held constant. The scatterer distribution was static and identical between pre- and post-heating states; no inter-state motion, mechanical vibration, convective flow, or decorrelation noise was introduced. The Field II results therefore represent an idealized upper-bound tracking condition intended to evaluate the composition- and temperature-dependent response, rather than the lower boundary of experimental performance.

Temperature-dependent SoS in gelatin and oil solution

Figure 3 summarizes the measured SoS–temperature relationships for the gelatin solution and oil over 24–40 °C. The fitted slopes show opposite signs, with a positive slope for the gelatin solution and a negative slope for oil (R 2 = 0.93 and 0.82, respectively). These linear fits provide the temperature-dependent relationships c g e l T and c o i l T used in the Wood’s equation and SoS simulations.

Linear regression of the measured castor-oil data yielded c (T) = −2.1 T + 1,557.0, with a slope of −2.1 m / s / °C (95% CI: − 2.9 to −1.3 m / s / °C). Joint propagation of the slope and intercept uncertainty through Wood’s mixture model preserved the monotonic composition-dependent strain trend. At T 0 = 25 °C and Δ T = 2 °C, the nominal strains for 0%, 10%, 20%, 30%, and 40% oil were approximately −0.19%, −0.14%, −0.09%, −0.04%, and 0.01%, respectively. For the 40% oil condition, the propagated 95% interval was approximately −0.03% to 0.04% and crossed zero, indicating uncertainty in the precise high-oil zero-crossing location. A boundary calculation over 37 °C–40 °C yielded a maximum absolute strain uncertainty of approximately 0.06%, supporting cautious interpretation near the upper measured temperature boundary.

Validation of cooling-model-based temperature estimation

The cooling-model-based temperature estimation was validated using 12 repeated thermocouple measurements acquired at four heating voltages as shown in Figure 4. The directly measured and model-fitted temperature rises at t 0 showed a strong linear association, with R 2 = 1.00, an RMSE of 0.24 °C, and a mean bias of −0.23 °C. The negative bias indicates a small systematic underestimation by the fitted model, while the voltage-dependent trend was preserved. Using the 60-s fit as the reference, shortening the fitting endpoint to 40 and 50 s produced mean absolute differences in ΔT (t 0) of approximately 0.041 °C and 0.017 °C, respectively. These differences were small relative to the measured temperature rises of approximately 0.7 °C–2.9 °C, indicating limited sensitivity to the selected cooling-fit duration.

Representative B-mode, displacement, and strain maps for different oil mass fraction

Figure 5 shows representative B-mode, displacement, and strain maps for 0% oil phantom. The 0% case exhibits strong, spatially coherent displacement and pronounced axial strain gradient.

Axial displacement profiles used for slope-based strain estimation

Figure 6 shows axial displacement profiles due to temperature changes induced SoS change, which are further used for thermal strain estimation. Each column corresponds to a different heating condition (increasing transmit power/temperature rise from power 1–7), while each row represents a different oil mass fraction (0%–40%). Displacement varies approximately linearly with depth across all conditions. A baseline offset between simulation and experiment is observed, this is because, in simulation, the SoS is changed across the entire field, so displacement accumulates along the full depth. In contrast, in experiments, heating is localized and only starts around 21 mm depth. This difference contributes to the offset between the simulated and experimental displacement profiles. Because thermal strain was estimated from the displacement gradient, a constant profile offset does not directly affect the slope estimate. However, the spatially nonuniform experimental temperature field may still contribute to residual differences in the fitted strain.

Thermal strain varies approximately linearly with ΔT in both simulation and experiment

Figure 7 shows the thermal strain as a function of temperature rise (ΔT) for different oil mass fractions (0%–40%), including both simulation (solid lines) and experimental measurements (dashed lines). For each oil mass fraction, thermal strain exhibits an approximately linear dependence on ΔT over the tested range. The linear fits yielded high coefficients of determination across all cases (simulation: R 2 = 0.975–0.996 and p < 0.001; experiment: R 2 = 0.868–0.968 and p < 0.005), validating the linear approximation over the investigated temperature range. Linear trends are observed consistently across both simulation and experiment, with strong agreement in slope for all oil mass fractions. Across all 35 matched simulation–experiment conditions, the overall RMSE was 0.016% strain, the MAE was 0.013% strain, and the mean simulation-minus-experiment bias was −0.001% strain. The RMSE corresponded to 7.8% of the total experimental strain range. Simulated and experimental strain values were strongly correlated (r = 0.95, R 2 = 0.89), indicating strong overall quantitative agreement with minimal systematic bias. The slope of strain versus ΔT (strain/°C), decreases monotonically in magnitude as oil mass fraction increases, transitioning from strongly negative values at 0% oil to near-zero or slightly positive values at 40% oil.

Additional processing-sensitivity analyses were performed by comparing strain estimates obtained using averages of three and five neighboring scanlines centered on the thermocouple position. Reducing the averaging width from five to three scanlines changed the estimated strain magnitude by approximately 4.2%–11.7% across the oil groups, while preserving the overall temperature- and composition-dependent trends. A direct comparison using the 10% oil datasets showed strong trend agreement between the slope-based and Savitzky–Golay estimates for both the matched 21–26.5 mm ROI (r = 0.98, R 2 = 0.96) and the localized 25–26 mm ROI near the thermocouple (r = 0.97, R 2 = 0.95). The Savitzky–Golay estimates generally exhibited smaller negative magnitudes, indicating that the two estimators were consistent in trend but were not numerically interchangeable. Across seven independent pre-heating acquisitions of the 10% oil phantom, the pooled 95th percentile of the absolute baseline strain was approximately 0.01%. Acquisition-specific 95th-percentile and 3 σ thresholds were of the same order, approximately 0.7–1.4 × 10 − 4.

Independent endpoint repeatability was evaluated from three measurements of newly prepared 0% and 40% oil phantoms. For the 0% oil control, the temperature rise and five-scanline strain were 0.1691 ± 0.0694 C ∘ and − 5.94 ± 4.99 × 10 − 5, respectively, at 10 V, and 1.9001 ± 0.0124 C ∘ and − 2.507 ± 0.058 × 10 − 3, respectively, at 40 V. For the 40% oil phantom, the corresponding values were 0.3039 ± 0.0118 C ∘ and − 1.118 ± 0.214 × 10 − 4 at 10 V, and 3.1017 ± 0.0223 C ∘ and 9.304 ± 0.398 × 10 − 5 at 40 V. The low-voltage measurements showed greater relative variability because the induced strain was close to the empirical baseline threshold, whereas the high-voltage endpoint measurements exhibited smaller standard deviations relative to their mean values.

Thermal strain varies approximately linearly with oil mass fraction in both simulation and experiment

Figure 8 shows the thermal strain as a function of oil mass fraction for four fixed temperature rises (ΔT = 1 °C–4 °C), including both simulation predictions and experimental measurements. For each heating level, thermal strain exhibits a strong linear dependence on oil mass fraction, with high coefficients of determination (simulation: R 2 = 0.998–0.999, p < 0.001; experiment: R 2 = 0.925–0.993, p < 0.01). The strain becomes progressively less negative with increasing oil mass fraction, approaching near-zero values around 35% oil mass fraction. Overall, the experimental results follow the same trend as the simulations, with slight deviations in magnitude.

The results of this study demonstrate that thermal strain in oil–gelatin phantoms can be quantitatively explained by temperature-induced change in SoS. By combining experimentally measured SoS–temperature relationships with mixture modeling based on Wood’s equation, the proposed framework enables prediction of thermal strain as a function of both temperature rise and oil mass fraction.

At 37 °C, the fitted SoS values of the gelatin solution and castor oil were 1,528.0 and 1,479.3 m/s, respectively. The predictions of Wood’s equation for the 10%, 20%, 30%, and 40% oil–gelatin mixtures were 1,522.6, 1,517.4, 1,512.3, and 1,507.3 m/s, respectively. As summarized in Table 1, increasing oil mass fraction shifted the modeled SoS from the water-rich gelatin value toward the range reported for lipid-rich tissues and tissue models. Published adipose-tissue measurements include a SoS of 1461.77 ± 26.11 m/s, while breast-tissue models use values spanning approximately 1,480–1,528 m/s (Fernandes et al., 2025; Simson et al., 2024; Bamber, 1983). Local liver SoS has also been reported to decrease with increasing hepatic fat content (Telichko et al., 2022). These comparisons support the biological relevance of the observed lipid-dependent acoustic trend. However, the present phantoms should be interpreted as controlled two-component models rather than exact replicas of atherosclerotic plaque, fatty liver, or other biological tissues, which additionally contain fibrous tissue, calcification, cellular structures, and spatial heterogeneity.

The relationship between temperature rise (ΔT) and thermal strain, as shown in Figure 7, demonstrates that strain varies approximately linearly with ΔT for each oil mass fraction. This observation is consistent with the established theoretical framework of thermal strain in ultrasound, where the apparent strain is governed by the temperature dependence of SoS. Specifically, under small temperature change, thermal strain can be approximated as a first-order function of temperature through the coefficient λ shown in Equation 1, which characterizes the temperature dependence of SoS (Seo et al., 2011). In this formulation, water-bearing tissues exhibit positive λ values, while lipid-bearing tissues exhibit negative λ values, leading to opposite strain responses under heating. For a temperature rise on the order of 1 °C, the expected thermal strain magnitude of water-based tissue is on the order of 10 –3 (≈0.1%), consistent with the magnitude (0.094%) of 0% oil mass fraction observed in this study. The approximately linear dependence of strain on ΔT observed in both simulation and experiment therefore agrees with this first-order approximation, confirming that temperature-induced SoS variation is the dominant mechanism governing thermal strain generation under mild heating conditions.

A key finding of this study is the strong linear dependence of thermal strain on oil mass fraction under fixed temperature conditions, as shown in Figure 8. For each ΔT, the strain–oil mass fraction relationship exhibits high linearity in simulation (R 2 > 0.99, indicating that oil mass fraction-dependent contrast can be quantitatively described within the investigated range (0%–40% oil mass fraction). The slope of each linear fit represents the compositional sensitivity of the system, defined as the change in strain per unit change in oil mass fraction. This definition follows classical measurement system theory, where sensitivity is defined as the ratio of output variation to input variation (Doebelin, 1990). The observed increase in slope magnitude with ΔT indicates that thermal excitation enhances sensitivity without altering the underlying linear relationship.

The sensitivity of measuring oil mass fraction using US-TSI is fundamentally governed by the signal-to-threshold relationship. As illustrated in Figure 8, for a fixed temperature rise, the thermal strain varies approximately linearly with oil mass fraction. In the conceptual schematics shown in Figure 9, the points where the strain curve reaches the detection threshold define the limits of detectable oil mass fraction (O i l min a n d O i l max ⁡. Increasing the temperature rise effectively increases the slope of the strain–oil mass fraction relationship, thereby increasing the strain magnitude for a given oil mass fraction. As a result, the intersections between the strain curve and the detection threshold (−| ε | min to +| ε | min) shift toward more extreme oil mass fractions, leading to an expanded detectable range, as indicated by the arrows. In the present system, the empirical strain detection threshold was estimated from pre-heating baseline measurements as approximately 1 × 10 − 4, or 0.01 % strain. Low- and intermediate-oil compositions generally produced strain magnitudes above this threshold, whereas values near the high-oil zero-crossing approached the baseline-noise and model-uncertainty range. Therefore, the method is more reliable for detecting composition-dependent strain differences than for estimating an exact oil fraction near the zero-crossing. A higher temperature rise (ΔT1), compared to a lower temperature rise (ΔT2), improves oil mass fraction detectability by increasing the signal relative to the detection threshold. At a temperature rise of approximately 2 ∘ C, the detectable oil-mass-fraction range should therefore be interpreted using the present system-specific threshold of 0.01 % strain rather than a universal criterion. The exact detectable bounds depend on both the strain–composition relationship and the experimental noise level.

The underlying physical mechanism is the opposite temperature dependence of SoS in gelatin and oil. As temperature increases, SoS in gelatin increases while that in oil decreases, resulting in divergence between the two components. This difference amplifies propagation-time mismatch under constant-speed beamforming, leading to increased apparent displacement and strain. This mechanism explains both the monotonic variation of strain with oil mass fraction and the amplification of compositional contrast with increasing temperature rise.

The close agreement between simulation and experiment further supports the validity of the proposed model. Experimental data closely follows the simulated trends, indicating that temperature-dependent SoS variation is sufficient to explain the observed strain behavior under controlled conditions. The remaining differences between simulation and experiment likely reflect experimental variability that was not represented in the idealized model. Because the Field II simulation used an identical static scatterer distribution before and after the SoS change, it did not include motion, signal decorrelation, electronic noise, or scatterer perturbation. The simulated results should therefore be interpreted as an idealized upper-bound tracking condition rather than a prediction of the experimental noise floor or lower detection limit. In practice, motion- and heating-induced decorrelation can increase displacement-estimation error and reduce detectability. Future simulations should incorporate controlled motion, decorrelation, and noise perturbations to better characterize the lower boundary of experimental performance.

Additional analyses were performed to assess the robustness of the experimental strain estimates and to examine the comparatively lower agreement observed in the 10% and 20% oil groups. Using three rather than five neighboring scanlines produced only modest changes in most groups and did not alter the main temperature- or oil-mass-fraction-dependent trends. The 20% group showed greater sensitivity to the lateral averaging width, consistent with increased local spatial variability in the measured displacement field, although the present analysis could not identify a single underlying cause. The slope-based and Savitzky–Golay methods also showed strong agreement in trend, despite differences in absolute strain magnitude. These findings indicate that lateral sampling, displacement-tracking variability, and the spatial mismatch between the point thermocouple measurement and the ROI-averaged strain estimate contributed to uncertainty in the fitted strain magnitude and linearity, particularly for the 10% and 20% groups. Nevertheless, the expected temperature- and composition-dependent trends were preserved. The slope-based strain should therefore be interpreted as an average displacement gradient over the 21–26.5 mm ROI rather than a local strain value at the thermocouple depth.

To assess whether the composition-dependent acoustic trend is expected to persist near physiological temperature, a model-based analysis was performed at a baseline temperature of 37 °C using the measured component SoS–temperature relationships. Because these relationships were characterized up to approximately 40 °C, the predicted SoS-dependent trend is supported within this temperature range. However, the analysis does not account for possible softening or mechanical changes of the solidified gelatin phantom during additional heating. The 37 °C results are therefore interpreted as predictions of the SoS-dependent sensitivity trend rather than a validation of the complete phantom response. The present framework was characterized only over the 24 °C–40 °C component-property range and under mild heating conditions. Extension above 40 °C, including thermal-therapy or ablation regimes, will require additional characterization of nonlinear SoS behavior, thermal expansion, and temperature-dependent mechanical or structural changes.

Several limitations should be considered. First, the oil–gelatin phantom was modeled as an acoustically homogeneous effective medium under the Wood-limit approximation. Microscopy showed that the droplets were generally on the order of tens of micrometers, but the size distribution was broad and localized clustering was present; therefore, strict microscopic homogeneity and the long-wavelength condition were not satisfied for every inclusion. Direct SoS measurements of the cured oil–gelatin mixtures were not obtained; therefore, the Wood-model predictions were not independently validated against mixture-level SoS measurements. Second, the analytical and Field II models included only the temperature-dependent SoS contribution. Although the available estimates indicate that this term was dominant under the investigated conditions, thermal expansion was not strictly negligible. Third, temperature was measured locally near the heating focus, whereas strain represented an average displacement gradient over the 21–26.5 mm ROI, introducing systematic uncertainty in the absolute strain magnitude. The cooling-model validation also showed a small systematic underestimation of temperature rise. In addition, the empirical strain threshold of approximately 1 × 10 −4 was specific to the present imaging and processing configuration. Finally, the phantom and static Field II model did not reproduce the structural heterogeneity, perfusion, motion, viscoelasticity, attenuation variability, or signal decorrelation encountered in biological tissue. Accordingly, the present study should be interpreted as a controlled assessment of composition-dependent SoS sensitivity rather than a demonstration of in vivo diagnostic performance. Future work will focus on spatially resolved temperature estimation, more realistic heterogeneous and dynamic simulations, and validation in temperature-stable tissue-mimicking materials and in vivo models.

In this study, we developed a quantitative framework linking oil mass fraction, temperature rise, and thermal strain in US-TSI. Measured SoS–temperature relationships showed strong linearity (gelatin: R 2 = 0.93; oil: R 2 = 0.82). Thermal strain exhibited a linear dependence on temperature rise (simulation: R 2 = 0.975–0.996; experiment: R 2 = 0.868–0.968) and on oil mass fraction (simulation: R 2 = 0.998–0.999; experiment: R 2 = 0.925–0.993), consistent with theoretical predictions. The strain magnitude decreased with increasing oil mass fraction, transitioning from negative to near-zero or slightly positive values. These results demonstrate strong agreement between simulation and experiment and highlight the potential of US-TSI for quantitative tissue characterization.