The Possibility to Estimate In Situ the Porosity of Permafrost: A Preliminary Study

It is well-known that only a part of the formation’s pore water turns into ice at 0 ^oC. With the additional lowering of the temperature, the phase transition of the water continues, however, at gradually decreasing rates. The quantity of unfrozen water is practically independent of the total moisture content for soil with the concrete physical-chemical parameters (Tsytovich, 1975). Frozen soil is the matter in which stresses and strains arise under the influence of an external load. These forces are not constant but vary with time. It gives rise to the relaxation of stresses and creeps (i.e., increased strains over time). These complex physical-chemical processes are called rheological ones. The vigorous development of the rheological processes in frozen soils is caused by the peculiarities of their internal relationships in which ice plays a significant role. Numerous laboratory experiments show that frozen formations’ mechanical, thermal, and rheological properties strongly depend on the moisture content. Hence, porosity and temperature are essential parameters of permafrost. It is known that the electric resistivities of frozen sediments are affected by the freezing-thawing transition to a greater extent than the seismic velocities. In the same interval of temperatures, seismic velocities may increase by 2 to 10 times in transition to a frozen condition, whereas the electrical resistivity may increase by 3·10 – 3·102 times (Hnatiuk and Randall, 1977; Dobinski, 2011; Eppelbaum, 2021). The position of the interface of the thawing-freezing transition can be determined with the application of the surface electric resistivity method and sonic logs in wells. For instance, the transition from higher resistivity and velocity values to lower ones can be considered the thawing radius's position. Thus, the position of the radius of thawing (because of well drilling or production) can be estimated by using different geophysical exploration methods. A method of estimating the refreezing time surrounding the wellbore thawed (during drilling) formations was suggested by Kutasov and Eppelbaum (2017, 2018). Only three temperature logs taken after the freezeback are completed and are needed to apply this method. The conducted numerical modeling indicates that the dimensionless time of refreezing can be expressed as a function of two dimensionless parameters: radius of thawing and latent heat density (Kutasov, 2006; Eppelbaum and Kutasov, 2018). From the last parameter, the porosity of formations can be estimated.

The only objective of this study is to offer a possible way to estimate permafrost's porosity in situ. This study is the first effort to determine this parameter in natural conditions. This parameter in situ estimation is of special importance when the permafrost parameters are constantly changing under the influence of global warming. To demonstrate the possibility of the new approach applicability of this new method, five cases of numerical modeling are presented.

Since deep wells drilling in permafrost areas usually use warm mud, a certain unknown degree of the formation thawing around the wells exists. The warm mud disturbs the borehole's temperature field, and as a result, the permafrost thaws. To calculate the static formation temperature and permafrost thickness, before conducting temperature logs, engineers must wait for some period after the entire completion of drilling. The duration of the refreezing of the layer thawed during drilling dramatically depends on the natural temperature of geological formation(s). Therefore, the rocks at the bottom of the permafrost zone freeze very slowly. A lengthy restoration period (up to ten years or more) is required to calculate the permafrost's temperature and thickness with the necessary accuracy.

As it was mentioned above, only a part of the formation's pore water changes to ice at 0^oC. The phase transition temperature interval exists from numerous laboratory and field experiments. With the subsequent lowering of the temperature, the water phase transition continues (Figure 1). The temperature interval of phase transition mainly depends on the mineralogical composition of the geological formations. Let us assume that the water-ice phase transition is completed at the time t_ep, and at least three temperature logs were taken after refreezing of thawed formations (Figure 1).

Kutasov and Eppelbaum (2003, 2017, 2018) have shown that the cooling process at t > tep (Figure 1) is like that of temperature recovery in borehole sections below the permafrost base (i.e., unfrozen formations). Let us assume that thermal recovery's starting point is t = tep. Thus, the thermal disturbance time is td + tep, where td is the time of drilling mud circulation at a given depth. It will be explained below that the subsequent borehole cooling can be approximated by a constant linear heat source (per unit of length) after refreezing.

Hence, a modified Horner equation (in some publications, this method is called KEM (for instance, Bassam et al., 2010; Liu et al., 2016)) can be used to predict frozen formations' temperature to estimate the formation temperature (Kutasov and Eppelbaum, 2017, 2018). Then

Then the values of shut-in temperatures can be determined (Figure 1):

Combining Eqs. (2) – (4), Eq. (5) is obtained for estimating the time of refreezing (t_ep).

For solving Eq. (5), Newton’s method was used (Grossman, 1977). In this method, a solution to an equation is obtained by defining a sequence of numbers that become successively nearer and nearer to the expected solution (Kutasov and Eppelbaum, 2017, 2018).

The parameter B is found from the following equation:

Furthermore, the temperature of formations can be obtained from Eqs. (2) – (4).

It was assumed that the heat flow influence from the thawed zone to the thawed zone – frozen zone interface could be neglected. The results of hydrodynamical modeling have established that this is an acceptable assumption (Eppelbaum et al., 2014). In this case, the Stefan equation – energy conservation condition at the phase change interface (r = h) can be applied

Assuming the semi-steady temperature distribution in the frozen zone (a conventional assumption): where r_if is the radius of thermal influence during the freezeback period. The ratio D_f = r_if/h was obtained from a numerical solution. A computer program (Kutasov, 1999) was used to find a numerical solution for a system of differential heat conductivity equations for frozen and thawed zones and the Stefan equation. It was found that where I_f is the dimensionless latent heat density, L is the latent heat per unit of mass, c_f is the specific heat of the formation, φ is the porosity, ρ_w is the water density, ρ_f is the formation density, r_o is the well radius, h is the radius of thawing, λ_f and a_f are the thermal conductivity and diffusivity of the frozen formations, respectively.

From Eqs. (7) – (9) and the condition H(t_ep) = 1, it was obtained (Kutasov, 2006) where t_epD is the dimensionless time of refreezing.

Taylor (1978) introduced a cylindrically symmetric source of thermal disturbance in a semi-infinite medium for the analysis of thermal borehole measurements. A numerical model has been developed to simulate the rising transient thermal regime in the mentioned model. The assumed medium is a permafrost sandstone formation. The model allows for a change of phase (ice-water, water-ice). In this model, only heat transfer by radial conduction is considered. The following parameters are introduced: the radius of cylindrical source r_o = 0.17m, the radial variable r, the thermal conductivity of frozen formation λ_f = 4.40 3.84 Wm^-1K^-1, the thermal conductivity of unfrozen formation λ_un = 3.84 Wm^-1K^-1 , specific heat c_f = 950, c_un =1138 Jkg^-1K^-1, the density of sandstone ρ_f = 2483 kg m^-3 , the density of water/ice ρ_w = 1000 kg m^-3, the porosity φ = 0.09, the density of water/ice ρ_w = 1000 kg m^-3, and latent heat L = 334960 J kg^-1 for the water/ice boundary. The latent heat density of the medium is χ = Lρ_wφ = 334960⋅1000⋅0.09 = 30⋅10⁶ Jm^-3. The duration of source disturbance is t_d. The temperature of 0 ^oC is assumed as a phase change. The calculated data in Table 1, Table 2, Table 3 and Table 4 are given for initial temperatures (T_f): -0.01^o, -5^o, -10^o C for source temperatures (T_w): 10 ^oC and 30 ^oC. In the construction of Table 1, Table 2, Table 3 and Table 4, dimensionless distance and dimensionless time were used:

Thus, t_d is the time of thermal disturbance, and t_s is the “shut-in” time. The radial temperature distributions during the thermal disturbance (T_d) and “shut-in” (T_s) periods are presented as follows:

In Table 1 and Table 2, the values of T_s = T_s (R = 1, t_sD) are used.

At this stage, to use Eqs. (1) – (6) the following parameters will be relaced by its dimensionless values

The results of calculations after Eqs. (1) – (6) (in the dimensionless units) are presented in Table 3 and Table 4 (columns 4-6).

As can be seen from Table 3 for case 1, the calculated dimensionless refreezing time varies in the restricted limits (18.5-21.3). At the same time, the results of numerical modeling provide corresponding values of 20-30 (Table 1, values in bold). For cases 2-5, the calculated refreezing times (Table 3 and Table 4, column 4) agree with the results of mathematical modeling (Table 1 and Table 2, values in bold).

The dimensionless thawing radius was defined as the position of the 0 ^oC isotherm and was found from T_d = f (R, t_dD) as 0 ^oC = f (R = H, t_dD). Here linear interpolation was used.

From Eq. (11) follows that the dimensionless latent heat density (I_f) can be determined as a function of dimensionless refreezing time and thawing radius. As at the solution of implicit Eq. (6), Newton’s method was again used to estimate the value of I_f from Eq. (11). The calculation results are presented in Table 3 and Table 4 (column 7).

And finally, the porosity is found form Eq. (10)

The calculation results are displayed in Table 3 and Table 4 (column 8).

It also should note that for cases 1, 2, 4, and 5, the calculated values of formation temperature (Table 3 and Table 4, column 6) are in excellent agreement with the assumed values (T_f = -5 ^oC and -10 ^oC) at numerical modeling. The difference T_f – T_fcal for case 3 can be explained by the low accuracy in determining the radius of thawing. As mentioned earlier, the linear interpolation method was used to estimate the thawing radius. Figure 2 also shows that the basic formula 1 (in dimensionless units) can be used to estimate the transient shut-in temperature.

The estimated porosity (ϕ) values are presented in Table 3 and Table 4 (column 8). For case 1, the ϕ values vary in narrow limits 0.086-0.095. The ranges for cases 2-5 are: 0.088-0.109; 0.081-0.011; 0.096-0.102; and 0.084-0.091, correspondingly. Thus, the average porosity values in all five cases are very close to the assumed numerical modeling value of ϕ = 0.09. The results of calculations shown in Table 3 and Table 4 (columns 4-8) testify that the basic formulas 1, 5, and 11 approximate the results of the numerical modeling with the necessary accuracy.

It should note that all input parameters (dimensionless heating and “shut-in” time, latent heat density, thermal properties of the formation a.o.) were specified in numerical modeling. Only the dimensionless thawing radius was defined as the position of the 0 ^oC isotherm and was found from tables T_d = f (R, t_dD) as 0 ^oC = f (R = H, t_dD). From the implicit Eq. (11) follows that the dimensionless latent heat density is a function of two parameters: dimensionless thawing radius (H) and dimensionless refreezing time (t_epD). Thus, to estimate the porosity of permafrost in field conditions (in-situ), it is essential to determine the values H and t_epD by temperature logging and other geophysical methods. For this reason, this work can be considered as a preliminary study. As it was mentioned earlier, the position of the interface of the thawing-freezing transition (radius of thawing) can be verified with geophysical methods – sonic logs and electric resistivity. Kutasov and Eppelbaum (2017, 2018) presented three examples of estimation of the refreezing time using temperature logging results. Furthermore, to validate the approach presented in this paper, the calculated porosity values should be compared with those obtained from cuttings and samples.