Development and testing of a mathematical model of the permafrost thawing processes during drilling of wells

A mathematical model conjugate heat exchange of the well consideration phase transitions and drilling fluid circulation during to drilling was developed. The method was tested on the experimental data. Good agreement with the experiment was shown. The calculation of the permafrost thaw around the well during drilling was carried out. Dependences of the thaw radius on the flow rate and properties of the drilling mud were obtained.


Introduction
Conditions of drilling oil and gas wells in the northern latitudes are determined by the permafrost presence. Permafrost are rocks that are constantly in conditions of negative temperatures. A numerous problems arise in the development of oil and gas condensate fields in the permafrost zone. Construction and operation of oil and gas facilities in permafrost is complicated by the problem of partial or complete thawing of soils near the objects of heat. There is an active thawing of surrounding rocks, resulting in subsidence, collapses, voids, which leads to a number of negative factors and even accidents in the drilling and operation of wells. This leads to long-term repair, well downtime and significant losses of drilling fluids.
Computational study of the rocks thaw around the wells in the drilling and in the exploitation is a subject of many papers [1][2][3]. However, despite the variety of investigate on the study of the thaw well processes at operation stage, there is no investigate on the study of the effect of drilling mud circulation at the construction stage. Although there is the formation of cavities, which leads to a violation of the well casing and bad cementing in the process of construction of wells. Therefore, the purpose of this study was to create and test a mathematical model that considered the processes of thaw rock at the well drilling during the circulation of the drilling fluid.

Mathematical model of drilling mud and fluids circulation in a well.
In the general case, a viscous fluid flow is described by a system of Navier-Stokes equations consisting of the mass conservation equation: , and the equations of motion or momentum conservation law: where v -is the fluid velocity vector, τ -is the tensor of viscous stresses, F -is the body forces vector, p -is the static pressure,  -is the fluid density. Since in most cases the drilling mud is non- Newtonian fluid, for simulating non-Newtonian flows we used well-known approach [4][5], in which the medium is considered as a nonlinear viscous fluid with the introduction of effective fluid viscosity    , which in general is dependent on shear rate. At that, the tensor of viscous stresses τ is determined as follows: The components of the strain velocity tensor D are of the form: shear rate  is the second invariant of the strain velocity tensor:

Mathematical model of phase transitions in the process of thaw / frost penetration of permafrost.
When thaw/frost penetration of permafrost is the process of melting or crystallization respectively. In this paper, the enthalpy -porosity formulation is used to model these processes. In this method the melt boundary is not explicitly tracked. Instead a value called the volume fraction of the fluid fraction is introduced, which indicates the volume fraction of the cells in the fluid state. The fluid fraction is calculated at each iteration based on the enthalpy balance. In the two-phase zone the fluid fraction has a value from 0 to 1. It is modeled as a "pseudo" porous medium where the porosity decreases from 1 to 0, which indicates the material solidification. When the material is completely solidified in the cell, the porosity becomes zero and therefore the velocity also drops to zero. The two -phase fluid-solid zone is considered as a porous zone with porosity equal to the volume fraction of the fluid fraction. The corresponding sources of energy and momentum flow (inflow) are added to the energy conservation and momentum conservation.
For solidification/melting problems the energy equation is written as: where H is the enthalpy,  is the density; v is the fluid velocity. Enthalpy is defined as: The latent heat content can be written in terms of the heat of the phase transition, L: The enthalpy-porosity approach simulates a two-phase region (partially hardened region) as a porous medium. The porosity in each cell is set equal to the fraction of the fluid phase in that cell. In fully solidified regions the porosity is equal to zero. The pulse drain due to the decrease in porosity in the two-phase zone has the following form: where  -volume of the fluid phase,  -numerical parameter (0.001) to prevent division by zero; mush A -constant of the two-phase region.

Boundary condition
The annual course of the average monthly temperature of the soil surface can be mathematically described, as well as the annual course of the air temperature, by a harmonic function: At the lower boundary of the computational area is set the value of the heat flow considered the temperature gradient:

The mathematical model test
The computational model was tested for several model tests. Well hydraulics were tested on experimental data from Escudier [8]. A calculation grid consisting of 20×120×3 (20 nodes in radius and 120 nodes in circumference and 3 in length of the channel) nodes were used to solve this problem. Figure 1 shows the distribution of the axial velocity contours and comparison with the experiment of the profile of the dimensionless axial component of the velocity in the cross section of the channel. The size was determined by the mean velocity U = 0.0728 m·s -1 . The distance between the cylinders, which is dimensioned by the width of the annular gap, is postponed along the x-axis. Evident the calculation is very well consistent with the experiment.  To test the model of permafrost thaw was used experiment from [9]. The model was tested on the annual temperature measurements of a drilling well in Alaska.
The two-dimensional problem was considered in our study. The width of the calculation area was 10 m, depth was 100 m. The calculation was carried out on a grid consisting of 20250 cells. The time step was 5000 seconds. The thermophysical properties of the soil specified were used for the calculation: soil density 2000 kg·m -3 , heat capacity 2000 J·(kg·K) -1 , thermal conductivity coefficient 0.8 W·(m·K) -1 , latent heat of melting 350000 J·kg -1 , freezing point 273.14 K, thaw temperature 273.16 K.
The boundary condition was the harmonic function of the temperature approximation on the soil surface at the upper boundary of the computational area: At the lower-the heat flux considered the geothermal gradient equal to a constant value of 0.06. Figure 2 shows the evolution of temperature at different soil depths over several years of observation. The graphs show good agreement between the simulation results and the experimental data.

Results and discussion
The calculations of the process of thaw permafrost during well drilling considered the circulation of the drilling mud were carryed out on the basis of the developed method. Real thermophysical parameters (density, specific heat capacity, thermal conductivity) of permafrost and usual parameters of the well (construction) and drilling (flow rate, temperature, rheological parameters of drilling mud) were used. Figure 3a shows the usual well geometry considered in the calculations. The calculation was carried out in an axisymmetric setting. The radius of the computational area was 100 m, the depth was 500 m. The grid consisting of 425,000 cells was used for the calculation. Figure 3b shows a fragment of the computational grid. The grid was thickened to the walls of the drill pipe and well.
A constant value of the drilling mud flow rate was set at the entrance to the well. The flow rate was varied from 10 to 80 l·s -1 . The solution temperature ranged from +20 to +5°C. The temperature of permafrost was chosen -2ºC.
The maximum time of the drilling in the calculation was taken 1 month. There were approximation that the drilling mud fed into the well permanently and drilling didn't stop. a b Figure 3. The well geometry (a) and a fragment of the grid around the drill pipe and the well wall (b).  Figure 4 shows usual velocity profiles in the drill pipe and annular gap, and the temperature distribution of the drilling mud and the rocks adjacent to the well. Apparently the temperature of the drilling mud in the drill pipe and the annular channel was almost uniform. The maximum temperature difference between the drilling mud in the center of the pipe and the well wall didn't exceed 1°C. a b c Figure 4. Velocity profile (a), the rate modulus (b) and temperature distribution (c) in the well at a depth of 300 m through the week from the start of drilling.
The temperature distribution over the radius for different drilling durations is shown in figure 5. At the boundary of the well wall and rock temperature changes sharply. Analysis of the results showed that during the drilling week with a usual flow rate of 30 l·s -1 the radius of the thaw zone was about 0.68 m. This was consistent with measurements of the size of the cavities after drilling in permafrost.  Further the effect of the drilling mud flow rate injected into the well during drilling was investigated. The mud temperature was set to +10°C. The quantity of heat entering the well with the drilling mud increased with the rise of drilling mud flow rate. The rate of the thaw process also increased. This data is shown in figure 6. The analysis of the results showed that with the increase of the drilling mud flow rate from 30 to 50 l·s -1 , the radius of the thaw zone increased from 0.68 to 0.93 m. This change was caused by an increase in the heat transfer coefficient from the drilling mud to the well walls with increasing flow rate.
The study results of the effect of drilling mud rheology on the rate of splitting process are shown in figure 7. Three different polymer-based drilling muds were considered, the rheology of which was described by the power model. It was shown as a result of modeling that with increasing of the consistency parameter K of the drilling mud the radius of the thaw zone was significantly reduced. This was due to a decrease in the heat transfer coefficient with an increase in the effective viscosity of the drilling mud. Thus the use of drilling muds with higher effective viscosity was preferable in terms of reducing the negative impact on the thaw processes during drilling. However it is important to understand that higher viscosity lead to additional heating of the drilling mud due to friction in the bit channels and the pump. Because it is necessary to select the optimal value of the effective viscosity. This requires further research. Figure 7. The temperature radius distribution at a depth of 300 meters in a week from the start of drilling at different drilling mud rheology at flow rate of 30 l·s -1 .

Conclusion
A mathematical model conjugate heat exchange of the well consideration phase transitions and drilling fluid circulation during to drilling was developed. The mathematical model takes the real properties of the drilling mud (viscosity, rheology), the actual parameters of the drilling (flow rate of mud and the rotation speed of the drill pipe) and easily adapts the geometry of the well. The developed method was tested on the available experimental data. It was shown good agreement with the experiment. The calculated study of the permafrost thaw around the well during drilling. The effect of the flow rate and rheological properties of the drilling mud on the dynamics of the advance of the front of thaw was investigated. It was shown that the control of rheological characteristics of drilling muds can significantly reduce the process rate of thaw during drilling in permafrost.