Heat Transfer in Unfrozen and Frozen Porous Media: Experimental Measurement and Pore‐Scale Modeling

In this work, we conducted a mechanistic study on pore‐scale mechanisms controlling heat transfer in partially saturated porous media at unfrozen and frozen conditions. Experimental measurement of effective thermal conductivity (ETC) of simulated partially saturated sediments was carried out to explore the effect of different parameters including pore and overburden pressures, temperature, and water/ice saturation on the pore‐scale mechanisms governing heat transfer in multiphase porous media. The experimental measurements show that the heat transfer is a complex phenomenon affected by several important pore‐scale mechanisms such as particle‐particle conduction and particle‐fluid‐particle conduction, which are governed by water content and distribution, packing structure, wettability characteristics of grains, coordination number, and physical contact among sediment particle. A numerical model was also developed for prediction of ETC using free‐energy lattice Boltzmann model and a space renormalization method. The model predictions were in good agreement with the experimental data, showing that the model is able to reliably estimate ETC with average relative deviations of less than 10%, as it appropriately incorporates the pore‐scale mechanisms influencing ETC. The numerical model predictions were also compared with those of six predictive models available in the literature, and root‐mean‐square errors were calculated to assess its accuracy against the existing models.


Introduction
Effective thermal conductivity (ETC) is a key physical property of heat transport in porous media. ETC studies of dry and fluid-saturated porous media are of utmost importance in a variety of geotechnical, geophysical, and geoenvironmental applications, such as underground transmission lines, geothermal energy resources, solar thermal storage facilities, radioactive waste disposal, geological CO 2 sequestration, and methane recovery from gas hydrate-bearing sediments (Abdulagatova et al., 2009;Li et al., 2013;van Antwerpen et al., 2010). In environmental science context, accurate quantification of ETC for soil and shallow unfrozen/frozen marine sediments could deliver essential information regarding the seafloor stability assessment, sedimentation, submarine slide formation, and global climate change. In particular, ETC plays a key role in geothermal energy engineering (e.g., ground source heat pumps, borehole thermal energy storage, and geothermal heated bridge/pavement) as it affects the kinetic behavior of working fluid temperature, which is essential for efficient design of the processes in this chain. A faster charging and discharging rate could be achieved by selecting a geological medium with higher ETC, while geological media with lower ETC values could be chosen for longer maintenance of the working fluid temperature. Therefore, reliable ETC values for soil are of particular importance for geothermal applications (Zhang & Wang, 2017). ETC is even more crucial in the study of the impact of environmental changes on stability of the natural gas hydrate-bearing sediments and evaluation of the the exploitation associated risks (Waite et al., 2009). For instance, dissociation of stable hydrate-bearing layers could be derived by withdrawal of deep, warm fluids via a conductor that travels across these layers, weakening seafloor sediment stability (Freij-Ayoub et al., 2007;Hovland & Gudmestad, 2001). Here, higher ETC could increase the risk. In contrast, ETC is required to be high enough for choosing materials used in nuclear waste disposal as it results in more efficient cooling and uniform temperature distribution and reduces thermally induced stress which, in turn, improves fatigue properties (Pinheiro et al., 2011). Returning to geological importance of ETC for soils, integrity of caprocks for stored CO 2 could possibly be weakened by thermal stress caused by injected fluid that is also highly dependent on ETC of caprock (Gor et al., 2013).
Thermal properties of porous media depend upon not only the intrinsic thermal properties and volume fraction of each constituent but also pressure, temperature, porosity, and packing structure . Focused on sediments, the other influencing factors include mineralogy, packing structure and gradation, and interparticle physical contact among grains (Dong et al., 2015;Ghanbarian & Daigle, 2016;Sadeghi et al., 2018;Zhang & Wang, 2017). At temperatures below the freezing point, transformation of water to ice (with lower density and higher intrinsic thermal conductivity) together with the presence of unfrozen water may even add more complexities in the system. Experimental studies have been conducted so far to determine ETC of fully and partially water saturated soil samples at unfrozen conditions. The results indicate that ETC of a given soil is largely influenced by quartz and water contents since quartz has the highest thermal conductivity among all the soil minerals and thermal conductivity of air is markedly lower than that of water (Barry-Macaulay et al., 2013;Côté & Konrad, 2005;Lu et al., 2019;Nikolaev et al., 2013;Tarnawski et al., 2009;Venuleo et al., 2016;Zhang et al., 2015;. There are also a number of work investigating the effect of the soil mineral components, particle size, shape, and gradation on its ETC (Campbell et al., 1994;Chen, 2008;Côté et al., 2011). Complexities associated with the presence of clay or other binders were also explored in laboratory (Abuel-Naga et al., 2008de Zárate et al., 2010;Yu et al., 2016). Clay appears to significantly improve the heat conduction by cementing the solid grains and reducing the thermal contact resistance (TCR). It was also found that ETC slightly depends on the temperature (Campbell et al., 1994;Hiraiwa & Kasubuchi, 2000;Liu et al., 2011;Lu et al., 2007). Similarly, studies were also conducted to determine the soil thermal properties at frozen conditions, mainly aiming at understanding the effect of climate change, seasonal freezing-thawing processes, and human activities on thermal response in permafrost regions (Domine et al., 2016;Du et al., 2020;He et al., 2018;Kojima et al., 2018;Overduin et al., 2006;Putkonen, 2003;Wang et al., 2017;. Despite considerable work conducted so far, it is still imperative to conduct experimental studies, particularly at pore scale, to shed light on the mechanisms controlling the heat transfer in unsaturated porous media and more importantly understand their interdependencies. Numerous predictive models have been proposed to estimate ETC of porous materials; however, a unified model or prediction procedure with universal applicability is yet to be developed (Dong et al., 2015). Depending on their principles and assumptions, available predictive models could be categorized into (i) mixing models, which determine ETC by combining the intrinsic thermal conductivity values of the coexisting constituents, typically as a function of volume fraction (Abdulagatova et al., 2009); (ii) empirical models, which build a relationship between ETC and volumetric water content, porosity, and type of the sediment grains by normalizing ETC over the difference between fully saturated and dry state ETC values; (iii) mathematical models, which are adopted from predictive models used to determine other physical properties of composite materials, such as electrical conductivity and hydraulic conductivity; (iv) volume fraction models, which estimate ETC with respect to level of their solid volume fraction (low, medium, and high volume fraction materials); (v) packing structure models, which are developed based on the different packing structures of spheres such as simple cubic, body centered cubic, and face centered cubic; and (vi) pressure-dependent models, which treat the sediment as a standard unit cell containing two contacting particles and take the effect of the compressive pressure into account in order to capture the effect of TCR between grains as the contact area and TCR was shown to markedly control the heat transfer through a porous medium (Garrett & Ban, 2011;Yun & Santamarina, 2007). These predictive models were critically reviewed elsewhere (Abdulagatova et al., 2009;Dong et al., 2015;Zhang & Wang, 2017).
Available predictive models are based on either based on empirical fits to experimental data or conceptualization of the multiphase porous system as a certain combination of series/parallel solid, air, and/or water blocks in a representative elementary volume. These models cannot realistically account for pore-scale mechanisms (fluid-fluid and rock-fluid interactions, pore-scale habit of each coexisting phase, and changes in the thermal state) controlling the heat transfer. Therefore, they need adjusting parameters whose values are usually a priori unknown, making the models dependent on the experimental data, and may accordingly restrict their applications (Sadeghi et al., 2018). In fact, a more comprehensive insight into these pore-scale mechanisms is required for further development of predictive models. Efforts have been made to account for 10.1029/2020WR027885

Water Resources Research
VASHEGHANI FARAHANI ET AL. these mechanisms by developing predictive models having physically meaningful parameters. One popular approach is to implement the pore-scale mechanisms into the unit cell. There are a number of models developed for predicting ETC of dry particulate beds in which a contact equation is employed to account for TCR (Bahrami et al., 2006;Dai et al., 2019;Garrett & Ban, 2011;Mo & Ban, 2017;Slavin et al., 2000;Weidenfeld et al., 2004). Some other unit cell-based models take into consideration the effect of the water content, hence applicable for unsaturated soils (Chu et al., 2019;Corasaniti & Gori, 2017;Ewing & Horton, 2007;Gori & Corasaniti, 2013;Likos, 2015). Use of such models is computationally affordable; however, their accuracy and application strongly depend on how they incorporate the pore-scale associated phenomena. More sophisticated approaches such as discrete element method and boundary element method have been also adopted recently to study the heat transfer in soil and establish more detailed predictive models (El Shamy et al., 2013;Feng et al., 2008;Yun & Evans, 2010;Zhou et al., 2007).
In this study, we report experimental measurements of ETC of simulated partially saturated sediments. Four sets of experiment were conducted to investigate the effect of pore pressure, overburden pressure, temperature, and degree of water/ice saturation on the pore-scale associated phenomena controlling heat conduction in a multiphase porous system. We also developed a numerical model for prediction of ETC, where free-energy lattice Boltzmann model (LBM) is first utilized to set the coexisting phases of water and gas (N 2 ) at equilibrium in a cubic domain, composed of a 3-D standard unit cell surrounded by an annular region filled with the gas. ETC of the cubic domain is then calculated via applying a space renormalization method on the system at equilibrium conditions. ETC of the standard unit cell is simply calculated by deducting the thermal conductivity of the annular region from ETC of the cubic domain based on the weighted averaging method. The model is able to account for several pore-scale mechanisms simultaneously, having more accurate predictions than the other existing models available in the literature.

Background
Heat transfer in granular materials is controlled by several particle-level mechanisms: (1) conduction within the mineral, (2) particle-particle conduction through the contact area, (3) particle-fluid-particle conduction, which is of particular relevance in partially saturated sediments, (4) radiation at interparticle contacts, (5) particle-fluid-particle conduction, (6) conduction in pore fluid, (7) convection in pore fluid, and (8) radiation from the particle surface into the surrounding medium (Yun & Santamarina, 2007). Figure 1 schematically illustrates these mechanisms for a cylindrical standard unit cell. Heat transfer could be considered conductive at naturally occurring porous sediments, particularly in fine sands and clays at low temperatures given negligible relative contribution of free convective and radiative mechanisms (Waite et al., 2009).
There are two types of thermal conductivity measurements including steady state methods and unsteady state methods. In this work, we employed transient hot wire method, which is an unsteady state procedure extensively used to measure the thermal conductivity of intact and reconstituted porous substances in a wide range of temperatures. In this method, thermal conductivity is obtained by a variation of the line source test method using a needle probe, which is required to have a large length to diameter ratio to comply with the line-source solution assumption, that is, an infinitely long, infinitely thin heating source (ASTM D5334-14, 2014). The needle probe consists of a heating wire, representing a perfect line source and a temperature sensor for measuring the temperature. The probe is inserted into the specimen, certain current and voltage are applied to the probe, and the temperature rise with time is recorded over a period of time. Thermal conductivity is obtained via analyzing the temperature time series data during the heating cycle (ASTM D5334-14, 2014). (1) particle conduction, (2) contact conduction, (3) particle-fluid-particle conduction, (4) particle-particle radiation, (5) particle-fluid conduction, (6) pore fluid conduction, (7) pore fluid convection, and (8) radiation into the surrounding medium.

Water Resources Research
VASHEGHANI FARAHANI ET AL.
In contrast with steady state methods, transient hot wire method could be conducted with substantial less required time for measurements (Abdel-Wahed et al., 1978;Jongwon et al., 2019). Moreover, the sample size is not critical as long as the diameter of the medium surrounding the needle is 10 times the diameter of the thermal needle probe (ASTM D5334-14, 2014). At a constant amount of heat flux, the temperature response of an ideal (zero mass) heater over a period of time could be expressed by Equation 1: where t, ΔT, Q, r, κ, λ, Ei, and t 1 stand for the measuring time, temperature change, applied heat per unit length of heater, distance from the heated needle, thermal diffusivity, thermal conductivity, exponential integral, and heating time, respectively (ASTM D5334-14, 2014). The exponential integral in Equation 1 could be approximated by the most significant terms of its series expansion, resulting in Equation 2: Therefore, λ could be calculated from the slope of the best fit line relating ΔT to ln(t). When analyzing the data, care must be taken to exclude the early and late portions of the test, representing transient conditions where the measurements were strongly influenced by terms ignored in Equation 2 and boundary conditions where the measurements were affected by the external boundary of the system, respectively.

Experimental Setup
ETC measurements were conducted using a 316 stainless steel cylindrical cell setup, which is schematically illustrated in Figure 2. The setup consists of a high-pressure stainless steel cell with an inner diameter of   -14, 2014). The setup can be placed horizontally or vertically on a pivot. An integral cooling jacket was fitted around the cell and connected to a cryostat (Grant LTC) to maintain the desired isothermal conditions during measurements. A calibrated platinum-resistance thermometer (Pt100, supplied by TC Ltd.) with an accuracy of better than ±0.1 K was also coated in the cell to measure its temperature. A movable piston driven by hydraulic pressure was placed at one end of the cell, used for compacting sediments at a desired overburden pressure. The overburden pressure was achieved by fluid injection behind the piston using a Quizix pump with dual cylinders (SP-5200, Vindum Engineering Inc., United States) and measured using a calibrated Druck pressure transducer with an accuracy of ±0.05 MPa. A linear variable differential transmitter is mounted on the piston rod to determine its position, which enables us to compute the instantaneous sample volume, hence the porosity of the system. There were two ports at the fixed endcap, one for the gas introduction and the other connected to another calibrated Druck pressure transducer for measuring the pore pressure. It should be noted that we inserted a fine mesh screen in the discharge point of each port to avoid intrusion of the sand particles and baffle the gas stream during injection to retain the initial homogeneous water distribution. The pressure and temperature data together with the piston movement were continuously monitored at regular time intervals, collected using a data acquisition module, and recorded on a PC with a LabVIEW software (National Instruments).
The thermal conductivity measurements were made using a calibrated purpose-built needle probe (TP08, Hukseflux, Netherlands), which can accurately measure the thermal conductivity at elevated pressures in compliance with the ASTM D5334-14 standard. In TP08, the reference junction of the thermocouple is located in the base. Hence, before loading the cell with the test specimen, the probe is mounted on the fixed endcap and its base fully covered with an insulation material to minimize the influence of heat exchange with the laboratory environment on the thermal conductivity measurements. The thermal conductivity data were all recorded in a measurement and control unit connected to the PC.

Materials
Research-grade N 2 (purity 99.995 vol.%) was supplied by BOC Ltd. Deionized water generated by an integral water purification system (ELGA DV 25) was used throughout the experiments. The simulated sediment used in this study was a well-characterized silica sand from Fife, Scotland (grain density: 2.64 g/cm 3 , mean size: 257 μm, and specific area: 0.059 m 2 /cm 3 ); a detailed analysis of which could be found elsewhere .

Procedure 3.3.1. Sample Preparation and Configuration
Each test specimen was made by mixing the silica sand and deionized water at a predetermined mass ratio with respect to the desired water saturation. The mixture was loaded within the cell that was washed, dried, positioned vertically, and equipped with the thermal conductivity probe instrumentation. While loading, the specimen was regularly tamped using a purpose-built pestle to attain uniform compaction and ensure full contact with the probe. After the specimen reached the desired height, the other endcap was installed, and the cell was rotated back to the horizontal position. The bath temperature was also set to the target temperature, and the cell was vacuumed. A hydraulic pressure was then applied behind the piston using the Quizix pump to further compact the sediment at the desired overburden pressure. The system was left for 24 hr at this condition to reach thermal equilibrium. N 2 was then injected into the cell until the pore pressure reaches the desired value while maintaining the effective overburden pressure. The system was left at this condition until the pressure and temperature of the system become stable and no piston movement is observed. Finally, the specimen underwent the ETC measurement.

Thermal Conductivity Measurements
Four sets of experiment were carried out in order to explore the effect of the system temperature (T), pore pressure (P Pore ), effective overburden pressure (P OB, eff. ), and water content (θ w ) on ETC of sediments. The experimental sets are summarized in Table 1. As observed, for each individual parameter, a base value was considered: T ¼ 0.5°C, P Pore ¼ 3.45 MPa, θ w ¼ 12.6 wt.% (S w ≈ 55%), and P OB, eff. ¼ 3.45 MPa. Each experimental set was carried out by changing two parameters while keeping the other two at their base 10.1029/2020WR027885

Water Resources Research
values. This enabled us to investigate the effect of each individual parameter on ETC as well as capture their interdependencies.
Most of ETC measurements were repeated for the elimination of the experimental data contingencies. Moreover, at frozen conditions, the total heat to the specimen was reduced to impose a smaller temperature gradient in order to minimize the possibility of fluid convection and phase change around the needle probe. After completion of ETC measurements, three sediment samples were taken from different sections in the cell and their water saturation was measured to check the homogeneity of the system. ETC was then obtained by analyzing the experimental data according to Equation 2.

Numerical Modeling
A numerical model was developed to predict ETC of partially saturated sediments under applied overburden pressure. The model uses a 3-D standard unit cell containing two monosized touching hemispheres in contact with the coexisting fluids, that is, water/ice and gas (N 2 ) at equilibrium conditions. The whole computational procedure could be divided into three steps:

Preprocessing
To begin with, a 3-D standard unit cell with 2r by 2r square cross-section area (r stands for the hemisphere radius) and h in height, containing two perfectly round hemispheres with a flat Hertzian contact interface, is constructed ( Figure 1). In this study, the simple cubic packing configuration was identified to be appropriate based on the porosity range of our specimens (Siu & Lee, 2000). The contact radius, r c , is calculated by the Hertz contact equation (Johnson, 1987): where F, v, and E are the contact force, Poisson's ratio, and Young's modulus. This makes it possible to account for the heat flow through the sphere-contact-sphere and accordingly incorporate TCR as a function of the effective overburden pressure, which could be considered as a distinguishing feature of our numerical model. According to Figure 1, the half spheres are in contact with static fluids filling the rest of unit cell. We employed free-energy LBM to set the two phases of water and gas (N 2 ). LBM is a computational fluid dynamics solver based on mesoscopic kinetic equations, which could be utilized to study fluid flow problems (Krüger et al., 2017). LBM has proven to be a powerful tool for solving problems at different length and time scales owing to its advantages in numerical simulation of single-phase and multiphase fluid flows in pore scale (Bakhshian et al., 2019;Chen et al., 2019;Norouzi et al. 2019;Zhang et al., 2019), particularly in complex geometries such as porous media (Bakhshian & Sahimi, 2016;Vasheghani Farahani et al., 2019;Vora & Dugan, 2019). Free-energy LBM as an option for simulating multiphase and multicomponent fluid flow is able to appropriately account for fluid-fluid interfacial tensions and solid surface contact angles. The free-energy LBM algorithm is implemented in OpenLB (http://www.openlb.net), an Open Source numerical framework for lattice Boltzmann simulations (Krause et al., 2019). Here, we briefly describe the free-energy LBM for a binary system in which the fluid-fluid surface tension and fluid-solid surface interactions (contact angle) could be derived and independently controlled. Further details of the algorithm implemented in OpenLB could be referred in the literature (Sadullah et al., 2018;Semprebon et al., 2016;Wöhrwag et al., 2018). A free-energy functional that models a system with two fluid components with concentration fractions of C 1 and C 2 (C 1 + C 2 ¼ 1) could be expressed by (Semprebon et al., 2016): where Ω is the system volume and α and κ are two tuning parameters for the interfacial width and surface tension. The interfacial width is expressed by α (typically chosen to be 1.0) and the surface tension (γ 12 ) by (Sadullah et al., 2018) The h i parameters in Equation 4 make it possible to quantify the fluid-solid surface energies, hence the contact angle of fluid 1 on a solid surface in the presence of fluid 2 (θ 12 ), which could be given by (Sadullah et al., 2018) where γ si is the surface tension between the fluid phase i and solid surface, which essentially includes contributions from the both phases, i and j, and could be expressed by (Sadullah et al., 2018) Therefore, the contact angle could be given by the following expression as a function of κ 1 , κ 2 , h 1 , and h 2 (Semprebon et al., 2016): Variable transformation from C 1 and C 2 to two equivalent order parameters of ρ ¼ C 1 + C 2 and ϕ ¼ C 1 − C 2 could be applied to satisfy the hard constraint of C 1 + C 2 ¼ 1. The continuum equations of motion of the system are the continuity, Navier-Stokes, and Cahn-Hilliard equations (Wöhrwag et al., 2018): where v, P, and η are the fluid velocity, fluid dynamic viscosity (which generally depends on the local order parameter ϕ), and pressure tensor, respectively. The thermodynamic properties of the system described in the free-energy model (Equation 4) enter the equations of motion via the chemical potential (μ ϕ ) and the pressure tensor. The lattice Boltzmann algorithm described by Semprebon et al. (2016) can be implemented to solve the equations of motion. For a binary system, two distribution functions, f i (x,t) and g i (x,t), corresponding to the density of the fluid ρ and the order parameter ϕ are required. The physical variables could be related to the distribution functions by The quantities e iβ correspond to the standard lattice velocities in the lattice Boltzmann (LB) method. The lattice structure we utilized here is D3Q19 with 19 velocities in three dimensions. A standard BGK (Bhatnagar-Gross-Krook) collision operator could be used for the collision step: 10.1029/2020WR027885

Water Resources Research
For the steaming step, we have where Δt stands for the lattice time step and f i eq (x, t) and g i eq (x, t) are the local distribution functions. The relaxation times τ and τ ϕ are related to the transport coefficient in the hydrodynamic equations, η and M ϕ : in which Γ ϕ is a tunable parameter that appears in the equilibrium distribution. The simulation parameters are tabulated in Table 2.
The LB simulation starts with setting up the geometry of the unit cell, followed by creating and initializing mesh with an optimum resolution in terms of the stabilization of the results and computational cost. This consists of classifying them with material numbers according to their type (1: fluid nodes, 2: solid boundary nodes, and 0: inner solid nodes). Fluid nodes are divided into two lattices, and their initial spatial distribution is specified based on the desired saturation. The lattice dynamics are then set for all nodes based on their material number, and the collision model along with the boundary behavior is characterized. The next step is the main LB loop where the collision and streaming are conducted, followed by the lattice coupling and boundary condition treatment in order to account for the fluid-fluid and solid-fluid interactions, respectively. The spatial distribution of the phases at equilibrium conditions are finally extracted to be used for the next step. If necessary, the volume treatment due to the expansion of the liquid phase at subzero temperatures is also performed.

Processing
As previously mentioned, we apply a 3-D space renormalization technique to predict ETC of partially saturated sediments. Similar technique has been applied to calculate ETC of tetrahydrofuran hydrate-bearing sands based on 2-D pictures taken by a digital electron microscope and shown to have strong stability and fault tolerance (Huang & Fan, 2005). In this method, a given volume is partitioned into several blocks, each with its own initial property. Based on certain rules, the adjoining eight blocks are considered as a renormalization group to form a new block. The process is repeated until the final block is obtained. We consider the following assumptions to obtain ETC: (i) heat flows in one direction (see Figure 3), (ii) thermal conductivity of blocks containing two phases is averaged based on their volume fractions, (iii) conductive heat resistance is analogous to electrical resistance, and (iv) heat transfer via convection and radiation mechanisms is ignored given their negligible contribution compared with the conduction mechanism (see Cortes et al., 2009).
According to Huang and Fan (2005), a four-block renormalization group (2-D) could be treated as a combination of eight heat conductors (or resistors) on the basis of the analogy between conductive heat transfer and electrical conduction. Similarly, as can be seen in Figure 3, the analogy between conductive heat transfer and electrical conduction enables us to treat each eight-block renormalization group as a combination of 20 heat conductors. The resistances of these heat conductors are

Water Resources Research
VASHEGHANI FARAHANI ET AL.
where k stands for thermal conductivity and subscripts refer to the adjoining blocks, as shown in Figure 3. The equivalent heat resistance and accordingly the equivalent thermal conductivity could be obtained for each renormalization group by applying Kirchhoff's current law for each individual node (see supporting information). Therefore, this is considered as the core rule of the renormalization process. In addition, it could be concluded that "n" stages of the renormalization is required to be conducted on a cubic domain with 8 n (2 n × 2 n × 2 n ) blocks, which necessitates the LB simulation to be performed on a cubic domain with the size of (2 n ) 3 . In fact, the parameter n is a nonnegative integer that determines the number of stages required for the renormalization process and essentially controls the resolution.
The LB simulation results at the preprocessing step consist of the status of each voxel (grain, water, or N 2 ) at equilibrium conditions. Therefore, the simulation domain is already partitioned and labeled. In this step, we assign the thermal conductivity of each individual voxel according to its label and implement the renormalization program in order to obtain the equivalent heat resistance. The thermal conductivity of the voxels located in the fluid-fluid interface is calculated by taking an arithmetic saturation-based average of thermal conductivities of water and N 2 .
It should be noted that we used the thermal conductivity data available in our in-house software for water/ice and nitrogen (Hassanpouryouzband et al., 2018;Hassanpouryouzband, Farahani, et al., 2019). The intrinsic thermal conductivity of sand particles was also obtained by running the numerical model for a case with parameters (T, P Pore , P OB, eff. , and θ w ) set as base values and comparing the predicted ETC with the experimental result, similar to the approach applied by Huang and Fan (2005). By employing this method, the thermal conductivity of sand was obtained as 6.86 W/m.K. The expansion of water due to transformation from liquid to solid (ice) at temperatures below the freezing point is accommodated as follows:

Water Resources Research
VASHEGHANI FARAHANI ET AL.
When performing the renormalization process, the change in dimensions of the renormalized block is taken into consideration as ETC is an intensive property, while resistance is extensive (Huang & Fan, 2005).

Postprocessing
Given that the LB simulation is required to be performed on a cubic domain and the renormalization process on (2 n ) 3 blocks (Huang & Fan, 2005), we calculated ETC of the cubic domain (k Cube ) comprised of the 3-D standard unit cell and an annular region (see Figure 4) and used weighted averaging method to obtain ETC of the unit cell (k Cell ), as illustrated in Equation 16. Note that the annular region essentially contains gas.
where V is the volume.

Results and Discussion
In this section, we present the experimental results along with the numerical model predictions. Moreover, we plotted our model predictions together with a number of existing models in the literature against the experimental results. This helps us to understand the accuracy of the newly developed model in comparison with the other models. Figure 5 shows the measured ETC of partially saturated sediment samples (12.6 wt.%) against the effective overburden pressure at three different pore pressures and a constant temperature (0.5°C). General trends at different pore pressures have been also added to assist with when discussing the thermal conduction characteristics. As expected, increasing the effective overburden pressure results in a higher ETC. This could be attributed to two main factors. The first is the higher particle-particle conduction through the contact area (lower TCR). Increasing the effective overburden pressure results in a higher contact area, as can be simply confirmed by Hertz contact equation (the inset in Figure 5). Given the intrinsic thermal conductivity of sand particles is considerably higher than that of the other constituents, a lower TCR and greater contribution from the particle-particle conduction to ETC is expected when the effective overburden pressure increases. This contribution essentially depends on the contact area, hence more influential at lower effective overburden pressure values. Increasing the effective overburden pressure may also increase the coordination number, which results in a higher number of particles in contact with the particle under consideration. The other factor to be considered is the higher particle-fluid-particle conduction. Increasing the effective overburden pressure results in redistribution of pore fluids enclosing the contact area. Assuming the particles to undergo elastic deformation, when the effective overburden pressure increases, the pore space around the contact areas becomes smaller with a higher capillary pressure. As a consequence, if there is any gas in these regions, it will be displaced by water as the wetting phase. Since these areas contribute critically to the heat transfer, increasing the effective overburden pressure may also facilitate the heat conduction through these critical pathways.

Dependence on Overburden and Pore Pressures
As can also be seen in Figure 5, ETC is directly proportional to the pore pressure of the system. Several pore-scale associated heat transfer mechanisms could cause this behavior. The first mechanism could be higher heat conduction and convection in the gas phase. Increasing the pore pressure results in a higher intrinsic thermal conductivity for the gas phase. At 0.5°C, the intrinsic value for the thermal conductivity of N 2 at 13.79 MPa is around 24% higher than 2.76 MPa. Apart from that, higher pore pressure can improve the convective heat transfer mechanism in the gas phase. The other mechanism is higher heat transfer at the fluid-fluid interface. In a multiphase system, increasing the pore pressure enhances the mass transfer at the fluid-fluid interface. This mass transfer is always associated with an energy transfer. Therefore, it is expected to have more contribution in the heat transfer from the gas-water interface at elevated pore pressures. The other interesting mechanism could be the redistribution of the water toward the smaller pores. Since gas is the nonwetting phase in the system, it prefers to occupy the largest pores whereby influencing the pore water network. In fact, when pore pressure increases, gas pushes water from the large pores to the small ones where the capillary pressure is higher. These small pores remarkably contribute to the heat transfer; hence, the presence of water with higher thermal conductivity in these pores is conducive toward higher ETC. Figure 6 illustrates the measured ETC of partially saturated sediment samples (12.6 wt.%) against the effective overburden pressure at different temperatures and constant pore pressure (3.45 MPa). As observed, the increase in the effective overburden pressure at both unfrozen and frozen conditions resulted in a higher ETC. Potential influential mechanisms are already discussed in section 5.1.1. Figure 5. ETC versus the effective overburden pressure (P OB, eff. ) at three different pore pressures (P Pore ); T ¼ 0.5°C and θ w ¼ 12.6 wt.%. Contact radius ratio is also plotted according to Hertz contact equation against the effective overburden pressure as an inset.

Water Resources Research
ETC values measured at frozen conditions are markedly higher in comparison with unfrozen conditions since on one hand the intrinsic thermal conductivity of ice is almost 4 times more than that of water and on the other hand the transformation of water to ice is associated with volume expansion in the water phase. Our measurements show that the transformation of water to ice resulted in the gas volumetric saturation (S g ) to decrease around 7% relative to its initial value. Hence, the results are presented in Figure 6 separately, according to the thermal state of the system, to allow for a clearer investigation.
Increasing the system temperature at unfrozen conditions results in a higher ETC (Figure 6a). This could be due to higher heat conduction and convection in the pore fluids. The intrinsic thermal conductivity of each constituent is directly proportional to the system temperature at the temperature range we dealt with. Therefore, a higher temperature may increase the contribution of heat conduction from each component in ETC. At 3.45 MPa, the intrinsic value for the thermal conductivity of N 2 at 10°C is not significantly different from 0.5°C. However, despite the pore pressure, temperature change also can influence the intrinsic thermal conductivity of water and sand particles. Increasing the system temperature could also improve the convective heat transfer mechanism in the pore fluids (i.e., gas and water), resulting in a higher ETC. The other factor that could play a role is the higher heat transfer at the fluid-fluid interface. Similar to what discussed earlier for the effect of the pore pressure, increasing the system temperature results in a higher mass and energy transfer in the gas-water interface and higher ETC, accordingly.
As shown in Figure 6b, further reducing the system temperature at frozen conditions resulted in a higher ETC, a behavior against what is observed at unfrozen conditions. This behavior could be attributed to the presence of unfrozen water at the critical pathways. It has been shown that the unfrozen water content at −5°C is higher than −10°C (Istomin et al., 2017). This unfrozen water with a significantly lower thermal conductivity than ice still contributes critically to the heat transfer because ice formation starts at larger pores first and pore fluid migrates toward the nucleating ice due to cryosuction (Hohmann, 1997). Therefore, at −5°C, ETC is more influenced by the unfrozen water content. It should also be noted that the transformation of unfrozen water to ice could result in further volume expansion and contribution in ETC. Moreover, heat conduction in ice is higher at lower temperatures. It has been shown that the intrinsic thermal conductivity of ice is negatively correlated with temperature (Slack, 1980). Therefore, it could be expected that the ice contribution in the heat transfer would be higher at −10°C relative to −5°C, resulting in a higher ETC.
As observed in Figure 6, enhancement of ETC as a result of increasing the effective overburden pressure is typically smaller at frozen conditions compared with unfrozen conditions. The average enhancement rates of ETC at unfrozen and frozen conditions are~0.0316 and 0.0146 W/(m.K.MPa), respectively. This behavior could be due to the fact that at frozen conditions, ice crystals could act as coarse grains pushing apart the sediment grains, hence enlarging the sediment pores, a phenomenon known as ice-forced heave, and eventually influence ETC negatively. A similar phenomenon caused by gas hydrates has been already observed to be markedly influential in permeation characteristics of hydrate-bearing sediments . Moreover, it can be seen that ETC enhancement is higher for an individual effective overburden pressure at unfrozen conditions. This behavior can clearly highlight the contribution of the heat conduction and convection in the pore fluids along with the heat transfer at the fluid-fluid interface in ETC at unfrozen conditions. Figure 7 shows the measured ETC of partially saturated sediment samples (12.6 wt.%) against the pore pressure at different temperatures and a constant effective overburden pressure (3.45 MPa). As expected, ETC increases as the pore pressure increases for both unfrozen and frozen sediments. We already discussed potential influential mechanisms in section 5.1.1. Dependence on the system temperature is also observed to be similar to what discussed in section 5.1.2. Here, the contribution of the heat conduction and convection in the pore fluids in conjunction with the heat transfer at the fluid-fluid interface to ETC is observed once again to be more influential at unfrozen conditions. We should note that care was taken not to be within the hydrate phase boundary of N 2 when performing the ETC measurements. Hence, it is observed that the ETC measurement was only conducted outside N 2 hydrate stability zones at both unfrozen and frozen conditions. Figure 8 illustrates ETC of partially saturated sediment samples against the pore pressure and volumetric water saturation (S w ) at a constant temperature (0.5°C) and effective overburden pressure (3.45 MPa). According to Figure 8a, ETC could be positively correlated with the pore pressure, which can be interpreted in the similar mechanisms previously discussed in section 5.1.1. We should also note that for fully saturated sediments we carried out the ETC measurements for more elevated pore pressures up to 20.68 MPa since water is incompressible and injection of a small amount of water could result in a drastic increase in the pore pressure.

Dependence on Water Content
Figure 8b enables us for further investigation regarding the effect of the water content on ETC. At a typical pore pressure, it can be seen that ETC increases considerably once the sediment experiences a small increase in the water content, gaining almost 45% of the heat transfer capability at fully saturated conditions during wetting from dry state to around 10% volumetric water saturation. This rapid increase in ETC is because the sediment experiences the pendular regime at which water forms a thin film around sediment particles and builds water bridges at the contacts, which results in a significant improvement in the connectivity of heat transfer pathways throughout the system (Dong et al., 2015). In fact, our sediment samples had very limited hydration stage at the early stage of wetting, as there was no fines content (clay and silt) in the system (Dong et al., 2015). This behavior is in accordance with the correlation proposed for the critical water content (θ c ) as a function of the clay content by Sadeghi et al. (2018). According to Figure 8b, the pendular regime for our sediment samples is ranging from S w ¼ 0 to 15%. At higher saturations (up to 70%), prevalence of the funicular regime is observed where ETC still increases gradually mainly due to the buildup of the pore water network. The effect of the pore pressure on ETC can also be observed to be more obvious at this stage given that the establishment of the pore water network is influenced by the pore pressure of the system, as

10.1029/2020WR027885
Water Resources Research discussed in section 5.1.1. At saturations higher than 70%, the sediment experiences the capillary regime at which the effect of the water content on ETC becomes relatively inconspicuous as increasing the water content does not alter the preferred heat flow pathways or their connectivity anymore. Figure 9 visualizes the cubic domain comprised of the standard unit cell and annular region. The optimum resolution of the numerical model, that is, the optimum renormalization rank in terms of the stabilization of the results and computational cost, was found to be 2 7 via running the model and processing the simulation results at five different resolutions (2 5 , 2 6 , 2 7 , 2 8 , and 2 9 ). This means that the model is required to undergo seven stages of the renormalization process to obtain corresponding ETC. Further details regarding the procedure followed to choose the optimum resolution could be found in the supporting information.

Numerical Model Predictions
As mentioned earlier, the unit cell contains two monosized touching hemispheres. The contact region set for the effective overburden pressure of 3.45 MPa is also shown in Figure 9a. It should be noted that the number of nodes in the contact region was determined according to the contact radius ratio (r c /r) with respect to the desired effective overburden pressure. Moreover, it can be seen that the solid boundary nodes (which lie at

10.1029/2020WR027885
Water Resources Research the solid-fluid interface) are visualized in red, a different color from the inner solid nodes (which are isolated hence inactive), as the boundary nodes should be treated differently in terms of the boundary conditions in the preprocessing step (the LB simulation) to eliminate unnecessary computations at inactive nodes (Sukop & Thorne, 2006). The no-slip boundary conditions were imposed for the outer walls of the domain (x, y, and z directions) by using the standard bounce-back method. Bounce-back boundaries with a controllable contact angle were also considered for the solid boundary nodes. Figures 9b and 9c also visualize the cubic domain from two different views, perpendicular to x and y axes, respectively. The porosity of the standard unit cell can be approximately obtained as 33%, which is quite close to the porosity of our sediment samples ranging in 38-40%.
The model was run to predict ETC with respect to the system temperature, pore pressure, effective overburden pressure, and water content at which we measured ETC. The predicted values are presented in Figure 10 in comparison with the experimental data. As could be seen, the model adequately accounts for the effect of temperature, pore pressure, overburden pressure, and water content when predicting ETC. Snapshots taken from the distribution of water at different saturations in equilibrium with gas are also presented in Figure 11 to help further explore how the model is able to capture pore-scale Figure 11. Water distribution at the equilibrium conditions in the cubic domain at volumetric saturations of (a) 8%, (b), 37%, (c) 57%, and (d) 89%. The rest of the cell (white space) is occupied by the gas phase. Water Resources Research mechanisms in predicting ETC. As observed, the model is able to consider the heat conduction through the particle-particle contact area specified according to the Hertz contact equation. Next, the particle-fluidparticle conduction is well taken into consideration, particularly at low water saturations (Figure 11a) where water as the wetting phase forms a thin film on the particles and builds capillary bridges enclosing the contacts, contributing significantly to the improvement of ETC. As discussed in section 5.1.1, the intrinsic thermal conductivity of the gas phase (N 2 ) is a function of the pressure, enabling the model to account for the effect of the pore pressure on ETC. Moreover, the fluid-fluid interfacial phenomena are accommodated at the preprocessing step where the volume fraction of water at the interface is determined via applying free-energy LBM. Average relative deviations (ARD) were also calculated using Equation 17 and provided in Table 3, evidencing an excellent agreement between the model predictions and experimental data.
At frozen conditions (Figures 10b and 10c), it is observed that the model predictions are all higher than the experimental data. This  Water Resources Research could be due to the procedure we followed to accommodate the transformation of water to ice in the processing step. In fact, the model neglects the unfrozen water content, which may critically contribute to ETC, and calculates the amount of the volume expansion according to Equation 15. Ice-forced heave is the other possible factor required to be further investigated experimentally.
We also compared the predictions with those of six existing models: Quadratic Parallel (Dong et al., 2015), Hashin-Shtrikman Lower Bound (HSL) (Dong et al., 2015), Hashin-Shtrikman Upper Bound (HSU) (Dong et al., 2015), Haigh (2012), He et al. (2017), and Sadeghi et al. (2018). When using these models to predict ETC, the instantaneous porosity and water saturation corresponding to each experimental data point were considered. Moreover, the volume expansion due to the transformation of water to ice was included when predicting ETC at frozen conditions. Figure 12 exhibits the comparison of the different modeling predictions of ETC with the experimental data. As expected, all predictions are within the Hashin-Shtrikman bounds. It can be seen that the other models except HSU and quadratic parallel underestimate ETC; some even fail to make predictions with a relative error of less than 30%. This could be due to the fact that the predictive models proposed so far are mixing models relying on simple mixing laws such as series, parallel, and quadratic parallel models, empirical models strongly depending on the experimental data such as the model developed by He et al. (2017), or mathematical models that were implicitly/explicitly developed such as Hashin-Shtrikman's model, Haigh model (Haigh, 2012), and GD model (Ghanbarian & Daigle, 2016) and its explicit form developed by Sadeghi et al. (2018). Therefore, they would not able to capture all complex pore-scale associated phenomena contributing to ETC.

Conclusion
In this work, we conducted a mechanistic study on the heat transfer in partially saturated sediments at unfrozen and frozen conditions. Four sets of experiment were carried out, and the effects of a variety of parameters on heat transfer mechanisms were investigated, including pore pressure, overburden pressure, temperature, and degree of water/ice saturation. A pore-scale numerical model was also developed based on the deployment of free-energy LBM and a space renormalization method for the prediction of ETC. The model prediction is in good agreement with the experimental data with an average RMSE of less than 0.1 W/m.K. The experimental and modeling results reveal that the heat transfer in porous media is complex phenomena and controlled by several key factors.
1. Heat transfer in porous media is largely affected by water content, packing structure, and wettability characteristics of grains as they influence the capillary pressure and preferential occupancy and hence the water distribution in pores. 2. Coordination number and physical contact among sediment particles substantially affect the heat transfer. In fact, ETC is remarkably influenced by the particle-particle conduction through the contact regions. The particle-fluid-particle conduction through the pathways enclosing the contact regions is the next important heat transfer mechanism which contributes to ETC, particularly at low water saturations. 3. Heat conduction and convection in the pore fluids should be taken into account to further improve the model performance, especially at elevated pore pressures and temperatures. 4. ETC at frozen conditions was affected by the unfrozen water content still enclosing the contact regions.
Ice-forced heave may also play a role in heat transfer at frozen condition. However, this phenomenon is to be further investigated. 5. Performance of the pore-scale model could be further improved for different packing structures such as body centered cubic and face centered cubic, which requires employing different contact equations appropriate with the packing structure and contact angles, and by coupled discrete element method-computational fluid dynamics methods. Such coupled methods could also assist with when investigating the effect of particle shape and size distribution. More accurate predictions could be 10.1029/2020WR027885

Water Resources Research
achieved at frozen conditions by accounting for the coexistence of ice and unfrozen water together with their pore-scale distribution and interfacial effects. This is of particular importance at subzero temperatures in the vicinity of the freezing point where the amount of unfrozen water is not negligible, and the ice crystals are floating in unfrozen water (pore-filling habit).