Well drilling in permafrost regions: dynamics of the thawed zone
Published: June 6, 2019
Latest article update: Aug. 6, 2023
In the cold regions, warm mud is usually used to drill deep wells. This mud causes formation thawing around wells, and as a rule is an uncertain parameter. For frozen soils, ice serves as a cementing material, so the strength of frozen soils is significantly reduced at the ice–water transition. If the thawing soil cannot withstand the load of overlying layers, consolidation will take place, and the corresponding settlement can cause significant surface shifts. Therefore, for longterm drilling or oil/gas production, the radius of thawing should be estimated to predict platform stability and the integrity of the well. It is known that physical properties of formations are drastically changed at the thawing–freezing transition. When interpreting geophysical logs, it is therefore important to know the radius of thawing and its dynamics during drilling and shutin periods. We have shown earlier that for a cylindrical system the position of the phase interface in the Stefan problem can be approximated through two functions: one function determines the position of the meltingtemperature isotherm in the problem without phase transitions, and the second function does not depend on time. For the drilling period, we will use this approach to estimate the radius of thawing. For the shutin period, we will utilize an empirical equation based on the results of numerical modelling.
Keywords
Radius of thawing, Stefan problem, freezeback period, permafrost temperature
The problem of phase change (Stefan problem) around a cylindrical source occurs in many engineering designs: oil/gas production/injection wells in permafrost areas, underground pipelines, steam production boreholes, melting of metals and storage of nuclear waste. In the cylindrical coordinates, an exact solution of the Stefan problem in the infinite domain exists only for two cases: subcooled liquid, which freezes while the solidified region remains at the fusion temperature; and a line source, which extracts energy at some constant rate per unit length (Carslaw & Jaeger 1959). Heat balance integral solutions were applied to the problem of finite superheat, when the initial temperature of the medium is lower than the fusion temperature (Tien & Churchill 1965; Sparrow et al. 1978; Lunardini 1988). Different numerical methods were used to solve the heat valance integral equations, but the results were essentially the same (Lunardini 1988). A coordinate transformation method reduced the problem with a variable phase change, such as a cylinder, to one with a constant phase change area (Lin 1971). An interesting effective thermal diffusivity concept was introduced (Churchill & Gupta 1977). It was assumed that the actual thermal diffusivity could be replaced by the effective thermal diffusivity, which includes the latent heat. The accuracy of these two methods is limited to certain ranges of dimensionless parameters (Lunardini 1988).
When wells are drilled through permafrost, the natural temperature field of the formations in the vicinity of the borehole is disturbed and the frozen rocks thaw for some distance from the borehole axis (e.g., Romanovsky et al. 2007; Eppelbaum et al. 2014). Making geothermal measurements to determine the static temperature of the formation and the permafrost thickness must be postponed for some period after completion of the drilling. This is the socalled restoration time. A lengthy restoration period of up to 10 years or more is required to determine the temperature and thickness of permafrost with sufficient accuracy (Lachenbruch & Brewer 1959; Melnikov et al. 1973; Shiu & Beggs 1980; Judge et al. 1981; Taylor et al. 1982; Clow 2014; Eppelbaum et al. 2014). The duration of the refreezing of the layer thawed during drilling is very dependent on the natural temperature of the formation; therefore, the rocks at the bottom of the permafrost refreeze very slowly (e.g., Dobinski 2011). The position of the interface of the thawing–freezing transition can be determined with resistivity and sonic logs. For example, the transition from higher resistivity and velocity readings to lower values can be considered as the position of the thawing radius. In our recent paper, we suggested a method to estimate how long it takes for formations thawed by drilling to refreeze (Kutasov & Eppelbaum 2017a, b). The method requires just three temperature logs taken after the freezeback is completed. Earlier we conducted numerical modelling and found that the dimensionless time of refreezing can be expressed as a function of two dimensionless parameters: the dimensionless radius of thawing and the dimensionless latent heat density (Kutasov 1999, 2006). Kutasov et al. (1977) proposed an effective approach for solving the problem of phase change (Stefan problem) around a cylindrical source of heat. It was shown that a known solution for a planar system can be utilized to obtain an approximate solution of the Stefan problem for a cylindrical source of heat with a constant temperature. Later (Kutasov 1998), we used an adjusted heating time concept (Kutasov 1987, 1999) to determine the position of a thawing temperature isotherm in the problem without the ice–water transition. The results of numerical solutions presented in Taylor (1978) were used to verify the results of our calculations.
It was shown that for a cylindrical system the position of the phase interface in the Stefan problem can be approximated through two functions. One function (r_{m}) determines the position of the meltingtemperature isotherm in the problem without phase transitions, and the second function Ψ does not depend on the time (Kutasov et al. 1977; Kutasov 1998, 2006; Eppelbaum et al. 2014). Recently, Wang et al. (2017) reported an interesting study in which they developed a coupled thermal model of the wellbore–permafrost system. It considered the latent heat of fusion, water migration and the change in thermal parameters. In this study, work done by Ramey (1962) and the temperature prediction model for incompressible fluids developed by Wu & Pruess (1990) and Hasan and Kabir (1994) were used to describe heat transfer in the boreholeformation system. Numerical solutions of moisture transfer problems for frozen soil were based on the equations given by Harlan (1973). For numerical simulations, Wang et al. (2017) selected a simulation well: well depth is 3100 m and thickness of permafrost is 750 m. The suggested coupled thermal model of the wellbore in permafrost regions allowed for the estimation of the volume of thawed permafrost (and the corresponding radius of thawing) and radial temperature of formations. It was also shown that the use of low thermal conductivity cements significantly reduces the volume of thawed permafrost. The most important result of the study by Wang et al. (2017) is the following: very low rate of water migration (only 4.72 × 10^{−5} m^{3}/s) (see the parameters for a simulated well, Table 1 in Wang et al. 2017). As a result, the heat transferred with the water migration is negligible compared with the heat caused by heat conduction. The simulation results (radius of thawing) can be validated against data obtained by geophysical logging (see also end of the “Time of the complete freezeback” section) and by temperature logs taken during shutin (the degree of thermal disturbance caused by drilling).
Table 1 Comparison of values of position of 0°C isotherm: r_{Dm} (Eqn. 6), R (numerical modelling; Taylor 1978).  
r_{Dm}  t_{D} = 100  t_{D} = 300  t_{D} = 500  t_{D} = 1000 
 
q  R  q  R  q  R  q  R 
 
2  0.761  1.968  0.797  1.977  0.810  1.983  0.826  1.984 

3  0.621  2.934  0.678  2.950  0.700  2.949  0.725  2.955 

4  0.523  3.893  0.594  3.918  0.621  3.924  0.653  3.926 

5  0.447  4.864  0.529  4.887  0.560  4.897  0.597  4.899 

6  0.386  5.836  0.476  5.857  0.511  5.854  0.551  5.878 

7  0.335  6.820  0.432  6.816  0.469  6.826  0.513  6.835 

8  0.292  7.805  0.393  7.806  0.433  7.794  0.479  7.826 

9  0.255  8.797  0.360  8.766  0.401  8.774  0.450  8.787 

10  0.223  9.793  0.330  9.753  0.373  9.741  0.424  9.753 

11  0.195  10.797  0.303  10.749  0.347  10.741  0.400  10.742 

12  0.170  11.826  0.279  11.735  0.324  11.721  0.379  11.694 

13  –  –  0.257  12.735  0.303  12.703  0.359  12.684 

14  –  –  0.238  13.683  0.284  13.673  0.341  13.652 

15  –  –  0.219  14.722  0.266  14.672  0.324  14.640 

16  –  –  0.202  15.739  0.249  15.695  0.308  15.64 

The objective of this study is to estimate the dynamics of the unfrozen zone (radius of thawing during drilling and shutin periods). To demonstrate the applicability of the suggested equations in estimating the formation temperature, the radius of thawing and the time of complete freezeback in a field case are presented.
The results of field and analytical investigations have shown that in many cases the effective temperature (T_{w}) of the circulating fluid (mud) at a given depth can be assumed constant during drilling or production (Lachenbruch & Brewer 1959; Jaeger 1961; Edwardson et al. 1962; Ramey 1962; Kutasov et al. 1966; Raymond 1969). Here we should note that even for a continuous mud circulation process the wellbore temperature is dependent on the current well depth and other factors. The term “effective fluid temperature” is used to describe the temperature disturbance of formations while drilling. We should note that the effective temperature T_{w} takes into account the changes in heat transfer during all period of mud circulation on the given depth. The good agreement between calculated values of transient (shutin) temperatures and the results of temperature surveys in numerous wellbores in permafrost regions confirms this assumption (Kutasov 1999; Kutasov & Eppelbaum 2003). Because of its low permeability, water migration is usually neglected in permafrost studies. Recent numerical modelling was used to confirm this assumption (Wang et al. 2017).
Lachenbruch & Brewer (1959) have shown that the wellbore shutin temperature mainly depends on the amount of thermal energy transferred to (or from) formations. Therefore, for every depth a value of T_{w} can be estimated from shutin temperature logs. Drawing from Kutasov’s earlier work (1977, 1998, 2006), we can write that
where r_{w} is the radius of the well and h_{0} is the radius of thawing. The function r_{m} determines the position of the meltingtemperature isotherm in the problem with phase transition. The function Ψ shows to what extent the melting process affects the position of the meltingtemperature isotherm and does not depend on time. From physical considerations, it follows that 0 < Ψ < 1. The function r_{m} is known, but is expressed through a complex integral. By introducing an adjusted dimensionless heating time (t_{D}*), we have found (Kutasov 1987) that the exponential integral (a tabulated function) can be used to approximate the function T_{D} = T(r_{D}, t_{D}):
where r is the cylindrical coordinate, r_{D} is the dimensionless radial distance, Ei(x) is the exponential integral, T_{w} is the temperature at the wall of the borehole, T_{m} is the temperature of melting, T_{f} is the initial temperature, T(r, t) is the radial temperature, T(r_{D}, t_{D}) is the dimensionless radial temperature, t_{D} is the dimensionless time and a_{1} is the thermal diffusivity of the frozen formation.
The correlation coefficient G(t_{D}) varies in the narrow limits: G(0) = 2 and G(∞) = 1. From Eqn. 2 (assuming that T = T_{m}) we can determine the position of the meltingtemperature isotherm. In our case, the melting temperature (for pure ice) is 0°C, and the following equation can be used to determine r_{m}
As can be seen from Table 1, Eqn. 6 approximates the position of 0°C isotherm with good accuracy.
Analysis of physical characteristics indicates that at the very small times the solutions for cylindrical systems approach those for plane systems. For this reason, we assumed that the function Ψ can be determined from the known solution for the plane Stefan problem (Carslaw & Jaeger 1959). In a plane system, the position of the solid–melted interface is
Here, subscripts _{1} and _{2} correspond to the frozen and thawed zones, respectively, l is the thermal conductivity, c is the specific heat, T_{s} is the surface temperature, r is the density, L is the latent heat per unit of mass, w is the ice mass content per unit of volume, Φ is the probability integral and Φ* = 1 – Φ.
Introducing the dimensionless parameters, we obtain
As we can see (Eqn. 7), the depth of thawing X is a product of two functions and the function l does not depend on time. The changes of thermal properties are due to thawing, and therefore, by using our assumption, we obtain from Eqn. 8 that
The results of a numerical solution (Taylor 1978) and results of calculations (after Eqns. 1–10) were compared (Kutasov 2006). It was found that both approximate solutions are in satisfactory agreement.
After the cessation of the drilling process, the radius of thawing and the radius of thermal influence will increase for a definite period of time Δt_{0} at the expense of the heat accumulated in the thawed zone. Correspondingly, the increase of the radius of thawing will be Δh_{0}. We suggested an empirical relationship to estimate the parameter Δt_{0} (Kutasov 1999: 181). Hydrodynamical modelling (Kutasov 1999) has shown that the maximum value of the thawing radius can be given by
where t_{c} is the drilling mud circulation time at a given depth.
As an example, the results of one iteration of the hydrodynamic modelling are presented in Fig. 1. The input parameters are the time of heating disturbance t_{c} = 2200 hr, the temperature of drilling mud is T_{m} = 8°C, the temperature of the formation is T_{f} = −2°C, the well radius is r_{w} = 0.1 m, thermal diffusivity and thermal conductivity of the thawed formation is a_{t} = 0.0030 m^{2}/h and l_{t} = 2.0 kcal/(h·m·°C) and the thermal conductivity of the permafrost (frozen area) is l_{f} = 2.5 kcal/(h·m·°C).
We obtained that h(t_{c}) = 1.12 m, Δt_{0} = 113 hours and
Below we will neglect the difference between h_{c} and h_{max}. Refreezing of the thawed zone starts at the moment of time t = t_{0} and ends at t = t_{ep} (Fig. 2).
We should note that only a part of the formation’s pore water changes to ice at 0°C. With further lowering of the temperature, phase transition of the water continues, but at steadily decreasing rates. The amount of unfrozen water is practically independent of the total moisture content for a given soil (Tsytovich 1975).
Time of the complete freezeback
It was assumed that the heat flow from the thawed zone to the thawed zone–frozen zone interface can be neglected. The results of hydrodynamical modelling have shown that this is a valid assumption (Kutasov 1999). In this case, the Stefan equation–energy conservation condition at phase change interface (r = h) is
Assuming the semisteady temperature distribution in frozen zone (a conventional assumption), we obtain
where r_{if} is the radius of thermal influence during the freezeback period. The ratio D_{f} = r_{if}/h was determined from a numerical solution. A computer programme was used to obtain a numerical solution of a system of differential equations of heat conductivity for frozen and thawed zones and the Stefan equation (Kutasov 1999, 2006). 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 porosity, r_{w} is the water density, r_{f} is the formation density, r_{w} is the well radius, and l_{f} and a_{f} are the thermal conductivity and diffusivity of frozen formations, respectively.
From Eqns. 14–16 and the condition h(t_{0}) = h_{max}, we obtained (Kutasov 1999, 2006)
where H is the dimensionless radius of the thawing during shutin.
From Eqn. 18 at h = r_{w}, we obtain the relationship for the duration of complete freezeback t_{cf} = t_{ep} – t_{0} (Fig. 2).
where t_{epD} is the dimensionless time of refreezing.
It is known that the electric resistivity of frozen sediments is affected by the freezing–thawing transition to a greater extent than the seismic velocities. Seismic velocities may increase by 2–10 times in transition to a frozen state, whereas the electrical resistivity may increase by 30 to 300 times in the same temperature interval (Hnatiuk & Randall 1977; Dobinski 2011). The dynamics of the thawed zone (the radius of thawing, Eqn. 18) while refreezing can, therefore, be monitored by geophysical methods (electric resistivity and seismic logs).
The input data and location of this well are presented in Table 2. This is a unique wellbore. Indeed, the first temperature log was taken after five days of shutin and the last after three years of shutin (Table 3). As we can see, the refreezing time increases with sufficiently high accuracy with the formation temperature increasing (Table 4). Figure 3 shows the results of calculations after Eqn. 1. It is clear that Eqn. 1 approximates the observed shutin temperatures with sufficiently high accuracy.
Table 2 Input data and location of the example well (USGS 1998).  
Site code  PBF 
Site name  Put River, N1, Alaska 
Latitude  70° 19′07′′N 
Longitude  148°54′35′′W 
Surface elevation (m)  8 
Casing diameter (cm)  51 
Hole depth (m)  763 
Date of drill start  02/09/70 
Drilling time (days)  44 
Number of logs  9 
Shutin time (days)  5–1071 
Table 3 Observed shutin temperatures (in °C) in well Put River N1, Alaska.  
z (m)  Shutin time (days) 
 
5  22  34  48  66  91  117  163  1071 
 
30.48  −0.400  −2.686  −4.793  −6.252  −7.040  −7.602  −7.970  −8.716  −9.167 

45.72  −0.300  −2.093  −4.507  −6.012  −6.910  −7.511  −7.900  −8.428  −9.052 

60.96  −0.250  −2.941  −4.911  −6.148  −6.950  −7.497  −7.860  −8.263  −8.957 

91.44  −0.300  −1.633  −4.101  −5.646  −6.590  −7.227  −7.620  −7.965  −8.771 

121.92  −0.210  −0.882  −2.565  −4.781  −6.060  —  −7.250  −7.624  −8.520 

152.40  −0.030  −0.976  −1.852  −3.173  −4.760  −5.880  −6.510  −7.026  −8.124 

182.88  0.020  −0.757  −1.217  −2.506  —  —  −6.140  −6.600  −7.619 

213.36  0.200  −0.490  −0.805  −1.528  —  −5.049  −5.680  −6.111  −7.144 

243.84  0.380  −0.433  −0.608  −0.950  −2.660  —  —  −5.521  −6.602 

274.32  0.640  −0.418  −0.555  −0.823  —  −3.186  —  −4.822  −6.029 

304.80  0.740  −0.379  −0.506  −0.682  −1.150  —  −3.610  −4.325  −5.462 

335.28  0.910  −0.325  −0.451  −0.577  —  −1.720  −2.840  −3.705  −4.935 

365.76  1.040  −0.322  −0.452  −0.579  −0.800  —  −2.290  −3.258  −4.454 

396.24  1.230  −0.354  −0.505  −0.644  −0.860  −1.248  −1.880  −2.677  −4.039 

426.72  1.220  −0.280  −0.415  −0.517  −0.630  —  −1.130  −1.796  −3.453 

Table 4 Results of calculations B and T_{f}.  
z (m)  t_{0}  t_{2}  t_{3}  T_{2} (°C)  T_{3} (°C)  T_{f} (°C)  B (°C) 
30.48  22  48  66  −6.252  −7.040  −9.101  2.289 
45.72  22  48  66  −6.012  −6.910  −9.248  2.621 
60.96  22  48  66  −6.148  −6.950  −9.029  2.352 
91.44  22  48  66  −5.646  −6.590  −9.016  2.797 
121.92  22  48  66  −4.781  −6.060  −9.316  3.830 
152.40  22  48  66  −3.173  −4.760  −8.763  4.804 
182.88  48  117  163  −6.140  −6.600  −7.609  1.885 
213.36  48  91  163  −5.049  −6.111  −7.182  2.034 
243.84  66  163  1071  −5.521  −6.602  −6.767  1.812 
274.32  91  163  1071  −4.822  −6.029  −6.190  1.401 
304.80  117  163  1071  −4.325  −5.462  −5.587  0.892 
335.28  117  163  1071  −3.705  −4.935  −5.069  0.970 
365.76  117  163  1071  −3.258  −4.454  −4.584  0.949 
396.24  117  163  1071  −2.677  −4.039  −4.186  1.088 
426.72  117  163  1071  −1.796  −3.453  −3.631  1.332 
The medium is a permafrost formation (sandstone). Only heat transfer by radial conduction is considered. The following parameters are introduced: the radius of the wellbore is 0.255, the thermal conductivity of frozen and unfrozen formations are l_{f} = 4.40 and l_{un} = 3.84 Wm^{−1} K^{−1}, specific heat c_{f}= 950 and c_{un} =1138 Jkg^{−1} K^{−1}, the density of sandstone is r_{f} = 2483 kg m^{−3}, the density of water/ice is r_{w} = 1000 kg m^{−3}, porosity is f = 0.09, and latent heat L = 334 960 J kg^{−1} for water/ice. The latent heat density of the medium is c = Lr_{w}f = 334 960×1000×0.09 = 30×10^{6} Jm^{−3}. The duration of the source disturbance is t_{c}. The phase change is assumed to be at 0°C. The temperature of drilling mud is assumed to be 8°C.
Recently, we suggested a new approach in predicting the undisturbed formations temperatures from shutin temperature logs in deep wells (Kutasov & Eppelbaum 2018). The main features of the suggested method are the following: the refreezing of the thawed formations (around the wellbore) is completed; the temperature logs are taken after refreezing and the starting point in the well thermal recovery is moved from the end of well completion to the moment of time when the first shutin temperature log was conducted. It is shown that after refreezing the further cooling of a well can be approximated by a constant (per unit of length) linear heat source. Hence, a modified Horner equation can be used for predicting the temperature of frozen formations for estimation of the formation temperature. A simple method to process field temperature data is presented. To demonstrate this approach, temperature shutin time data for four depths from four wells in Alaska were successively used (Kutasov & Eppelbaum 2018).
We did not have an access to the drilling journal of this borehole to estimate the time of thermal disturbance, t_{c} (at a given depth), caused mainly by circulation of the drilling mud. For this reason, we consider thermal disturbance to start at the moment of time when the bit reached a given depth and the end of thermal disturbance is the moment when the drilling operations are completed. Then
where t_{tot} is the total drilling time, H_{t} is the vertical well depth and z is the current depth. Now we can consider that the period of time t*_{c} = t_{c} + t_{0} as a new “thermal disturbance” period (t_{0} = t_{s}_{1}). Here t_{s}_{1} is the shutin time for the first temperature log (Fig. 4). Therefore, for logs 2 and 3 the modified Horner equations are
From Eqns. 21 and 22, we can estimate the parameter B for a given depth and the formation temperature:
The results of calculations after Eqns. 22 and 23 are presented in Table 4.
The radius of thawing in the drilling period
Results of calculations after Eqns. 1, 2, 6–9 are presented in Table 5.
Table 5 The dimensionless radius of thawing, I_{a} = 0.7285, I_{ka} = 0.9780.  
z (m)  t_{D}  θ  I_{f}  Ψ  H 
30.48  80.92  1.1380  1.327  0.8382  3.19 
45.72  79.24  1.1560  1.327  0.8417  3.16 
60.96  77.55  1.1290  1.327  0.8365  3.18 
91.44  74.18  1.1270  1.327  0.8361  3.15 
121.92  70.82  1.1640  1.327  0.8433  3.07 
152.40  67.45  1.0950  1.327  0.8297  3.13 
182.88  64.08  0.9510  1.327  0.7985  3.30 
213.36  60.72  0.8980  1.327  0.7860  3.35 
243.84  57.35  0.8460  1.327  0.7731  3.39 
274.32  53.98  0.7740  1.327  0.7541  3.47 
304.80  50.62  0.6980  1.327  0.7325  3.56 
335.28  47.25  0.6340  1.327  0.7128  3.64 
365.76  43.88  0.5730  1.327  0.6927  3.70 
396.24  40.51  0.5230  1.327  0.6751  3.75 
426.72  37.15  0.4540  1.327  0.6488  3.84 
where t_{D} is the dimensionless drilling mud circulation time.
It is interesting to note that the dimensionless radius of thawing varies in narrow limits (Table 5). This can be explained by a combination of two factors: (1) the radius of thawing increases with depth (due to temperature increase with depth); (2) the time of mud circulation reduces with depth and the radius of thawing reduces with depth. The time of complete freezeback was calculated after Eqn. 19 (Table 6).
Table 6 The time of the complete freezeback.  
z (m)  H  T_{f} (°C)  t_{cf} (days) 
30.48  3.190  −9.101  7.46 
45.72  3.160  −9.248  7.19 
60.96  3.180  −9.029  7.48 
91.44  3.150  −9.016  7.33 
121.92  3.070  −9.316  6.68 
152.40  3.130  −8.763  7.45 
182.88  3.300  −7.609  9.74 
213.36  3.350  −7.182  10.71 
243.84  3.390  −6.767  11.72 
274.32  3.470  −6.190  13.57 
304.80  3.560  −5.587  16.01 
335.28  3.640  −5.069  18.66 
365.76  3.700  −4.584  21.54 
396.24  3.750  −4.186  24.46 
426.72  3.840  −3.631  30.02 
Drilling intermediate and deep wells in permafrost areas usually includes a warm mud application with unknown dynamics of the formation thawing around the wells. A new method based on the phase change (Stefan problem) around a cylindrical source is proposed. This method allows for the estimation of the radius of thawing during drilling and shutin periods. Determining formation temperature and estimating the time of complete freezeback are illustrated with an example of field case.
The authors would like to thank two anonymous reviewers, who thoroughly reviewed the manuscript. Their critical comments and valuable suggestions were helpful in preparing this article.