1 Introduction

Thermal conductivity of porous rocks depends on a large variety of proper to rock parameters as well as external influences. Thus, it can generate difficulties in determining accurate thermal behavior of the rock. The rock parameters which influence the thermal conductivity are principally the porosity, the microstructure [1] and the mineral composition. However, these parameters, in turn, can be impacted by external influences such as temperature and pressure.

An accurate determining of the thermal conductivity is crucial in oil and gas engineering or in geothermal application. For example, during thermal EOR or geothermal application, the porosity and/or the microstructure of the sedimentary rocks can vary due to the increase of temperature and pressure, and this modification must be quantified to be accounted for the thermal behavior of rocks.

Many efforts have been done to estimate the thermal conductivity of sedimentary rocks in parallel to the experimental methods for its determination. These estimations have always been the subject of intensive studies, and a lot of data [2, 3] are obtained as well as models and methodologies to characterize the thermal conductivity of rocks [4–7]. Moreover, this type of estimation is well-known by other research communities. Indeed, we find the same formal analogy between Fourier, Ohm’s law, Darcy’s laws and thermal conductivity. For example, considering Darcy’s laws, the same problem is well-known and termed “upscaling” [8, 9] and consists in computation of the effective permeability considering a heterogeneous rock. Classically, the upscaling process can be related to percolation theory [10], which describes connectivity of objects within for example, a porous structure. We can also determine effects of this connectivity on macroscale properties such as thermal conductivity [11]. In particular, the fundamental contributions of Torquato who proposed strategies via rigorous microstructure-property relations [1, 12].

Finally, many technics are based on a porosity dependence and a link between the conceptual thermal conductivity of the non-porous rock, λ_R, and the thermal conductivity of the fluid saturated the porous rock λ_f. These technics are simple to implement especially when there is no precise information about the microstructure.

In order to predict accurately the thermal efficiency of the geothermal installation or the oil recovery of a thermal EOR process, such as, for example, Steam Assisted Gravity Drainage (SAGD), very often the engineers invoke numerical simulations. Numerous reservoir simulators [13–16] allow to estimate the thermal conductivity as function of porosity, but these solutions are often based on a mixing laws which are quite simplistic models.

The purpose of this paper is to propose a better methodology to predict a thermal conductivity of reservoir rocks depending on the porosity for a reservoir simulator which satisfy the following conditions

• easy to implement in a software, precise and accurate,

• it should have finite values to avoid numerical problems such as non-convergence,

• have ideally correct limits: $\underset{\phi \to 1}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{f}}$ and $\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$,

• have a limited number of parameters in order to minimize a required number of laboratory experiments.

First for all, we perform a review of many means or methods available in the literature, starting by the classical mixing laws and finishing by empirical and theoretical methods. For each method, the limits for φ → 1 and φ → 0 are examined and the corresponding parameters are detailed.

Then, we compare these methods to experimental data. Two experimental datasets are considered, water and air-saturated rock of the same properties for the both cases (sandstone). The purpose is to match these two datasets using the same model. Indeed, in reservoir simulation (oil and gas or geothermal), various fluids may flow in the same rock (oil, gas, CO₂, water, steam). Thus, it is crucial to propose a model which yields correct thermal conductivity independently on the fluids flowing in the pore volume.

2 Thermal conductivity of rock function of porosity: a review

The heat transfer in the reservoir is mainly dominated by conduction. Thermal conductivity is the property that quantifies the ability of a geological formation for heat transfer. There are three types of methods to define this conductivity of multicomponent systems, mixing laws, empirical models and theoretical models. The following denotations are used

• λ_R, the rock thermal conductivity,

• λ_f, the fluid thermal conductivity,

• φ, the rock porosity.

2.1 Mixing laws

The mixing laws are very common and are also used for Ohm’s law or Darcy’s laws effective parameters determination [8, 9]. But they don’t consider the structural characteristics of rocks and, therefore, their applications can be constrained.

First, consider the weighted arithmetic and harmonic means, reviewed by Zimmerman [5] which are also known as parallel, or linear, and series models, respectively. They correspond to the classical Wiener’s bound used to compute effective permittivity.

The weighted arithmetic or so-called parallel model is written as

$${\lambda }_{\mathrm{eff}}\left(\phi \right)=\phi {\lambda }_{\mathrm{f}}+\left(1-\phi \right){\lambda }_{\mathrm{R}}.$$(1)

This model is available in the most of the reservoir simulators such as Eclipse [14], Intersect [15] and Stars [16].

The weighted harmonic or so-called series model, is given by

$${\lambda }_{\mathrm{eff}}\left(\phi \right)={\left(\frac{\phi }{{\lambda }_{\mathrm{f}}}+\frac{1-\phi }{{\lambda }_{\mathrm{R}}}\right)}^{-1}.$$(2)

It can be noticed that the parallel model yields the highest values of thermal conductivity and the series model yields the lowest ones [4]. It was also shown by Wiener, that these two models represent upper and lower bounds on the effective conductivity [17].

Then, the mixing law is reviewed by Beck 18] as the weighted geometric mean. It has no physical background but seems to work better than the parallel and series means; it is written as

$${\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}{\left(\frac{{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{R}}}\right)}^{\phi }.$$(3)

This mean is also used by [5] and [18]. It can be noticed that this mean is also available in Stars [16]. In this reservoir simulator, the formula is given in the logarithm form,

$$\ln {\lambda }_{\mathrm{eff}}\left(\phi \right)=\ln {\lambda }_{\mathrm{R}}+\phi \ln \left(\frac{{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{R}}}\right)=\left(1-\phi \right)\ln {\lambda }_{\mathrm{R}}+\phi \ln {\lambda }_{\mathrm{f}}.$$(4)

There is another mixing law, proposed by Hashin and Shtrikman [19], which yields the bounds (HS bounds) always tighter than the Wiener bounds

$${\lambda }_{\mathrm{eff}, \min }\left(\phi \right)={\lambda }_{\mathrm{f}}+\frac{3{\lambda }_{\mathrm{f}}\left({\lambda }_{\mathrm{R}}-{\lambda }_{\mathrm{f}}\right)\left(1-\phi \right)}{3{\lambda }_{\mathrm{f}}+\left({\lambda }_{\mathrm{R}}-{\lambda }_{\mathrm{f}}\right)\phi },$$(5)

$${\lambda }_{\mathrm{eff}, \max }\left(\phi \right)={\lambda }_R+\frac{3{\lambda }_R\left({\lambda }_f-{\lambda }_R\right)\phi }{3{\lambda }_R+\left({\lambda }_f-{\lambda }_R\right)\left(1-\phi \right)}.$$(6)

In [20], the mean of both bounds is often used as best approximation of rock bulk thermal conductivity:

$${\lambda }_{\mathrm{eff} }\left(\phi \right)=\frac{1}{2}\left({\lambda }_{\mathrm{eff}, \min }\left(\phi \right)+{\lambda }_{\mathrm{eff}, \max }\left(\phi \right)\right).$$(7)

Finally, a dispersive or extended Maxwell model proposed by Zeb et al. [21] seems to have a solid physical basis [18]. The thermal conductivity of the porous rock is written as

$${\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}\frac{\left(\frac{{2\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}+1\right)-2\phi \left(\frac{{\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}-1\right)}{\left(\frac{{2\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}+1\right)+\phi \left(\frac{{\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}-1\right)}.$$(8)

However, it can be easily demonstrated that this formulation corresponds to the upper Hashin–Shtrikman bound λ_{eff, max}(φ). Therefore, in this paper, only HS bounds model is referred.

To conclude, all presented mixing laws are tend to expected values at the porosity limits, i.e., $\underset{\phi \to 1}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{f}}$ and $\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$. Moreover, no additional parameter is required to set these models. Notice that Hashin and Shtrikman bounds are classical and appear also in the computation of effective permeability of heterogeneous porous media [22].

2.2 Empirical models

This section addresses a review of the existing empirical models to determine the thermal conductivity of the porous rock. These models determine the thermal conductivity through the application of a regression technics. Generally, they are limited to a particular type of rocks. The purpose is to demonstrate the link in-between different models, empirical or mixing laws, as well as verify that the conductivity values at the porosity limits satisfy the expected values, λ_f for φ → 1 and λ_R for φ → 0.

First, Asaad [23] proposes a thermal conductivity function obtained by the mean which is very similar to weighted geometric model

$${\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}{\left(\frac{{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{R}}}\right)}^{\mathrm{c\phi }},$$(9)

where c is an empirical exponent and its value can be fitted using experimental data. When c = 1, this model becomes identical to the weighted geometric mean model. Moreover, ∀c > 0 only one of the conditions at the porosity limits is satisfied

$$\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}.$$(10)

Then, in [24] and [7], Sugawara and Yoshizawa propose a model based on an adjustable parameter A

$${\lambda }_{\mathrm{eff}}\left(\phi \right)=\left(1-A\right){\lambda }_{\mathrm{R}}+A{\lambda }_{\mathrm{f}},$$(11)

where

$$A=\left[2^n{\left(2^n-1\right)}^{-1}\right]\left[1-{\left(1+\phi \right)}^{-n}\right].$$(12)

with n (> 0) is an empirical exponent which needs to be fitted using experimental data. Notice that if n = −1, one obtains the weighted arithmetic mean (parallel model). Moreover ∀n > 0,

$$\underset{\phi \to 1}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{f}}\hspace{1em}\mathrm{and}\hspace{1em}\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}.$$(13)

Zeb and Maqsood [4] used exponential decay trial to predict thermal conductivity of consolidated porous media at room temperature and normal pressure,

$${\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}{\mathrm{e}}^{-\mathrm{z\phi }\frac{{\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}},$$(14)

where z is an empirical exponent which can be fitted using experimental data of thermal conductivity λ_eff and corresponding values of φ and λ_R as proposed by Zeb and Maqsood [4]. This means that the thermal conductivity of the non-porous rock is not given and is not a constant. Notice that this conductivity, λ_R, is very difficult to obtain because the existence of a non-porous rock with the same mineral composition as a porous one is almost impossible. However, for a reservoir simulator, only correct limits and a limited number of undefined parameters are identified as most important.

Here, $\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$ is directly verified. For φ = 1, the required limit condition is satisfied only if $z=\frac{{\lambda }_f}{{\lambda }_R}\ln \frac{{\lambda }_f}{{\lambda }_R}$, which yields the weighted geometric mean.

A model proposed by Zeb et al. [21] for prediction of thermal conductivity of porous consolidated igneous rocks is written as

$${\lambda }_{\mathrm{eff}}\left(\phi \right)={\left(\frac{1}{{\lambda }_{\mathrm{R}}}+\frac{\mathrm{m\phi }}{{\lambda }_{\mathrm{f}}}\right)}^{-1}$$(15)

where m is an empirical coefficient and should be fitted using experimental data of thermal conductivity λ_eff and corresponding values of φ and λ_f. It means that the thermal conductivity of the fluids is not given and is not a constant.

For the zero porosity, the model verifies $\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$, as required. However, for φ = 1 the condition $\underset{\phi \to 1}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{f}}$ is fulfilled only if $m=\frac{{\lambda }_{\mathrm{R}}-{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{R}}}$, i.e. when the model degenerates toward the weighted harmonic mean. This could represent a problem for a reservoir simulator since λ_R is not easy to determine contrary to λ_f which can be measured.

Veerendra and Chaudhary model for porous consolidated materials [4, 25] for ${\lambda }_{\mathrm{R}}>{\lambda }_{\mathrm{f}}$ is written as

$${\lambda }_{\mathrm{eff}}\left(\phi \right)=\left(1-\phi \right){\lambda }_{\mathrm{H}}+\phi {\lambda }_{\mathrm{L}}-\psi {\left(\frac{{\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}\right)}^{1/3},$$(16)

where ${\lambda }_{\mathrm{H}}={\lambda }_{\mathrm{R}}{\mathrm{e}}^{\beta \phi }$ and ${\lambda }_{\mathrm{L}}={\lambda }_{\mathrm{f}}{\mathrm{e}}^{-\beta \left(1-\phi \right)}$ with $\beta =\left(\frac{{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{R}}}-1\right)$. The term $\psi {\left(\frac{{\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}\right)}^{1/3}$ is to account for a high thermal conductivity ratio $\frac{{\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}$, ψ is an empirical coefficient which is fitted using experimental data. According to Veerendra and Chaudhary, there is a corrective term, $\phi \left(1-\phi \right)\sqrt{{\lambda }_{\mathrm{L}}{\lambda }_{\mathrm{H}}}$, which may be either added or subtracted from the precedent equation. For the porosity limits, φ = 0 and φ = 1, the required conditions are verified if $\psi$ = 0.

2.3 Theoretical models

This section addresses a review of the existing theoretical models to determine the thermal conductivity of the porous rock. These models are based on physical consideration such as pore structures.

The first theoretical model to predict the thermal conductivity of fluid saturated rocks is the Krupiczka model [26]. It is a semi-empirical equation derived from the numerical calculation of heat transfer through a bundle of cylinders for the effective conductivity of a packed bed of spheres, therefore, it is indexed in theoretical models,

$$\frac{{\lambda }_{\mathrm{eff}}\left(\phi \right)}{{\lambda }_f}={\left(\frac{{\lambda }_R}{{\lambda }_f}\right)}^{a_1+a_2\log \phi +B \log \left(\frac{{\lambda }_R}{{\lambda }_f}\right)}.$$(17)

were the coefficients $a_1=0.28$, $a_2=-0.757$ and $B=-0.057$ are obtained for this particular case of packed bed of spheres [26]. It should be noted that the required condition at the zero porosity limit, ${\underset{\phi \to 0}{\lim }\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$, cannot be satisfied for this model.

The last considered here model is based on the concept of the porous rocks represented as a composite material with the spheroidal non-connected pores. This formulation represents some mathematical difficulties due to arbitrary aspect ratio of the spheroidal inclusions. The solution has been proposed by Fricke [27] for the electrical conductivity [5] and can be successfully applied to the thermal conductivity. Thus, for small porosity, Fricke showed that the effective conductivity can be given by

$$\frac{{\lambda }_{\mathrm{eff}}\left(\phi \right)}{{\lambda }_{R }}=1-\beta \phi$$(18)

with

$$\beta =\frac{1-r}{3}\left(\frac{4}{2+\left(r-1\right)M}+\frac{1}{1+\left(r-1\right)\left(1-M\right)}\right),$$(19)

where $r=\frac{{\lambda }_{\mathrm{f}}}{{\lambda }_{\mathrm{R}}}$ and M is a factor that depends on the pore aspect ratio $\alpha$. For the two regimes of oblate (α < 1) and prolate (α > 1) spheroids, this factor is given by

$$\begin{array}{cc}\alpha <1:& \begin{array}{cc}M=\frac{\left(2\mathrm{\theta }-\sin 2\mathrm{\theta }\right)}{2\tan \mathrm{\theta }{\left(\sin \mathrm{\theta }\right)}^2},& \mathrm{\theta }={\cos }^{-1}\alpha \end{array}\\ \alpha >1:& \begin{array}{cc}M=\frac{1}{{\left(\sin \theta \right)}^2}-\frac{{\left(\cos \theta \right)}^2}{2{\left(\sin \theta \right)}^3}\ln \frac{1+\sin \theta }{1-\sin \theta },& \theta ={\cos }^{-1}\frac{1}{\alpha }\end{array}\end{array}.$$(20)

There are three limiting cases which are distinguished by Zimmerman [5] as thin cracks for α → 0, spherical pores for α = 0 and needle-like pores for α → ∞, with

$$\begin{array}{cc}\alpha \to 0:& \beta \to \frac{(1-r)(1+2r)}{3r}\\ \alpha =1:& \beta \to \frac{3(1-r)}{2+r}\\ \alpha \to \infty :& \beta \to \frac{(1-r)(5+r)}{3(1+r)}\end{array}.$$(21)

To consider rock with any value of porosity and, therefore, to consider realistic rocks we must evoque the “effective medium theory” for taking into account the interactions of neighboring pores as explained by Zimmerman [5]. Then, applying the effective medium theory proposed by Maxwell [28], the actual Fricke model [27] can be written as

$$\frac{{\lambda }_{\mathrm{eff}}\left(\phi \right)}{{\lambda }_{R }}=\frac{\left(1-\phi \right)\left(1-r\right)+\mathrm{r\beta }\phi }{\left(1-\phi \right)\left(1-r\right)+\beta \phi },$$(22)

where r and β are defined as previously. This model is tested in the next section.

The limit for the zero porosity verifies $\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$ and for φ = 1, the model tends to λ_f as required. Moreover, Fricke’s equation always satisfies the Hashin–Shtrikman bound.

3 Comparison with laboratory data

To test all these models, the published experimental data is used from samples of sandstone collected in five wells [3] in Perth Basin, Western Australia and Soultz-sous-Forets Basin, astern France using [2] to test thermal conductivity calculation from P-wave velocity. Two sets of data are available

• Measured thermal conductivity for dry sample, λ_f = λ_air,

• Measured thermal conductivity for saturated sample, λ_f = λ_water.

According to Clauser and Huenges [29] and using Esteban et al. [3], for water and air under room conditions (1 bar and 22 °C), the corresponding thermal conductivity is defined as

$${\lambda }_{\mathrm{air}}=0.025 \mathrm{W}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1},$$(23)

$${\lambda }_{\mathrm{water}}=0.6 \mathrm{W}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1}.$$(24)

These values are used together with the laboratory data and summarized in Table 1.

The purpose is to match different models to the experimental data. Independently on the fluid (air or water), the value of the thermal conductivity of the rock λ_R is assumed to be a unique constant within the same model. Thus, the both datasets, air and water, are considered together to search λ_R. Therefore, the matching methodology is constrained as follows.

First, an Objective Function (OF) based on least squares errors is built

$$\mathrm{OF}= \sum _{i=1}^N{\left(\frac{{\lambda }_i^{\mathrm{Air}}-{\lambda }_{\mathrm{eff}}\left({\phi }_i,{\lambda }_{\mathrm{air}},{\lambda }_{\mathrm{R}}\right)}{{\lambda }_i^{\mathrm{Air}}{\omega }_i}\right)}^2+\sum _{i=1}^N{\left(\frac{{\lambda }_i^{\mathrm{water}}-{\lambda }_{\mathrm{eff}}\left({\phi }_i,{\lambda }_{\mathrm{water}},{\lambda }_{\mathrm{R}}\right)}{{\lambda }_i^{\mathrm{water}}{\omega }_i}\right)}^2,$$(25)

where N is a number of experimental entries. Each term is weighted by ω_i which is assumed to be equal to one for all the experimental data. To impose the convergence at the two known points λ_Air and λ_water for φ = 1, their weight is set to ω = 0.1 (Tab. 1).

In order to get optimal model parameters to approach the data, the OF is minimized using the Powell’s dog leg minimization method.

Then, the standard deviation is used to estimate the veracity of the results

$$\mathrm{SDE}=\sqrt{\frac{1}{2(N-1)}{\sum }_{i=1}^N\left({\left({\lambda }_i^{\mathrm{Air}}-{\lambda }_{\mathrm{eff}}\left({\phi }_i,{\lambda }_{\mathrm{Air}},{\lambda }_{\mathrm{R}}\right)\right)}^2+{\left({\lambda }_i^{\mathrm{Water}}-{\lambda }_{\mathrm{eff}}\left({\phi }_i,{\lambda }_{\mathrm{Water}},{\lambda }_{\mathrm{R}}\right)\right)}^2\right)}.$$(26)

In order to bound a maximal thermal conductivity, the study of Orlander et al. [30] for several sandstone outcrops originate from Fontainebleau (France), Castlegate (USA), Bentheim (Germany), Obernkirchen (Germany) and Berea (USA), is used. A maximal thermal conductivity for sandstone obtained by Orlander et al. [30] is near 8 Wm⁻¹ K⁻¹ for φ → 0. However, to release the parameter restrictions, the bounds for λ_R are set to

$${\lambda }_{R }\in \left[{\lambda }_{\mathrm{Air} },16\right].$$(27)

3.1 Analysis of mixing laws

The mixing laws are given as simple means, which depend only on the thermal rock conductivity λ_R. As mentioned before, the measure of λ_R is very difficult, thus, for each model, λ_R is obtained in order to match as close as possible the experimental data.

In Figure 1, the resulting mixing laws are compared to the experimental data. The precision of each model is characterized by the corresponding OF and the standard deviation errors (SDEs) given in Table 2 together with resulting thermal rock conductivity.

The worst fit is obtained with the weighted harmonic mean. Moreover, the optimal value of λ_R corresponds to the set maximal bound which implies a bad convergence of the minimization.

Then, the Hashin and Shtrikman maximal and minimal bounds and mean models also yield high OF and SDE. The corresponding fitting curve for minimal bound matches with the data for the saturated sample, while for the non-saturated sample the minimal bound is found extensively low compared to the data. Moreover, the obtained optimal value of λ_R is higher than value proposed by Orlander et al. [30]. The maximal bound match with non-saturated sample data. The Hashin and Shtrikman mean model yields small SDE, however the objective function keeps high.

The weighted arithmetic mean has also a high error compared to the experimental data. This model is a linear correlation and it is not adapted for the used set of experimental results which has a non-linear trend. It can be noted that the SDE is very close to those of the Hashin and Shtrikman mean model.

Finally, the best fit is obtained with the weighted geometric mean model. It yields the smallest objective function and the standard deviation for the thermal conductivity of rock λ_R = 6.74 $\mathrm{W}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$. Moreover, this model yields the required bounds for porosity limits, $\underset{\phi \to 1}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{f}}$ and $\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$.

3.2 Analysis of empirical and theoretical models

The empirical and theoretical models are more complex as the mixing laws. Additionally to the thermal conductivity of rock, these models have other fitting parameters. Similarly to the mixing laws, these parameters can be determined through minimizing the objective function given before or defined by the pore structure. The minimization parameters should satisfy the following conditions.

• The Asaad model [23] tends to the weighted geometrical mean when c = 1. Moreover, it was demonstrated that the weighted geometrical mean yields the best fit with the experimental data. Therefore, the empirical exponent c of the Asaad model is assumed to be closed to 1.

• The parameter n of the Sugawara and Yoshizawa model [7, 24] is strictly positive and found in the range from 0.01 to 10.

• Zeb and Maqsood model [4] found the values of the empirical exponent z between 0.1 and 5.

• Zeb et al. [21] found the values of the model parameter m in the range [0.1, 5].

• The parameter ψ of the Veerendra and Chaudhary model is to account only for a high thermal conductivity ratio $\frac{{\lambda }_{\mathrm{R}}}{{\lambda }_{\mathrm{f}}}$. Moreover, in order to maintain the condition on the thermal conductivity at porosity limits, ${\underset{\phi \to 1}{\lim }\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{f}}$ and ${\underset{\phi \to 0}{\lim }\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_R$, it should be close to zero.

• For the Krupiczka model, we can think that the parameters values must be closed to those found by the author, a₁ = 0.28, a₂ = −0.757 and B = −0.057.

• The parameter M of Fricke model is a factor which depends on the aspect ratio of the pore. For spherical pore, M = 0. Consider that for the studied sample of sandstone the pores are nearly spherical. Thus, the parameter M should be found close to zero.

Table 3 summarizes the resulting optimal parameters as well as the corresponding objective function values and the SDEs. For all the models, the minimization parameters are obtained in agreement with the required conditions. The fits obtained for each empirical model are compared to the experimental data in Figure 2.

The worst results are obtained with the Zeb and Maqsood model and followed by Zeb, Gurmani, Ali and Maqsood and Veerendra and Chaudhary models. Similarly to the weighted harmonic mean, the Zeb, Gurmani, Ali and Maqsood model yields the thermal rock conductivity equal to the imposed upper bound, λ_R = 16, which means the poor convergence of the solution. Then, the Sugawara and Yoshizawa and Krupiczka models don’t yield a good approximation as well. Moreover, at zero porosity limit, the Krupiczka model has not a finite solution, therefore it is not of interest for the present study. Indeed, a not finite solution can bring a divergence of the reservoir simulator. For these five models, it seems impossible to match the data using common minimization parameters for saturated (water) and non-saturated (air) cases.

Finally, the Asaad and Fricke models give the best match with the experimental results. Applying the optimal value of the Asaad model parameter c = 0.99992, this model degenerates to the weighted geometric mean. The errors obtained with Fricke model are insignificantly higher than those of Asaad and geometric mean models. However, it is the only model which has a strong physical background compared to the others mostly pure mathematical models.

3.3 The empirical model with parameter dependent on saturated fluids

Consider four empirical models which are revealed unsatisfactory for the data approximation using unique parameters independent on the fluid, the Sugawara and Yoshizawa, Zeb and Maqsood, Zeb, Gurmani, Ali and Maqsood and Veerendra and Chaudhary models. The Krupiczka model is not considered here since it doesn’t yield the finite thermal conductivity values at zero porosity limit. Suppose that the fitting parameters are fluid dependent and λ_R keeps a unique constant for water and air curves. Thus, for the corresponding objective function there is one minimization parameter more than before. The results are summarized in Table 4.

As previously, the Zeb, Gurmani, Ali and Maqsood model is not correctly converged since the resulting value of λ_R is found on its upper limit, λ_R = 16. The Veerendra and Chaudhary models are found the same for fluid dependent parameters as for unique ones (Fig. 3).

Considering the fluid dependent parameters for Sugawara and Yoshizawa and Zeb and Maqsood, allow to obtain the better approximation of the data. The corresponding objective functions and the SDEs are as small as for the geometric mean (Tab. 4).

The dependency on type of fluid increases the number of the model parameters. Moreover, in a reservoir simulation, the algorithm must choose the good fitting parameter function of saturated fluids in each cell. It becomes especially problematic when the rock is saturated with many fluids (steam and water, for example). The fitting parameter is available only for one fluid. Therefore, when the cell is saturated of several fluids, the correlation is not available anymore. Thus, the choice of the fluid dependent parameters makes sense only if the rock is saturated of one fluid and don’t evaluate with time (change porosity, saturation, etc...).

4 Conclusion

The analysis of various models is summarized in Table 5.

Combining the model accuracy and required for the numerical simulator conditions, the Fricke and Asaad models reveal to be the most appropriate. Both of them yield a good precision (low OF and SDE), they have finite bounds for porosity limits (φ = 0, 1) and require only one additional to λ_R fitting parameter.

For both models, at zero porosity, the condition $\underset{\phi \to 0}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{R}}$ is always satisfied, and for porosity φ = 1, the required condition $\underset{\phi \to 1}{\lim }{\lambda }_{\mathrm{eff}}\left(\phi \right)={\lambda }_{\mathrm{f}}$ is also fulfilled.

It should be noted that the Asaad model is comparable to the geometric mean, but an additional fitting parameter makes the Asaad model more flexible which may be interesting if dealing with another rock type, for example.

Finally, it should be kept in the mind that the Fricke model has a strong physical background compared to other purely mathematical models. In view of the results, the Fricke model seems to be the most suitable: correct limits, only one parameter, very good matching and strong physical background.

References

© C. Preux & I. Malinouskaya, published by IFP Energies nouvelles, 2021

All Tables

All Figures