Staged Effective Medium Modeling and Experimental Validation for Rock Thermal Conductivity

1. Introduction

The thermal conductivity of rocks is defined as the ability of the material to conduct heat, or the rate of heat transfer through a unit area under a unit temperature gradient (W m⁻¹ K⁻¹), which is a key parameter in geoscience and engineering that governs heat transfer. (Carslaw and Jaeger, 1959). Thermal conductivity controls geothermal gradients, energy transport in reservoirs and repositories, and heat dissipation in tunnels, foundations, and buried infrastructure (García-Noval et al., 2024; Preux and Malinouskaya, 2021). Geothermal energy relies on extracting subsurface heat; the thermal conductivity of the host rocks directly determines the efficiency of heat flow toward production wells and thus the sustainability and power potential of geothermal systems (Ehleringer, 1989; Clauser & Huenges, 1995; Raymond et al., 2022; Mokheimer et al., 2019). Typical thermal conductivities of sedimentary formations span approximately 1–6 W m⁻¹ K⁻¹, reflecting substantial lithological variability: from high-conductivity quartz-rich rocks (~7.7 W m⁻¹ K⁻¹) to low-conductivity clay-rich units (~2 W m⁻¹ K⁻¹) (Horai, 1971; Brigaud & Vasseur, 1989). The commonly used thermal conductivity of natural quartz is ~7.7 W m⁻¹ K⁻¹ (Caluser and Huenges 1995; Brigaud & Vasseur, 1989). However, when field and lab data are extrapolated to zero porosity in quartz-rich sandstones, it typically falls below this value due to the influence of lower thermal conductivity minerals, thermal contact resistance at grain boundaries, the assumption of perfect grain packing at zero porosity and the exclusion of the influence of micropores and cements. Such variability introduces significant uncertainty into geothermal and heat-flow modeling, as even modest deviations in thermal conductivity can markedly influence subsurface temperature predictions. Consequently, recent investigations underscore the need to accurately constrain thermal conductivity distributions to improve the reliability of modeled geothermal gradients and subsurface thermal regimes (Hermans et al., 2014; Raymond, 2018).

Laboratory thermal conductivity measurements are broadly classified into steady-state and transient-state methods, depending on whether a constant or time-dependent temperature gradient is imposed across the sample (Carslaw & Jaeger, 1959; Parker et al., 1961). Steady-state techniques, such as the divided-bar and guarded-hot-plate methods, determine thermal conductivity once thermal equilibrium is achieved, yielding high accuracy but requiring longer equilibration times and precise control of boundary conditions (Bullard, 1939; Sass et al., 1971). In contrast, transient methods, including the line-source (needle probe) and laser-flash techniques, infer conductivity from the dynamic temperature response to a heat pulse or step input, allowing for faster measurements and application to a broader range of materials (Jaeger, 1958; Parker et al., 1961). Steady-state methods are limited by sample geometry constraints and susceptibility to heat losses, while transient methods are affected by uncertainties in contact resistance, heat capacity estimation, and assumptions of homogeneity (Popov et al., 1999; Hammerschmidt & Meier, 2006).

In situ determinations can be performed using downhole thermal-conductivity logging tools (Beck et al., 1971; Burkhardt et al., 1995; Hyndman et al., 1979; Kuriyagawa et al., 1983); however, these measurements are typically spatially discontinuous, technically demanding, and economically restrictive. In shallow subsurface investigations, the thermal response test (TRT) has become a widely adopted technique for assessing the effective thermal properties of borehole environments (Gehlin, 2002). For deeper boreholes, optical frequency-domain reflectometry (OFDR) enables direct, high-resolution profiling of thermal conductivity along the borehole wall (Lehr & Sass, 2014). Measurements in unconsolidated sediments require controlled variation of water content and the application of incremental compaction pressures to account for the strong dependence of heat transfer on fluid saturation and packing density (Sass and Stegner, 2012). These methods have led to wide variations in the values of thermal conductivity measurements, which inhibit model calibration and increase prediction errors in modeling subsurface heat flow (Dalla Santa, 2020). Often datasets are obtained using diverse experimental setups and methodological approaches (Dalla Santa, 2020). Also datasets are reported inconsistently, as some studies report only median values, while others provide full ranges or extreme bounds (minimum and maximum) with an uncertainty range of approximately 5–50% (Dalla Santa, 2020). In addition, this dispersion in thermal conductivity values reflects the influence of the microstructural and saturation properties of the rock (Schön, 2011). These differences in data reporting and measurement techniques preclude rigorous statistical analysis.

The bulk thermal conductivity of a rock is governed by its mineralogical composition, porosity, pore fluid characteristics, and the geometric configuration of its pore network (Sun et al., 2017). Sandstones, shales, and carbonates with comparable bulk compositions may nonetheless exhibit distinctly different thermal conductivities due to variations in pore structure and connectivity (Waqas et al., 2022; Hu et al., 2024). In particular, the presence of anisotropic pore systems and microcrack networks can enhance or inhibit heat transfer depending on their orientation relative to the thermal gradient (Popov et al., 1999; Fuchs et al., 2013). Furthermore, fluid saturation exerts a strong influence, as the substitution of air by water or brine substantially increases the effective thermal conductivity through improved heat transfer across pore spaces (Brigaud & Vasseur, 1989; Clauser & Huenges, 1995; Mielke et al., 2017). The application of induced stress closes microcracks and increases grain-to-grain contact areas, enhancing heat transfer efficiency and improving thermal conductivity (Abdulagatova et al., 2009a). Temperature exerts an opposite effect on thermal conductivity: higher temperatures enhance phonon scattering, intensifying lattice vibrations and disrupting the orderly transfer of thermal energy through the crystal lattice, thereby diminishing effective heat transport (Labus et al., 2018; Vosteen & Schellschmidt, 2003). In porous, fluid-saturated rocks, temperature-induced changes in fluid properties (e.g., viscosity and density) may partially counteract this reduction by improving conductive pathways within the pore network (Popov et al., 1999; Clauser & Huenges, 1995). Consequently, accurate geothermal modeling requires applying temperature-dependent corrections to measured thermal conductivity to obtain precise estimates of subsurface heat flow. These interrelated thermal effects highlight the importance of integrating mineralogical, textural, geomechanical, and petrophysical parameters to predict the thermal properties of geological formations in geothermal systems and basin-scale heat-flow models (Vosteen & Schellschmidt, 2003; Fuchs et al., 2015; Labus et al., 2018).

A substantial body of research has employed theoretical and empirical mixing laws to estimate the effective thermal conductivity of rocks from petrophysical parameters. In two-phase systems comprising a solid matrix and pore space, the simplest theoretical limits are provided by the Wiener bounds (Wiener, 1912), which represent end-member limits in which heat flow is either parallel or perpendicular to bedding. For isotropic media, the Hashin–Shtrikman (HS) bounds provide the most restrictive theoretical limits on effective thermal conductivity, offering a rigorous framework for evaluating intermediate behavior between them (Hashin and Shtrikman, 1962). Between the Wiener and Hashin–Shtrikman bounds, various rock physics models have been applied to estimate an “effective” thermal conductivity of rocks. The Maxwell–Eucken models provide two complementary formulations for heat transfer in rocks, distinguishing whether the solid matrix or the pore fluid forms the continuous phase containing isolated spherical inclusions (Maxwell, 1873; Eucken, 1940; Levy, 1981). Both models rely on dilute-sphere assumptions in which inclusions are non-interacting and randomly distributed, making them most valid at low porosity where phase interactions remain negligible. Although simple, they offer useful first-order estimates of effective thermal conductivity and underpin more advanced effective-medium theories and empirical corrections commonly used in geothermal and basin-scale heat-flow modeling (Brigaud and Vasseur, 1989; Fuchs et al., 2013). Bruggeman’s self-consistent formulation is a widely used symmetrical effective medium theory that predicts composite thermal conductivity by treating both the solid and fluid phases as statistically equivalent constituents embedded within an infinite surrounding medium (Bruggeman, 1935; Myers, 1989). This approach incorporates higher inclusion concentrations and phase-interaction effects, yielding more realistic bulk thermal-conductivity estimates in heterogeneous or isotropic materials (Fuchs et al., 2013; Vosteen & Schellschmidt, 2003). Differential Effective Medium theory models composite behavior by incrementally introducing inclusions into a uniform host and updating the effective properties at each step. This enables modeling the thermal conductivity of mixtures more realistically than static mixing laws (Berryman, 1995; Norris, 1985). Although DEM models often yield improved predictions for rocks with moderate porosity or anisotropy, it may diverge from observed thermal conductivities when microstructural heterogeneity or fabric orientation strongly control heat transport (Fuchs et al., 2013; Mavko et al., 2020). The geometric mean model (Woodside and Messmer, 1961) estimates bulk thermal conductivity as the geometric average of the matrix and pore conductivities, weighted by their volume fractions, assuming an isotropic, randomly distributed two-phase system. Although it lacks explicit microstructural parameters, it often matches the performance of more complex EMT formulations for consolidated rocks with moderate porosity. This typically predicts sandstone thermal conductivity within about ten to twenty percent (Brigaud and Vasseur, 1989). Other studies show that weighted geometric means often provide the best empirical fits to laboratory data (Fuchs et al., 2018). More sophisticated models yield only marginal improvements (Preux and Malinouskaya, 2021). Zimmerman (1989) attempted to modify the HS upper bound by incorporating a textural parameter, the pore aspect ratio, which directly influences the distribution of heat flow paths. However, most reservoir simulators still use very simple laws that fail to account for the pore structure in reservoirs (Preux and Malinouskaya, 2021). Pimienta et al. (2014) developed two Effective Medium Theory–based models, Model 2 and Model 3, to predict thermal conductivity in porous and microcracked rocks by explicitly incorporating crack density and aspect ratio into the upscaling scheme. Model 2 embeds compliant cracks into a high-pressure host medium, assuming all cracks are closed, reducing assumptions about mineralogy and pore shape and enabling more reliable inversion of crack parameters as pressure decreases. In contrast, Model 3 assumes a monomineralic matrix with spherical, stress-insensitive pores and uses the Mori–Tanaka homogenization scheme, introducing additional assumptions and generating larger uncertainties in predicted crack geometry and thermal conductivity. The authors also show that although thermal conductivity can be reasonably inferred from elastic wave velocity, mineral composition, and porosity, the predictions remain limited by assuming a simplified pore geometry, idealized microstructural assumptions, and the inherent non-uniqueness of the inverse crack-parameter problem. Thes models offer useful theoretical constraints on thermal conductivity, but none adequately capture the complexity of natural rocks (Hu et al., 2024; Preux & Malinouskaya, 2021). Numerous studies have used indirect methods for estimating thermal conductivity from petrophysical parameters such as bulk density, porosity, and mineralogical composition (Fuchs et al., 2013; Goutorbe et al., 2006; Gross & Combs, 1976; Hartmann et al., 2005). These regression-based empirical formulations are inherently limited, as their validity is generally restricted to the specific lithologies and environmental conditions for which they were calibrated. Comprehensive summaries are provided by Fuchs et al. (2013) and Schön (2015).

Elastic-wave velocity in rocks exhibits similar sensitivity to mineral composition, porosity, fluid saturation, and fabric anisotropy as thermal conduction, and seismic velocity measurements are often employed as indirect proxies for estimating thermal conductivity (El Sayed & El Sayed, 2019). Empirical relationships between thermal conductivity and P-wave velocity (Vp) have been established for specific lithologies, indicating that for dry sandstones and volcanites, thermal conductivity can be estimated from Vp with an accuracy of approximately ±0.5 Wm⁻¹K⁻¹ (Mielke et al., 2017). They also showed that increasing water saturation leads to simultaneous increases in both thermal conductivity and P-wave velocity, with porosity exerting a dominant influence on the strength of this relationship. However, Mielke et al. (2017) reported weak correlations between thermal conductivity and P-wave velocity in clastic sediments, as evidenced by low coefficients of determination (R² = 0.04–0.48) across sandstone types. This indicates that acoustic velocity alone provides only limited predictive capability for thermal conductivity, reflecting the complex influence of porosity, mineralogy, and microstructural heterogeneity. Subsequent investigations have proposed linear and non-linear empirical correlations between these properties. Such models remain largely statistical in nature and lack explicit physical grounding (Gegenhuber & Schoen, 2012). Several studies have explored the interdependence of thermal conductivity and elastic-wave velocity through their mutual sensitivity to porosity, mineralogy, fluid saturation, and confining stress (Esteban et al., 2015; Orlander et al., 2018). These approaches rely primarily on direct regression between velocity and thermal conductivity.

2. Materials and Methods

A custom-designed, calibrated, divided-bar 1-D system for high-precision thermal conductivity measurements of small cylindrical rock discs, equipped with multiple PID-controlled guard heaters was developed. The samples are confined between two heated reference samples in a 1-D axial setup, which minimizes any buoyancy-driven circulation. The thermal conductivity apparatus and pulse-transmission ultrasonic setup were then used to measure thermal conductivity and acoustic velocities under increasing isostatic stress (200-2000psi) for Fontainebleau sandstones, and laboratory-constructed sand packs (of varying grain size and mineralogy) in both dry and saturated states.

2.1. Sample Description and Preparation

Fontainebleau Sandstone

The Rupelian Fontainebleau Sandstone, was deposited in a marine shoreface subaerial-eolian environment, is quartz-rich (quartz >99%), and known for its high degree of cementation, and uniform grain fabric (Thiry et al., 1998). Three cylindrical cores (2.54 cm in diameter and ~6 cm in length) were drilled from large rectangular blocks using a water-cooled diamond-tipped coring bit. The samples were precision-ground and polished at both ends (right cylindrical geometry) to ensure good thermal contact for thermal conductivity measurements and reliable coupling of acoustic transducers to the samples. The core plugs were vacuum-oven-dried at 60 °C. The dry weight, bulk volume, and dry bulk density were measured using calipers and a precision balance. Grain density and porosity were determined using a nitrogen expansion (Boyle’s Law) porosimeter. The samples were first placed in a vacuum chamber for one hour, followed by the injection of degassed water. The sample then remained at 800 psi pressure for 12 hours to ensure full saturation.

The sand-pack specimens were prepared using Brazos River sand, which is quartz-rich and characterized by subangular grains.

Table 1. Sand Pack Properties.

Table 1. Sand Pack Properties.

Three narrow grain-size fractions (45–75 μm, 180–212 μm, and 355–425 μm) were prepared by sieving to ensure consistent particle-size distributions. Samples were sleeved in Viton (2.5 cm inner diameter and ~3 cm length). The initial porosity was almost independent of grain size.

A series of isostatic compaction experiments was performed on unconsolidated Brazos sand specimens with a 355–425 $\mathit{\mu }$m grain-size distribution, prepared as both monomineralic quartz packs and mixed-mineral frameworks. Both thermal conductivity and ultrasonic wave velocities were measured. We calculated the final porosity (${\mathsf{\varphi }}_{\mathbf{f}\mathbf{i}\mathbf{n}\mathbf{a}\mathbf{l}})$ by:

$${\mathsf{\varphi }}_{\mathbf{f}\mathbf{i}\mathbf{n}\mathbf{a}\mathbf{l}}=1-{\displaystyle \frac{{(1-\mathsf{\varphi }}_{\mathbf{i}\mathbf{n}\mathbf{i}}){\mathbf{L}}_{\mathbf{i}\mathbf{n}\mathbf{i}}}{{\mathbf{L}}_{\mathbf{f}\mathbf{i}\mathbf{n}\mathbf{a}\mathbf{l}}}}$$

(1)

where ${\mathit{\varphi }}_{\mathit{i}\mathit{n}\mathit{i}}$ is initial porosity (fraction), ${\mathit{\varphi }}_{\mathit{f}\mathit{i}\mathit{n}\mathit{a}\mathit{l}}$ is porosity after compaction, ${\mathit{L}}_{\mathit{i}\mathit{n}\mathit{i}}$ is initial sample length and ${\mathit{L}}_{\mathit{f}\mathit{i}\mathit{n}\mathit{a}\mathit{l}}$= final sample length (after compaction). The assumptions are that the sample is isotropic and the grain volume is constant.

Thermal Conductivity (λ) Setup, Calibrations, and Measurements

A steady-state divided-bar apparatus, integrated with a servo-controlled pressure vessel, was used to measure thermal conductivity. The apparatus combines two subsystems: (1) a divided bar thermal conductivity assembly with actively tuned guard heaters across a 1-inch diameter sample column and (2) a triaxial stress geomechanical press employing a programmable isostatic pressure (+1% error) to the sample stack (Figure 2). The core sample was placed between two high-thermal conductivity aluminum discs (Figure 2). A constant heat flux was applied to the top disc by an electric heater, while the bottom disc was held at a lower reference temperature (Figure 2). The temperature drop across the sample was measured by precision RTD sensors embedded in the aluminum discs, and a one-dimensional Fourier’s law (steady state) was used to compute the sample’s thermal conductivity (Figure 2). This divided-bar stack consists of (in order): a heat source (temperature-controlled bath or heater), an upper heat flow meter (aluminum disc), the rock sample, a lower heat flow meter, and a heat sink (International, 2013). The stack was surrounded by thermal insulation to minimize convective heat losses. Guard heaters were applied to eliminate radial heat flow. Measurements were made on both dry and brine saturated samples.

The improved divided-bar apparatus uses three independently PID (Proportionality-Integral-Derivative)-controlled guard heaters (two around the reference sections and one on the sample) to minimize radial heat losses. The steady-state condition for the divided bar apparatus for symmetric heat flux is:

$$\mathit{e}\left(\mathit{t}\right)=\nabla {\mathit{G}}_{\mathit{t}\mathit{o}\mathit{p}}-\nabla {\mathit{G}}_{\mathit{b}\mathit{o}\mathit{t}\mathit{t}\mathit{o}\mathit{m}}$$

(2)

where $\mathit{e}\left(\mathit{t}\right)$ is the input error to the PID controllers (Figure 2). The $\nabla {\mathit{G}}_{\mathit{t}\mathit{o}\mathit{p}}$ and the ${\nabla \mathit{G}}_{\mathit{b}\mathit{o}\mathit{t}\mathit{t}\mathit{o}\mathit{m}}$ terms are the measured temperature gradients. The PID control algorithm dynamically adjusts heater power to maintain uniform gradients between the reference sections. Simply, it measures the imbalance in heat flow symmetry (Figure 2). The time-discretized PID for the top and bottom sections is based on the PID control output ($\in \mathit{o}\mathit{r}\mathit{u}\left(\mathit{t}\right))$ and is adjusted to supply heat flux as:

$$\in =\mathit{u}\left(\mathit{t}\right)={\mathit{K}}_{\mathit{p}}\mathit{e}\left[\mathit{K}\right]+{\mathit{K}}_{\mathit{I}}{\mathit{T}}_{\mathit{S}}{\sum }_{\mathit{I}=0}^{\mathit{k}}\mathit{e}\left[\mathit{i}\right]+{\mathit{K}}_{\mathit{d}}{\displaystyle \frac{\mathit{e}\left[\mathit{k}\right]-\mathit{e}\left[\mathit{k}-1\right]}{{\mathit{T}}_{\mathit{s}}}}$$

(3)

where $\mathit{e}\left[\mathit{K}\right]$ is the current error at time step $\mathit{k}$, i.e., the difference between the target gradient (setpoint) and the measured gradient, ${\mathit{K}}_{\mathit{p}}\mathit{e}\left[\mathit{K}\right]$ is the proportional term, which applies a correction proportional to the current errors (in response to deviations), ${\mathit{K}}_{\mathit{I}}{\mathit{T}}_{\mathit{S}}\sum _{\mathit{I}=0}^{\mathit{k}}\mathit{e}\left[\mathit{i}\right]$ is the integral term, which accumulates past errors to eliminate steady-state offsets, ${\mathit{T}}_{\mathit{s}}$ is the sampling time (time interval between control updates), and ${\mathit{K}}_{\mathit{d}}{\displaystyle \frac{\mathit{e}\left[\mathit{k}\right]-\mathit{e}\left[\mathit{k}-1\right]}{{\mathit{T}}_{\mathit{s}}}}$ is the derivative term, which anticipates future error trends by estimating the rate of change of error, helping to dampen oscillations and overshoot. The proportional, integral, and derivative gains $({\mathit{K}}_{\mathit{p}},{\mathit{K}}_{\mathit{I}},{\mathit{K}}_{\mathit{d}})$ tune the responsiveness, stability, and smoothness of temperature control. The middle (sample) section’s heater power is modeled and balanced in the divided-bar thermal conductivity apparatus using Newton’s law of cooling to account for radial heat loss:

$$\in =\mathit{\alpha }\left({\mathit{T}}_{\mathit{c}}-{\mathit{T}}_{\mathit{e}\mathit{n}\mathit{v}}\right)=\mathit{\alpha }∆\mathit{T}$$

(4)

where $\in$ is the heat flux lost from the sample to the environment, not the main conductive heat passing through the stacked setup, $\mathit{\alpha }$ is the heat transfer coefficient, ${\mathit{T}}_{\mathit{c}}$ is the surface temperature of the middle section, ${\mathit{T}}_{\mathit{e}\mathit{n}\mathit{v}}$ is the ambient temperature. The term $\mathit{\alpha }∆\mathit{T}$ quantifies radial heat loss, which is used to maintain a constant axial heat flow throughout the entire stack. The proportionality that expresses how radial heat losses from each section of the divided-bar apparatus relate to their local temperature differences (ΔT) and surface heat transfer coefficients (α), is described as:

$${\mathit{\alpha }}_{\mathit{t}\mathit{o}\mathit{p}}∆{\mathit{T}}_{\mathit{t}\mathit{o}\mathit{p}}:{\mathit{\alpha }}_{\mathit{m}\mathit{i}\mathit{d}}∆{\mathit{T}}_{\mathit{m}\mathit{i}\mathit{d}}:{\mathit{\alpha }}_{\mathit{b}\mathit{o}\mathit{t}\mathit{t}\mathit{o}\mathit{m}}∆{\mathit{T}}_{\mathit{b}\mathit{o}\mathit{t}\mathit{t}\mathit{o}\mathit{m}}={\in }_{\mathit{t}\mathit{o}\mathit{p}}:{\in }_{\mathit{m}\mathit{i}\mathit{d}}:{\in }_{\mathit{b}\mathit{o}\mathit{t}\mathit{t}\mathit{o}\mathit{m}}.$$

(5)

In other words, sections that are hotter (larger ΔT) or less insulated (larger α) will lose more heat radially. Assuming similar heat loss characteristics for each section Equation 5 is simplified by canceling the terms: ${\mathit{\alpha }}_{\mathit{t}\mathit{o}\mathit{p}},{\mathit{\alpha }}_{\mathit{m}\mathit{i}\mathit{d}}\mathit{a}\mathit{n}\mathit{d}{\mathit{\alpha }}_{\mathit{b}\mathit{o}\mathit{t}\mathit{t}\mathit{o}\mathit{m}}.$ The power balance equation that expresses the heat flux equilibrium in the mid-section power (${\mathit{P}}_{\mathrm{mid}}$) scaled by their temperature ratios gives:

$${\mathit{P}}_{\mathit{m}\mathit{i}\mathit{d}}={\in }_{\mathit{m}\mathit{i}\mathit{d}}=\left({\in }_{\mathit{t}\mathit{o}\mathit{p}}{\displaystyle \frac{∆{\mathit{T}}_{\mathit{m}\mathit{i}\mathit{d}}}{∆{\mathit{T}}_{\mathit{t}\mathit{o}\mathit{p}}}}+{\in }_{\mathit{b}\mathit{o}\mathit{t}\mathit{t}\mathit{o}\mathit{m}}{\displaystyle \frac{∆{\mathit{T}}_{\mathit{m}\mathit{i}\mathit{d}}}{∆{\mathit{T}}_{\mathit{b}\mathit{o}\mathit{t}}}}\right)/2.$$

(6)

This corrects for heat losses in the mid-section by maintaining constant axial heat flux through the top and bottom standards.

The measurement stack is illustrated in Figure 1.

Figure 1. Cartoon of the experimental sample and reference stack. The Aluminum disks allow accurate temperature measurements to be made.

Figure 1. Cartoon of the experimental sample and reference stack. The Aluminum disks allow accurate temperature measurements to be made.

Figure 2. Temperature distribution along the z axis of the divided bar apparatus. The PID loops ensure that the gradients across the top and bottom standards are identical to ensure axial heat flow.

Figure 2. Temperature distribution along the z axis of the divided bar apparatus. The PID loops ensure that the gradients across the top and bottom standards are identical to ensure axial heat flow.

Figure 3. Viton sleeved high-pressure thermal conductivity measurement stacked setup with embedded reference discs, thermocouples, and heater wiring (left) and thermal insulated stacked setup (right).

Figure 3. Viton sleeved high-pressure thermal conductivity measurement stacked setup with embedded reference discs, thermocouples, and heater wiring (left) and thermal insulated stacked setup (right).

A fused quartz standard sample was used during PID-controlled heater calibration. Figure 2. shows the time dependent evolution of the temperature for the four thermocouples (T₁–T₄) monitored over a 6000-second interval, where T₁ (blue) corresponded to the top aluminum section maintained at approximately 65 °C, T₂ (solid black) to the upper reference standard near 48 °C, T3 (black dashed) to the sample region around 41 °C, and T₄ (red) to the bottom aluminum section stabilized near 25 °C. All temperatures reach a steady state between 4000 and 5000 seconds after an initial transient, indicating excellent thermal stability (Figure 4). The small fluctuations (< 0.01 °C) suggest that the PID controllers are maintaining a nearly constant temperature with minimal oscillation, which is a key condition for accurate steady-state thermal conductivity measurement. The zoomed-in view of T1 between 4000 and 5000 seconds shows variations within ±0.0025 °C, or a total range of ΔT1 $\le$ 0.005 °C. These small spikes correspond to micro-adjustments by the PID heater to maintain the specified boundary conditions. The nearly flat profile confirms that thermal equilibrium has been reached, and sensor noise dominates over real thermal drift. The temperature gradient across each section, is held nearly constant, with the upper gradient (∇T1) staying around 875 K/m and the lower gradient (∇T2) showing a slight downward drift but remaining within 845–880 K/m, which represents a change of 0.243 °C between the top and bottom reference sections, oor approximately a 2% error (Figure 4). The measured error of 2.0% indicates high gradient accuracy, reflecting only a minor mismatch between the top and bottom temperature gradients and excellent thermal symmetry. These results confirm that the PID-controlled divided-bar system achieved thermal equilibrium, stable gradients, and minimal thermal gradient asymmetry.

Prior to conducting sample measurements, the apparatus was systematically calibrated using reference materials of known thermal conductivity, namely fused silica (quartz) and Teflon, to evaluate stress- and temperature-sensitivity. The sample–meter-disc interfaces were subjected to a controlled axial preload of approximately 200 psi isostatic stress to establish consistent and reproducible initial thermal contact resistance. Figure 5 illustrates the relationship between heater voltage (control power) and top-surface temperature for both the top and bottom guard heaters during guard heater calibration, comparing fused quartz and Teflon (reference standards). Both materials exhibit a linear increase in voltage with temperature, confirming predictable, stable PID control behavior. The slopes (~0.017–0.020 V/°C) indicate that fused quartz requires slightly higher guard power than Teflon to maintain equivalent gradients, consistent with its higher thermal conductivity (~1.4 W/m·K) compared to Teflon (~0.25 W/m·K). The R² > 0.97 for all regressions indicates excellent linear correlation, validating the proportional control tuning between heating zones. The minor offset between the top and bottom heater voltages confirms balanced but distinct zone responses, as required for thermal symmetry.

The stress dependent measurements were performed at controlled confining pressures (0-1500 psi) using argon as the confining gas. During each test, a small differential (deviatoric) axial stress of ~50 psi was applied to ensure mechanical stability and the accuracy of the axial LVDTs. Data acquisition was conducted at multiple pressure increments after confirming steady-state heat flow. Figure 5 shows the small gradual increase in thermal conductivity with temperature, indicating excellent calibration stability. It also shows a moderate increase in thermal conductivity with pressures up to 1500 psi, consistent with reduced interface resistance at high pressure. These trends confirm that fused quartz is an excellent calibration reference, owing to its low temperature-pressure sensitivity and reproducible response. Figure 6 shows results for Teflon, indicating a minor increase in thermal conductivity with temperature (from ~0.247 W/m·K to ~0.275 W/m·K), though with slightly larger deviations, reflecting greater sensitivity to temperature fluctuations and a lower thermal diffusivity than quartz. It also shows that Teflon’s measured thermal conductivity increases only slightly with increasing confining pressure, likely due to reduced interfacial resistance between elements of the stack. The combined results verify that the divided-bar system maintains a consistent heat-flux balance and predictable thermal performance.

Figure 5. Guard Heater Voltage versus top Temperature (3-Zone Control) for fused quartz (left) and Teflon(right).

Figure 5. Guard Heater Voltage versus top Temperature (3-Zone Control) for fused quartz (left) and Teflon(right).

Figure 5. Fused Quartz calibration with respect to temperature (left) and confining pressure (right).

Figure 5. Fused Quartz calibration with respect to temperature (left) and confining pressure (right).

Each specimen comprising Fontainebleau sandstone and synthetic sandpacks was precisely machined to a length of 1.27 cm and a diameter of 2.5 cm, then mounted in the divided-bar assembly to obtain stress-dependent thermal-conductivity measurements under controlled isostatic loading conditions.

Figure 6. Results of the Teflon calibration with respect to temperature (left) and confining pressure (right).

Figure 6. Results of the Teflon calibration with respect to temperature (left) and confining pressure (right).

Ultrasonic P-wave (compressional) and S-wave (shear) velocities were measured using a stress- and temperature-controlled isostatic cell (Figure 8). Each sample (2.54 in diameter and 5.08 cm in length) was enclosed within a heat-shrink Viton jacket containing O-ring seals and a lateral pore-line tap, and positioned between top and bottom broadband transducers (T₁ and T₂). The jacketed assembly was secured to the sample holder, with the pore line routed to the cell head for fluid control (Figure 7). The top transducer serves as the source, and the bottom transducer as the receiver, both connected to the external control system using 4-pin Lemo feedthroughs. The confining cell was filled with mineral oil and pressurized using an ISCO pump to apply controlled isostatic stress. The Ultrasonic measurements were acquired using the pulse-transmission method, in which a one MHz piezoelectric P/S-wave transducer generates an acoustic pulse that propagates through the sample. The transmitted signal is recorded at the opposite end and digitized using a Tektronix TBS-2000 series oscilloscope. Signal selection between P- and S-wave modes was controlled electronically via the acquisition switch panel. Velocities were measured at multiple confining pressures under both vacuum-dry and water-saturated conditions.

3. Results

3.1. Thermal Conductivity and Acoustic Data for Fontainebleau Sandstone

The dry thermal conductivity (${\mathsf{\lambda }}_{\mathrm{d}\mathrm{r}\mathrm{y}}$) of Fontainebleau sandstone exhibited a nonlinear increase with confining pressure, with the lowest-porosity sample showing the highest thermal conductivity (Figure 8.). The initial lack of stress dependence is attributed to the low confining pressure of the measurements. The water-saturated measurements exhibited a similar nonlinear dependence (Figure 8. right). Figure 9 shows that both compressional (P-wave) and shear (S-wave) velocities exhibit a similar increase in velocity with confining pressure. The cross plots shown in Figure 10 demonstrate the similarity of these effects for the high porosity samples.

Figure 8. Measurements of dry thermal conductivity (λ_dry) and water-saturated thermal conductivity (λ_sat) with confining pressure for Fontainebleau sandstone.

Figure 8. Measurements of dry thermal conductivity (λ_dry) and water-saturated thermal conductivity (λ_sat) with confining pressure for Fontainebleau sandstone.

Figure 9. Confining pressure dependence of (a) P-wave velocity (Vₚ) and (b) S-wave velocity (Vₛ) for Fontainebleau sandstone and 10.64% porosity under dry (red circles) and water-saturated (blue squares) conditions.

Figure 9. Confining pressure dependence of (a) P-wave velocity (Vₚ) and (b) S-wave velocity (Vₛ) for Fontainebleau sandstone and 10.64% porosity under dry (red circles) and water-saturated (blue squares) conditions.

Figure 10. Measured thermal conductivity of dry and saturated sand packs of different grain sizes as a function of confining pressure. Larger grain sizes have lower thermal conductivities.

Figure 10. Measured thermal conductivity of dry and saturated sand packs of different grain sizes as a function of confining pressure. Larger grain sizes have lower thermal conductivities.

3.2. Thermal Conductivity and Acoustic Data for the Sand Packs

In the dry condition (Figure 10) the thermal conductivity of the sand packs increases with decreasing grain size and confining pressure. Water saturation increases the thermal conductivity to nearly twice the dry value. In Figure 11 the shear wave velocities for the sand packs increase linearly with confining pressure. In Figure 12 the three sand packs also exhibit an approximately linear increase in P-wave velocity with increasing confining pressure. The compressional velocity decreased with increasing grain size. Under saturated conditions, P-wave velocities increase, elevating velocities across all grain sizes compared to the dry state.

The thermal conductivity and acoustic velocity measurements cross plot (Figure 13) for the sand packs reveal a grain-size–dependent response. Water saturation increases thermal conductivity to nearly twice the dry value. These behaviors are consistent with the measured pore-volume reductions during isostatic compression as shown in Table 3.

Overall, the behavior across grain sizes reflects a balance between grain contact density (favoring fine grains for λ and V_p) and grain contact stiffness (favoring coarse grains for V_s). All the measurements correlate with the degree of compaction recorded in the pore-volume change measurements recorded in Table 3. The extremely high correlation coefficients demonstrate the close connection between these measurements.

4. Discussion

4.1. Formulation of the Staged Differential Effective Medium (SDEM) Model

The Staged Differential Effective Medium (SDEM) framework was introduced by Myers (1989) for conductive inclusions and later extended by Myers and Hathon (2009), to model elastic moduli. The method introduces an interpolation parameter $L$, where $L=0$ corresponds to the Reuss bound and $L=1$ describes the Voight bound. The staged approximation, allows the model to track rock properties starting from a critical porosity to the zero-porosity endpoint (Nur et al., 1998). Porosity is modified in sequential integration steps. The first stage uses $L=0$ (suspension) up to the critical porosity, after which a value of $L>0$ is applied to reflect the properties of the framework grains as they become load bearing. This SDEM model incorporates critical porosity behavior, and maintains consistency with Gassmann fluid substitution. L corresponds to measurable microstructural attributes, such as contact length and grain-fabric evolution (Myers and Hathon 2009). This staged formulation allows the predicted properties to evolve with changing porosity, rather than imposing a single mixing relation across the full porosity range (Simone et al., 2020; Villarroel et al., 2022). Starting from a percolating phase the ordering of the integration steps is determined by the length scales present in the model. The ordering goes from small length scales to large.

Thermal conductivity, like bulk modulus, depends on grain-to-grain contact geometry and is similarly bounded by the Reuss (series or iso-heat flux) and Voigt (parallel or iso-temperature gradient) limits. Substituting thermal conductivity for bulk modulus in the SDEM model is justified because both are derived from an application of mixture theory. The thermal conductivity, $\mathsf{\lambda }$, is then given by:

$$\mathsf{\lambda }={\mathsf{\lambda }}_{\mathrm{h}}{\displaystyle \frac{1+\left(\frac{{\mathsf{\varphi }}_0-\mathsf{\varphi }}{{\mathsf{\varphi }}_0}\right)·\left(\frac{{\mathsf{\lambda }}_{\mathrm{i}}-{\mathsf{\lambda }}_{\mathrm{h}}}{{\mathsf{\lambda }}_{\mathrm{h}}}\right)·\mathrm{L}}{1-\left(\frac{{\mathsf{\varphi }}_0-\mathsf{\varphi }}{{\mathsf{\varphi }}_0}\right)·\left(\frac{{\mathsf{\lambda }}_{\mathrm{i}}-{\mathsf{\lambda }}_{\mathrm{h}}}{{\mathsf{\lambda }}_{\mathrm{i}}}\right)·\left(1-\mathrm{L}\right)}}$$

(7)

where ${\varphi }_0$ is initial porosity, $\mathsf{\varphi }$ is the current porosity, ${\mathsf{\lambda }}_{\mathrm{i}}$ is the thermal conductivity of the inclusion (W/m·K), ${\mathsf{\lambda }}_{\mathrm{h}}$ is the host thermal conductivity (W/m·K), and L is the microstructural fitting parameter. The two-stage model for the effective saturated thermal conductivity (${\mathsf{\lambda }}_{\mathbf{s}\mathbf{a}\mathbf{t}})$ in the Reuss average form is:

$${\displaystyle \frac{{\mathsf{\lambda }}_{\mathrm{s}\mathrm{a}\mathrm{t}}}{{\mathsf{\lambda }}_{\mathrm{i}}-{\mathsf{\lambda }}_{\mathrm{s}\mathrm{a}\mathrm{t}}}}={\displaystyle \frac{{\mathsf{\lambda }}_{\mathrm{d}\mathrm{r}\mathrm{y}}}{{\mathsf{\lambda }}_{\mathrm{i}}-{\mathsf{\lambda }}_{\mathrm{d}\mathrm{r}\mathrm{y}}}}+{\displaystyle \frac{{\mathsf{\varphi }}_{\mathrm{c}}}{\mathsf{\varphi }}}·{\displaystyle \frac{{\mathsf{\lambda }}_{\mathrm{L}}}{{\mathsf{\lambda }}_{\mathrm{i}}-{\mathsf{\lambda }}_{\mathrm{L}}}}$$

(8)

$${\mathsf{\lambda }}_{\mathrm{L}}={\left({\displaystyle \frac{{\mathsf{\varphi }}_{\mathrm{c}}}{{\mathsf{\lambda }}_{\mathrm{c}}}}+{\displaystyle \frac{1-{\mathsf{\varphi }}_{\mathrm{c}}}{{\mathsf{\lambda }}_{\mathrm{i}}}}\right)}^{-1}$$

(9)

Equations 8. And 9. Give a convenient form for calculating the influence of varying the properties of the saturated phases.

4.2. Modeling Fontainebleau Sandstone

This SDEM model enables the calculation of thermal conductivity of saturated samples (${\mathsf{\lambda }}_{\mathrm{s}\mathrm{a}\mathrm{t}}$) from dry thermal conductivity measurements, analogous to a Gassmann’s fluid substitution.

Figure 14. shows the thermal conductivity data measured at 1250 psi. The SDEM model uses a quartz matrix thermal conductivity of 7.69 W/m·K for the host phase (${\mathsf{\lambda }}_{\mathrm{h}}$), and a critical porosity of 40%. For both the dry and saturated thermal conductivity data, the best fit is to a single L value of 0.32.

As shown in Figure 15, plotting the thermal conductivity and the bulk modulus on a normalized plot allows a direct comparison between the two measurements. The thermal and elastic properties exhibit similar trends governed by a similar model parameter. The experimental points for the two measurements cluster between the SDEM model curves for L = 0.3–0.5, demonstrating that both heat flow and elastic deformation respond to similar but not identical physical properties. This cross-property calibration capability, demonstrates that the model can be used to jointly interpret both types of measurements. This underscores the SDEM model’s value as a unified framework without assuming pore geometries such as ellipsoids or cracks as assumed in prior effective medium theories (Pimienta et al., 2014; Mavko et al., 2020).

Figure 14. Measured thermal conductivity of dry and saturated Fontainebleau sandstone plotted against porosity along with SDEM model predictions. Solid curves represent SDEM models for various geometric parameters L and critical porosity ϕ_c = 40%. The squares are the predicted saturated values calculated by fluid substitution from the dry measurements. The best fit for all the measurements gives L=.32.

Figure 14. Measured thermal conductivity of dry and saturated Fontainebleau sandstone plotted against porosity along with SDEM model predictions. Solid curves represent SDEM models for various geometric parameters L and critical porosity ϕ_c = 40%. The squares are the predicted saturated values calculated by fluid substitution from the dry measurements. The best fit for all the measurements gives L=.32.

Figure 15. Normalized properties for dry thermal conductivity (λ) and bulk moduli (K) in SDEM model for Fountainbleu sandstone.

Figure 15. Normalized properties for dry thermal conductivity (λ) and bulk moduli (K) in SDEM model for Fountainbleu sandstone.

4.3. SDEM Modeling of Sand Packs

Figure 16. Modeled L-curves for compaction of coarse quartz sand for dry and saturated conditions. Fluid substitution model (square symbols) is included to the right.

Figure 16. Modeled L-curves for compaction of coarse quartz sand for dry and saturated conditions. Fluid substitution model (square symbols) is included to the right.

Figure 17. SDEM-fitted L-curves for progressive compaction of coarse quartz sand. The best fit is for L=.12.

Figure 17. SDEM-fitted L-curves for progressive compaction of coarse quartz sand. The best fit is for L=.12.

To examine the impact of porosity changes under compaction on the thermal conductivity and bulk modulus samples were compressed isostatically. As shown in Figure 17. the large grain size sand packs were compressed isostatically in steps of 3500psi up to 10,500psi.

The SDEM modeling results demonstrate that the compaction of the coarse quartz sand (355–425 µm) increases both thermal conductivity (λ) and bulk modulus (K). The model captures the porosity-dependent trends (Figure 13 and Figure 14, and Table 4). A single parameter, L ~ 0.12, consistently fits both λ and K for both the dry and saturated samples of the sand-packs. As porosity decreases, λ and K increase along a constant L value. At the highest compaction level (Stage 4), the model shows a modest change in the associated L value, a discrepancy likely arising from the significant grain crushing that occurs at this high stress. The experimentally measured ${\lambda }_{sat}$values for the coarse sand closely match the SDEM model predictions, confirming that saturation effects are well described by fluid substitution rather than requiring a change in the geometric parameter L. These results show that the SDEM framework tracks the evolution of both the thermal conductivity and bulk modulus during compaction. This offers a unified tool for interpreting both the thermal and elastic responses.

4.3.1. Comparison of SDEM Models with Other Models

In this section the result of comparing the following models is compared.

Wiener bounds:

$$\mathsf{\lambda }={\mathsf{\varphi }}_{\mathbf{i}}{\mathsf{\lambda }}_{\mathbf{i}}+{\mathsf{\varphi }}_{\mathbf{h}}{\mathsf{\lambda }}_{\mathbf{h}}$$

(10)

$$\mathsf{\lambda }={\left({\displaystyle \frac{{\mathsf{\varphi }}_{\mathrm{i}}}{{\mathsf{\lambda }}_{\mathrm{i}}}}+{\displaystyle \frac{{\mathsf{\varphi }}_{\mathrm{h}}}{{\mathsf{\lambda }}_{\mathrm{h}}}}\right)}^{-1}$$

(11)

Hashin–Shtrikman upper and lower bounds:

$$\mathsf{\lambda }=\mathsf{\lambda }\mathbf{ᵢ}+{\displaystyle \frac{1-{\mathsf{\varphi }}_{\mathbf{i}}}{\frac{1}{{\mathsf{\lambda }}_{\mathit{h}}-{\mathsf{\lambda }}_{\mathbf{i}}}+\frac{{\mathsf{\varphi }}_{\mathbf{i}}}{3{\mathsf{\lambda }}_{\mathbf{i}}}}}$$

(12)

$$\mathsf{\lambda }={\mathsf{\lambda }}_h+{\displaystyle \frac{{\mathsf{\varphi }}_{\mathrm{i}}}{\frac{1}{{\mathsf{\lambda }}_{\mathrm{i}}-{\mathsf{\lambda }}_h}+\frac{{\mathsf{\varphi }}_h}{3{\mathsf{\lambda }}_h}}}$$

(E13)

Maxwell–Eucken (Clauser and Huenges, 1995):

$$\mathsf{\lambda }={\mathsf{\lambda }}_{\mathit{h}}{\displaystyle \frac{{\mathsf{\lambda }}_{\mathbf{i}}+{2\mathsf{\lambda }}_{\mathit{h}}+2{\mathsf{\varphi }}_{\mathbf{i}}\left({\mathsf{\lambda }}_{\mathbf{i}}-{\mathsf{\lambda }}_{\mathit{h}}\right)}{{\mathsf{\lambda }}_{\mathbf{i}}+{2\mathsf{\lambda }}_{\mathit{h}}-{\mathsf{\varphi }}_{\mathbf{i}}\left({\mathsf{\lambda }}_{\mathbf{i}}-{\mathsf{\lambda }}_{\mathit{h}}\right)}}$$

(14)

Self-consistent / Bruggeman (EMT) (Berryman, 1995; Zimmerman, 1990):

$$0={\mathsf{\varphi }}_{\mathit{h}}{\displaystyle \frac{{\mathsf{\lambda }}_{\mathit{h}}-\mathsf{\lambda }}{{\mathsf{\lambda }}_{\mathit{h}}+\mathsf{\lambda }}}+{\mathsf{\varphi }}_{\mathbf{i}}{\displaystyle \frac{{\mathsf{\lambda }}_{\mathbf{i}}-\mathsf{\lambda }}{{\mathsf{\lambda }}_{\mathbf{i}}+2\mathsf{\lambda }}}$$

(15)

Geometric mean:

$$\mathsf{\lambda }={{\mathsf{\lambda }}_{\mathbf{i}}}^{{\mathsf{\varphi }}_{\mathbf{i}}}{{\mathsf{\lambda }}_{\mathit{h}}}^{{\mathsf{\varphi }}_{\mathit{h}}}$$

(16)

Differential Effective Medium:

$${\displaystyle \frac{\mathbf{d}\mathsf{\lambda }}{\mathbf{d}\mathsf{\varphi }}}={\displaystyle \frac{{\mathsf{\lambda }}_{\mathbf{i}}-\mathsf{\lambda }}{1+\frac{{\mathsf{\lambda }}_{\mathbf{i}}-\mathsf{\lambda }}{3\mathsf{\lambda }}}}$$

(17)

Where: λ_i is the inclusion thermal conductivity, λ_h is the host thermal conductivity, and ϕ_i and ϕ_h is the porosity of the inclusion and host respectivly.

Figure 18 illustrates how the different models (SDEM, Weiner Average, HS average, GeoMean, DEM, EMT, and Maxwell Euken models) compare against measured dry and saturated thermal conductivity. The Voigt-type formulations (SDEM, EMT, and Maxwell Euken models) systematically overpredict λ, while the Reuss-type forms (Weiner Average, HS average) underpredict across all porosities. The geometric-mean curve using the fitted SDEM (L = 0.32) provides the closest match to the dry experimental data. In the saturated case, the models follow similar predictions. The Maxwell–Eucken and EMT produce the highest ${\lambda }_{sat}$, while the Reuss-type and HS averages still give the lowest. All measurements lie comfortably within these two bounds, but all the test models consistently overestimated ${\lambda }_{sat}$. The SDEM with the same L = 0.32 yields predictions that again cluster closely around the measured data.

Figure 19. shows a comparison between the various models and the measured dry thermal conductivity. The expected results are obtained, parallel-type bounds (e.g., Maxwell–Eucken and EMT) produce the highest λ values and consistently overpredict the measurements. Intermediate models including Wiener averages, HS averages, geometric mean, and DEM, also span the central region but still exhibit recognizable biases to too high a prediction. In contrast, SDEM with L = 0.12 yields a prediction that closely matches the measured trend. The plot also shows a comparison between the measured and calculated values for the saturated thermal-conductivity for all of the models. All the models again predict too high a thermal conductivity lying well above the y = x reference line. The geometric-mean model again occupies a mid-range position, capturing the trend but overpredicting the value. The direct SDEM predictions (L = 0.12) track the measurements consistently. Their points all cluster near the reference line, indicating that a single L parameter can accurately describe both dry and saturated measurements.

Figure 19. Comparison between measured and calculated dry (λ_dry,) and water-saturated (λ_sat,) thermal conductivity of the sand-packs using different predictive models. The diagonal dashed line represents the 1:1 reference line.

Figure 19. Comparison between measured and calculated dry (λ_dry,) and water-saturated (λ_sat,) thermal conductivity of the sand-packs using different predictive models. The diagonal dashed line represents the 1:1 reference line.

5. Conclusions

The paper introduces a novel high-pressure divided-bar apparatus, designed to measure thermal conductivity under controlled confining stresses up to 2000 psi while simultaneously acquiring ultrasonic P- and S-wave velocities. This dual-measurement capability provides insight into the coupled evolution of thermal and elastic properties during compaction. New experimental datasets on Fontainebleau sandstone and Brazos sand-packs reveal how compaction driven increases in grain coordination and contact area enhance both thermal conductivity and bulk modulus. These results constitute one of the first systematic demonstrations of how grain size, mineralogy, and progressive compaction effect thermal conductivity in unconsolidated sediments. By correlating thermal conductivity with acoustic velocities measured at identical stress states, the study shows that both properties reflect a similar evolving rock frame. This insight enables the use of sonic-log data to infer thermal conductivity profiles, once the relationship is established.

An important contribution is the development and validation of a Staged Differential Effective Medium (SDEM) model that accurately predicts effective thermal conductivity (ETC) from the dry to fully saturated limit. Across diverse lithologies, the SDEM model reproduces observed thermal conductivity trends with a single parameter, L, that remains valid under both dry and saturated conditions. This enables a Gassmann-style fluid-substitution capability for thermal conductivity, establishing SDEM as a powerful and generalizable tool for modeling thermal conductivity. The classical mixing laws and conventional effective-medium formulations show systematic biases and fail to capture detailed mineralogical and fabric effects.

In summary this work contributes:

• An advanced experimental platform for high-accuracy thermal measurements.
• New data linking compaction mechanics to thermal and elastic transport.
• A validated, physically informed predictive mixing law.

These findings strengthen the theoretical foundation of rock thermal transport, improve prediction accuracy beyond classical approaches, and open new pathways for integrating thermal, mechanical, and acoustic velocity observations in subsurface characterization.