Anatomical Asymmetry in Adipose Tissue Thickness Alters Heat Transfer During Superficial Thermotherapy: A Computational and Experimental Study

1. Introduction

Superficial thermotherapy is a non-invasive therapy for musculo-skeletal pain, sports injuries and circulatory disorders due to low cost, and ease of application [ref]. Thermotherapy causes vasodilation, increases local blood flow, improves metabolism and promotes tissue healing and pain relief. However, the amount of heat transferred to deeper tissues depends on complex interactions between tissue layers, thermal properties, exposure time and physiological heat dissipation mechanisms; thus the prediction of internal heat distributions is difficult [ref]. Peripheral Artery Disease (PAD), which in recent epidemiological studies has been called lower extremity artery disease (LEAD) [1], is a condition that is frequently treated with this type of therapy. In 2015 around 234 million people had PAD and more than 230 million adults worldwide are affected [2,3]. Furthermore, the global incidence of PAD increased by 72% between 1990 and 2019 [4]. Superficial thermotherapy is a common method for pain relief and local blood flow improvement in LEAD patients by means of controlled heat application to the skin surface within a temperature range of 38 °C to 45 °C. Within this therapeutic range heat increases local metabolism, enhances blood flow and promotes relief of pain and restoration of function [5]. However, long exposures to temperature above 45 °C may cause adverse effects like protein denaturation and thermal tissue damage [6]. Thus, precise control of both treatment temperature and exposure time is crucial for the safety and efficacy of thermotherapy, prompting the development of clinically optimized treatment protocols [7,8].

This is because direct measurement of temperature inside biological tissues is invasive and often impractical in a clinical environment. Therefore the thermal response to superficial thermo-therapy is often indirectly measured through empirical measurements of the surface temperature [6]. Therefore, computational modeling has become a useful means of investigating heat transfer within subcutaneous tissues. Finite Element Models (FEM) based on the Pennes bioheat equation are able to simulate transient temperature distributions taking into account heterogeneous tissue properties, blood perfusion, metabolic heat generation and realistic boundary conditions. They provide a mathematical framework for the evaluation of treatment parameters, optimization of therapeutic regimens, and prediction of tissue thermal responses under conditions which would be difficult or impossible to evaluate experimentally [9]. Although widely used, many computational bioheat models assume symmetric tissue structures or simplified anatomical geometries. This may lead to an underestimation of the impact of patient-specific anatomical heterogeneity on heat transfer.

Heat transfer in biological tissues has been extensively studied using both experimental and computational approaches. FEM have become the predominant tool for analyzing temperature propagation during thermotherapy because they enable controlled evaluation of tissue properties and treatment parameters. Nevertheless, many existing studies rely on simplified symmetrical geometries or homogeneous tissue representations, while only a limited number incorporate anatomically realistic models or experimental validation. More importantly, the influence of anatomical asymmetry associated with adipose tissue thickness on heat transfer has received little attention. Table 1 summarizes representative studies addressing computational and experimental modeling of heat transfer during thermotherapy.

As summarized in Table 1, superficial thermotherapy is affected by several physiological and anatomical variables such as body composition, tissue properties and the heating modality used. The thermal insulating effect of adipose tissue limits the depth of heat penetration and consequently the therapeutic effect. Similarly, computational studies have shown that FEM are a useful tool to predict temperature distributions under different treatment conditions and anatomical scenarios. However, there are still several limitations. Most mathematical models are based on simplified symmetrical geometries, homogeneous tissue distributions or single patient-specific anatomies not intended to evaluate the influence of anatomic variability. Thus, although the insulating effect of adipose tissue has been experimentally acknowledged, the influence of systematic changes in adipose tissue thickness, and the consequent anatomical asymmetry, on heat transfer during superficial thermotherapy has been limited. This limitation is particularly relevant since standardized thermotherapy protocols are usually prescribed without taking into account the patient morphology, although the adipose tissue thickness varies substantially among individuals and even between the contralateral anatomical regions. Therefore, the common assumption of symmetric tissue geometries can lead to inaccurate prediction of heat propagation and, consequently, variability in the delivered thermal control.

In this work, by using a combined computational and experimental approach, we aim to study how the anatomical asymmetry induced by the thickness of adipose tissue influences the heat propagation during superficial thermotherapy. An anatomically representative FEM of the lower limb was developed to simulate varying thicknesses of adipose tissue under clinically relevant heating conditions. The remainder of this paper is organized as follows. The computational model, the experimental protocol and the numerical methodology adopted to evaluate the heat transfer under different anatomic conditions are described in Section 2. Section 3 refers to computational and experimental results and their comparative analysis. Section 4 concludes with a discussion of the main findings, limitations of the proposed approach and directions for future research, and Section 5 summarizes the principal conclusions of this study.

2. Materials and Methods

This section describes the experimental and computational methodology to study the asymmetrical effect of adipose tissue thickness on heat transfer during superficial thermotherapy. Firstly, a FEM of the lower limb was reconstructed and anatomically representative and tissue-specific thermophysical properties was assigned. Following, a parametric analysis was performed to simulate the transient heat transfer under clinically relevant thermotherapy conditions with different thicknesses of adipose tissue and different hot-pack temperatures. Finally, the quantitative results were experimentally validated using ex vivo tissue under controlled heating conditions.

2.1. Anatomical Model Construction

The two-dimensional (2D) anatomical model was developed based on a representative of a cross-sectional illustration of the lower limb taken from a human anatomy atlas [15], with realistic morphological asymmetrical proportions of the main tissue layers. The model retains the multilayer organization of the lower limb and the irregular internal distribution of the modeled tissues (see Figure 1). The relative position of bone, muscle, adipose tissue and skin was defined according to anatomical references and geometrical configuration of the model, as summarized in Table 2.

The two-dimensional (2D) anatomical model was constructed based on a representative cross-sectional illustration of the lower limb derived from a human anatomy atlas [15], preserving realistic morphological proportions of the principal tissue layers. As shown in Figure 1, the model preserves the multilayer organization of the lower limb and the irregular internal distribution of the modeled tissues. The relative arrangement of bone, muscle, adipose tissue, and skin was defined according to anatomical references and the geometrical configuration of the model, as summarized in Table 2.

The asymmetrical 2D geometry was designed using SOLIDWORKS 2025 SP2.0 (Dassault Systèmes, Vélizy-Villacoublay, France). To maintain the morphology and relative spatial arrangement of each tissue layer, the boundaries of the primary anatomical structures including bone, muscle, adipose tissue, and skin were manually traced with polyline-based tools. For the FEM implementation, each tissue layer was defined as an independent geometrical domain, thus the tissue specific thermal properties could be assigned. The model was exported using the Drawing Exchange Format (DXF), a standard computer-aided design format that allows interoperability between CAD and FEM environments while preserving the geometric structures. After importing the DFX model into COMSOL Multiphysics®, the multilayer geometry allowed to heterogeneously assign thermal properties and to independently control each domain during the multiphysics simulations. The proposed modeling approach offers sufficient flexibility to systematically vary the thickness of adipose tissue, while keeping other anatomical structures unchanged, thereby enabling a controlled evaluation of the effect of anatomical asymmetry on heat transfer in superficial thermotherapy. Though the external shape of the model preserves an anatomically realistic lower limb morphology, the tissue organization can be characterized by the normalized radial distribution summarized in Table 2. These radial boundaries ranges give an equivalent description of the multilayer anatomical structure for computational purposes, and allow to reproduce the proposed model. Importantly, this representation is not based on concentric symmetry, but rather represents an anatomically-informed approximation where the tissue regions are defined according to the observed geometry, while maintaining their relative proportions within the computational domain.

2.2. Bioheat Transfer Model

The transient form of the Pennes bioheat equation was used to describe the heat transfer in the multilayer model considering the combined effects of thermal conduction, blood perfusion, and metabolic heat generation. As context, this formulation has been used in computational studies of thermoregulation, hyperthermia, thermotherapy and bi–medical heat transfer [6,16,17] due to its ability to describe the main mechanisms responsible for heat transport in perfused biological tissues. In our study, Equation 1 was used to study the effect of the variations in adipose tissue thickness during superficial thermotherapy.

$${\rho }_i{c}_i{\displaystyle \frac{\partial Ti}{\partial t}}=\nabla ·\left(k_i\nabla T_i\right)+w_{b_i}{\rho }_{b_i}{c}_{b_i}\left(T_{a_i}-T_i\right)+q_{{met}_i}$$

(1)

where $i$ denotes the tissue layer, $T_i$ is the local tissue temperature (°C), ${\rho }_i$ is the tissue density (kg/m³), $c_i$ is the specific heat capacity (J/(kg·K)), $k_i$ is the thermal conductivity (W/(m·K)), $w_{b_i}$ is the blood perfusion rate (kg/(m³·s)), $q_{me{t}_i}$ is the metabolic heat generation (W/m³), $T_{a_i}$ is the arterial blood temperature (°C), ${\rho }_{b_i}$ is the blood density (kg/m³), and $c_{b_i}$ is the specific heat capacity of blood (J/(kg·K)).

The implementation assumes that each anatomical layer behaves as a homogenous and isotropic medium with constant thermal properties during the simulation. The interfaces between the neighboring tissues were considered with thermal contact, ensuring the continuity of temperature and heat flux at the boundary between the border domains. Also, the temperature of arterial blood and the metabolic heat generation were assumed constant during the heating period, while the thermophysical properties were considered temperature independent.

The thermal properties assigned for the skin, adipose tissue and muscle domains are given in Table 3. The parameters, i.e. thermal conductivity, density, specific heat capacity and blood perfusion, were acquired from the study carried out by Flores-Cuautle et al. [18], ensuring that the heat transfer properties were physiologically appropriate for each tissue layer.

The bone domain was included to maintain a realistic spatial distribution of the lower limb geometry. However, no thermophysical properties have been assigned to this region since the superficial heating conditions considered in this work are not expected to induce physiologically relevant temperature variations in the cortical bone. Therefore, the bone was considered as an inactive domain in the thermal analysis so that the surrounding tissues preserved the anatomical structure without affecting the temperature field obtained significantly.

2.3. FEM and Simulation Conditions

Transient heat transfer analysis was employed to study the effect of temperature propagation through the asymmetrical multilayer model. The properties and boundary conditions were defined independently for each anatomical region. The independent variables studied were hot-pack temperature, treatment time and lower-limb morphology characterized by variations in adipose tissue thickness. This allowed us to assess the effect of asymmetry on temperature propagation in a controlled manner, keeping the thermal properties and boundary conditions the same for all simulated cases.

The hotpack was modeled as an external heat source applied over a localized region of the skin surface corresponding to the anatomical location of the vastus lateralis muscle as shown in Figure 2. The contact length was fixed at 25 cm, the dimensions of the heating device used in the experimental validation. Two clinically relevant hot-pack temperatures, 38 °C and 44 °C, were evaluated according to reported thermal safety thresholds for superficial thermotherapy [20].

The temperature of the heat source was evaluated at 38 °C and 44 °C, which remains invariant throughout the simulation interval. The condition is expressed mathematically as Equation 2 [19]

$$T\left(r,t\right)=T_{hp} for r\in {\Gamma }_{contact},t>0$$

(2)

where $T\left(r,t\right)$ is the temperature as a position function (r) and time (t), $T_{hp}$ is the constant temperature (38 °C and 44 °C) on the surface, ${\Gamma }_{contact}$ is the specific contour where the contact occurs. The condition $t>0$ indicates that it applies from the moment the process begins.

Initially, all tissue domains were set to the baseline temperature of 37 °C and the treatment time was limited to 15 min, which is in line with conventional superficial thermotherapy protocols. The external limb perimeter was altered based on standardized anthropometric measures reported in the National Institute of Standards and Technology (NIST) [16] to accommodate anatomical variability. The analysis considered two representative morphologies, related to Medium (M) and Extra Large (XL) body sizes, with lower-limb perimeters of 52.4 cm and 62.3 cm, respectively. The adipose tissue layer was scaled to different sizes during the scaling procedure, while the remaining anatomical structures and material properties were kept intact. This made it possible to quantify the insulating impact of adipose tissue thickness on heat propagation. The parameters are listed in Table 4.

The parametric simulations were followed by a mesh convergence analysis to ensure numerical accuracy. The model was discretized with an unstructured triangular mesh with local refinement around the heating region and the tissue interfaces, where the highest thermal gradients were expected. The mesh convergence was evaluated from the stabilization of the temperature region. The criteria for mesh convergence were considered as less than 5% variation between two consecutive mesh refinements in certain control locations. The ten meshing configurations examined varied from 5,441 to 42,181 finite elements. Additional mesh refinement showed insignificant variations in the computed temperature field, indicating mesh-independent numerical solutions.

After spatial discretization, the governing equation can be expressed as a system of first-order differential equations, where the temporal evolution of the temperature field is governed by the balance between thermal capacity, conduction, and heat sources, as shown in Equation 3.

$$M{\displaystyle \frac{dT}{dt}}+KT=F$$

(3)

where $M$ represents the thermal capacity matrix, $K$ the conduction matrix, and $F$ the contribution of internal and external heat sources, including perfusion, metabolic heat, and the applied thermal load. The temporal discretization was performed using an implicit time integration scheme to ensure numerical stability under transient conditions.

The transient thermal problem was solved by the implicit time dependent finite element solver implemented in COMSOL Multiphysics® which ensures stable numerical integration for parabolic heat transfer problems over the whole simulation interval. All the simulations were performed on an Intel® Core™ i7-3630QM workstation with a clock speed of 2.40 GHz and 8 GB RAM.

2.4. Ex Vivo Experimental Validation

Ex vivo porcine tissue layers were used in a controlled setting to validate the simulation results. Porcine tissue was chosen because of its multilayer anatomical structure and thermophysical properties that are similar to human skin, subcutaneous fat and skeletal muscle, and is widely accepted as a surrogate for thermal studies involving superficial heat transfer [17]. The aim of the experimental study was to compare the transient temperature response predicted by the computational model with the direct temperature measurements taken within the muscle layer after superficial thermotherapy.

Before the experiments, the thickness of the skin, adipose tissue and muscle layer were measured and compared with the corresponding sizes of the computational model. The layer thicknesses are presented in Table 5 and comparing both configurations.

Temperature measurements were taken placing Type J thermocouples in the tissue to a depth of about 25 mm, ensuring that the sensing junction reached the muscle layer after reaching skin and adipose tissues. The thermocouples were connected to a National Instruments NI-9219 data acquisition unit and tissue temperature was recorded at 2-min intervals. Figure 3 presents the diagram of the experimental setup.

Prior to each experiment, the hot-pack was immersed in a thermostatically controlled water bath by 10 min to reach a temperature of 44 °C. The hot-pack was applied to the surface of the tissue immediately after heating to replicate the prescribed thermal boundary condition used in the FEM simulations. A small temperature drop may be seen when going from the water bath to the tissue surface, but this was minimized by limiting the time of handling and the high heat capacity of the hydrogel pack.

Hence, the heating protocol, the measurement position, the treatment duration and the initial conditions were matched with those defined in the simulations to ensure consistency between the computational and experimental analyses. Finally, the accuracy of the proposed bioheat model was evaluated by comparing the experimental temperature profiles with the simulated transient responses.

3. Results

3.1. Numerical Model Setup and Convergence

A mesh convergence analysis was carried out to find a discretization that gave mesh-independent solutions and kept the computational cost reasonable. The temperature at the muscle boundary, which is the main region of interest, was monitored to assess ten meshes ranging from 5,441 to 42,181 elements. Also, as the mesh density increased, the predicted muscle temperature gradually converged, but the computational cost increased proportionally. The predicted temperature in the meshes considered differed only by 0.06 °C, indicating that further mesh refinement led to insignificant changes in the numerical solution. The solution time was 9 s for the coarsest mesh and 24 s for the finest mesh, respectively. The results are summarized in Table 6.

These results indicate that mesh M6 was selected for all subsequent simulations and has 22,183 elements. The mesh provided a balance between computational efficiency and numerical accuracy. The calculated muscle temperature was within 0.01 °C of the finest mesh, but the computational cost was reduced by ~38%.

3.2. Temperature Distribution and Thermal Response

The spatial temperature distribution in the multilayer tissue model was evaluated using the isothermal contours after 15 min of superficial thermotherapy. The thermal response of the M-size anatomical model exposed to hot-pack temperatures of 38 °C and 44 °C is shown in Figure 4.

At an application temperature of 38 °C (Figure 4a) the temperature elevation was mostly limited to the superficial tissue layers with only limited heat penetration to deeper regions. The result was a marginal temperature increase of the muscle boundary, being about 0.14 °C from the initial condition.

A higher hot-pack temperature of 44 °C (Figure 4b) resulted in a more widespread thermal diffusion into the tissue. The resulting thermal gradients penetrated deeper into the multilayer structure, increasing the temperature of the muscle boundary by ~1 °C after 15 min of treatment. The isothermal contours showed a smoother temperature gradient from the skin surface to the muscle layer than the 38 °C situation.

The comparison between both heating conditions shows that an increase in hot-pack temperature results in increased heat penetration and higher temperatures in the deeper tissue layers considering the evaluated conditions.

3.3. Effect of Adipose Tissue Thickness on Heat Transfer

To evaluate the effect of anatomical variability on heat transfer during superficial thermotherapy, the M-size and XL-size anatomical models were treated under the same conditions (44 °C hot-pack for 15 min) and their thermal responses were compared. The only geometric difference between the models was the adipose tissue layer thickness, which allowed to isolate the influence of this anatomical asymmetry on heat propagation.

The time evolution of the temperature at the muscle boundary is shown in Figure 5. In the M-size model, the temperature at the muscle boundary increased by about 1.03 °C after 15 min of heating, indicating that the thermal stimulus applied reached effectively the deeper layers of tissue. Conversely, the temperature increase of the XL-size model was much lower, just 0.21 °C under the same boundary conditions. In this configuration, temperature gradients were mainly confined to the superficial tissues and did not propagate much towards the muscle layer.

Comparison shows that geometric differences related to adipose tissue thickness significantly alter the predicted thermal response regardless of the heating conditions are identical. The time to reach the same temperature increase in the XL-size model was much longer than in the M-size model. From the results it was extrapolated that the duration of the treatment of the XL-size configuration should be longer than 1 h to increase 1 °C at the muscle boundary under the same thermal conditions.

3.4. Experimental Validation of the Finite Element Model

The M-size model was selected for the experimental validation since this model showed the highest temperature increase at the muscle boundary in the numerical simulations. This set-up provided the most appropriate conditions for comparison between the predicted thermal response and experimental measurements made from the ex vivo porcine tissue.

Figure 6 (a) shows the temporal evolution of the experimentally measured temperatures at various depths in the tissue. A progressive increase in temperature in the muscle layer was observed during the heating period, with the final temperature increase being approximately 1.16 °C after 16 min of thermal stimulation.

Figure 6 (b) compares the experimental measurements with the simulation results directly on the muscle boundary. Both datasets showed a similar temporal evolution with a continuous increase of the temperature during the treatment period. The average difference between the simulated temperature and the experimental temperature was 0.81 °C indicating a proper agreement between the computational predictions and the experimental observations.

The experimental results confirm the ability of the proposed FEM to reproduce the transient thermal response of the multilayer tissue under the evaluated heating conditions. The good agreement between numerical predictions and experimental measurements supports the computational framework used for investigating the effect of anatomical asymmetry on heat transfer in superficial thermotherapy.

4. Discussion

This paper demonstrates that anatomical variability, especially in adipose tissue thickness, has a significant impact on heat transfer in superficial thermotherapy. Both numerical simulations and experimental validation demonstrated a significantly different thermal response for the two anatomical configurations under the same heating conditions. In the M-size model, after 15 min heating with a 44 °C hot-pack, the temperature at the muscle boundary was raised around 1 °C, while in the XL-size model only a temperature of 0.21 °C was achieved. These results suggest that geometrical differences related to the thickness of the adipose tissue significantly impact the heat propagation, indicating that standardized thermotherapy protocols based on the assumption of anatomically symmetric conditions may not result in equivalent thermal doses in patients with different body morphologies.

These findings are consistent with previous studies that have demonstrated less deep-tissue heating in subjects with greater adiposity. Petrofsky et al. [5] reported that the transfer of heat to deeper tissues was slowed by an increased thickness of the subcutaneous fat layer. Kominami et al. [13] reported that traditional hot-pack therapy increased superficial temperature mainly and had little effect on the underlying muscle. However, unlike these studies, which were mainly experimental, the present work combines FEM with experimental validation to quantify the effect of systematic changes in adipose tissue thickness on heat propagation. This combined framework offers a controlled setting to assess anatomical asymmetry, allowing comparisons that would be challenging to evaluate experimentally while preserving physiological relevance through ex vivo validation.

The observed behavior is interpreted from a heat transfer point of view as a consequence of the relatively low thermal conductivity of adipose tissue, which limits the heat propagation into deeper regions and promotes the heat retention in the superficial layers. Thus, the increased thickness of the adipose tissue decreases the effective thermal gradient to the muscle, leading to the attenuated thermal response observed in the XL-size model. These results are consistent with previous work that has found subcutaneous fat to be an important contributor to thermal diffusion during superficial heating [10]. The simulations results also demonstrate the increase of insulating effect with the increase of adipose tissue thickness, emphasizing the significance of incorporating patient-specific anatomical features in computational bioheat models.

These results have one important clinical implication, the duration of treatment. Extrapolation of the XL-size model suggests that under the same thermal conditions, more than one hour of continuous heating would be required to increase the temperature by 1 °C at the muscle boundary. These exposure times are not recommended clinically, since prolonged application of temperatures around 44–45 °C increases the risk of the skin to burns and thermal injury [20]. Thus, extending the treatment time cannot be considered as a feasible way to compensate the lower heat penetration into subjects with thicker adipose layers. These results instead indicate the need to consider patient morphology in the selection of heating conditions as part of treatment protocols to achieve consistent therapeutic outcomes.

Experimental validation further supports the predictive capability of the proposed computational framework. The numerical and experimental temperature profiles exhibited similar temporal behavior with an average difference of about 0.81 °C at the muscle boundary. This difference can be explained by several factors. First, in the numerical model, an idealized constant temperature boundary condition is used, while the passive hydrogel hot-packs cool down gradually during the treatment because heat is transferred to the tissue. This assumption is an upper-bound heating scenario that simplifies the analysis, while permitting comparison between anatomical configurations. More importantly, the ex-vivo experimental model lacks active blood perfusion, a known physiological mechanism that influences heat transport in living tissues. Perfusion was included in the numerical simulations using the Pennes bioheat formulation, but the absence of perfusion in the experimental specimens will inevitably change the thermal dynamics. Finally, the small differences between the anatomical dimensions of the porcine specimens and the computational geometry and the measurement uncertainties may also contribute to the observed deviations. Notwithstanding such limitations, the excellent agreement in temporal trends indicates that the presented model reflects the major governing mechanisms of heat transfer during superficial thermotherapy.

To summarize, this study shows that anatomical asymmetry due to differences in adipose tissue thickness has a strong influence on heat transfer during superficial thermotherapy. The results highlight the limitations of standardized treatment protocols that do not consider patient morphology and the importance of anatomically informed computer models in predicting patient-specific thermal responses. In addition to the validation of the proposed finite element framework, this work provides quantitative evidence that systematic variations in the tissue geometry alone can significantly modify the delivered thermal dose, supporting the development of personalized thermotherapy protocols based on individual anatomical characteristics. A comparison between the present and previous experimental and computational studies is summarized in Table 7.

5. Conclusions

In this study, we conducted a combined computational and experimental study of heat transfer during superficial thermotherapy. A finite element model of the lower limb with anatomically representative multilayer tissues was used. The proposed computational framework was validated by ex vivo porcine experiments and demonstrated comparable temporal temperature trends between numerical predictions and experimental measurements.

The results indicate that the thickness of the adipose tissue significantly impacts the propagation of heat to the deeper tissues. The M-size model under the same heating conditions showed an increase of about 1 °C at the muscle boundary after 15 min heating with a hot-pack of 44 °C, while the XL-size model showed only an increase of 0.21 °C. These results show that anatomical asymmetry due to differences in adipose tissue thickness can significantly impact the delivered thermal dose, suggesting that standardized thermotherapy protocols based on uniform anatomical assumptions may not offer equivalent therapeutic effects across patients.

To summarize, this study provides quantitative support for considering patient morphology in the design of superficial thermotherapy protocols. The proposed FEM framework can assess the heat propagation under different anatomical configurations and can be regarded as a useful tool for the development of anatomically informed and patient-specific treatment strategies. Future work will include the incorporation of temperature-dependent blood perfusion models, transient hot-pack cooling and in vivo validation to further improve the physiological realism and clinical applicability of the proposed methodology.