Coupled THM Modelling of Wellbore Stability with Drilling Unloading, Fluid Flow, and Thermal Effects Considered

Both overbalanced drilling and underbalanced drilling will lead to the change of pore pressure around wellbore. Existing research is generally based on hydraulic-mechanical (HM) coupling and assumes that pore pressure near the wellbore is initial formation pressure, which has great limitations. According to the coupled theory of mixtures for rockmedium, a coupled thermal-hydraulicmechanical (THM) model is proposed and derived, which is coded with MATLAB language and ABAQUS software as the solver. Then the wellbore stability is simulated with the proposed model by considering the drilling unloading, fluid flow, and thermal effects between the borehole and the formation.The effect of field coupling on pore pressure, stress redistribution, and temperature around a wellbore has been analyzed in detail. Through the study of wellbore stability in different conditions, it is found that (1) for overbalanced drilling, borehole with impermeable wall is more stable than that of ones with permeable wall and its stability can be improved by reducing the permeable ability of the wellbore wall; (2) for underbalanced drilling, the stability condition of permeable wellbore is much higher than that of impermeable wellbore; (3) the temperature has important influence on wellbore stability due to the variation of pore pressure and thermal stress; the wellbore stability can be improved with cooling drilling fluid for deep well. The present method can provide references for coupled thermal-hydraulic-mechanical-chemical (THMC) process analysis for wellbore.


Introduction
Boreholes, as the access for the development of oil, gas, and geothermal energy, as well as deep geological storage of CO 2 , experience instability phenomena such as sloughing and borehole wall fracturing during drilling process, which arise from removal of the original supporting rock and the interaction between the drilling fluid and formation [1][2][3]. Wellbore stability is an important factor considered in drilling engineering, which seriously affects the quality of the well and even leads to drilling failure. Wellbore instability always has been a major problem in petroleum geomechanics and engineering.
The causes of wellbore instability can be classified into either mechanical (i.e., failure of the rock around the borehole because of high stresses, low rock strength, or inappropriate drilling practice) or chemical effects which arise from damaging interaction between the rock and the drilling fluid [4][5][6]. In recent years, wellbore stability analysis has mainly focused on the chemical interaction between rock and drilling fluid. It is generally believed that the wellbore instability by hydration also is caused by the stress around wellbore exceeding a certain limit value due to the decreasing of rock strength [6][7][8]. Therefore, the accurate calculation of the stress state near wellbore is of great significance for judging whether the wellbore is unstable or not. With the development of analysis models for wellbore stability, it was found that the formation fluid, drilling fluid properties, and temperature have a great influence on the wellbore stability [9][10][11]. Mechanical behavior, pore pressure, and heat flux all 2 Mathematical Problems in Engineering influence wellbore stability, while instances of instability are a result of a combination of coupled effects. Analysis involves studying the interactions among changes of related effective stress, temperature, and pore pressure during drilling process. In the process of overbalanced drilling, the pressure of drilling fluid column is greater than formation pressure which would cause dynamic and static water loss of drilling fluid into formation. The stress state near a wellbore is affected seriously due to drilling fluid into formation. For air or foam drilling (underbalanced), the drilling fluid column pressure can be lower than formation pressure, which may cause formation fluid flow into wellbore, resulting in the change of pore pressure and stress state near wellbore. In geothermal environments, drilling fluid circulation will lead to relatively greater disturbance near the well, resulting in higher thermal stress, which will affect wellbore stability. More and more analytical models on wellbore stability issues have focused on coupling mechanisms, but these studies generally assume that the pore pressure around wellbore is the initial pore pressure [12][13][14]. When the formation is drilled, the equilibrium state is disturbed, causing the stress redistribution around borehole. The drilling is essentially an excavation unloading process [15]. For numerical analysis of wellbore stability, most studies have focused on HM coupling, elastic behavior, and overbalanced drilling, without considering the influence of drilling unloading process, stress redistribution, filter cake quality, and permeability change of the borehole wall [16,17], which is not consistent with the real case.
In this article, thermal effect and drilling unloading process will be introduced into wellbore stability analysis in a numerical modeling framework to make such analyses more general. According to the previous studies, the THM coupled mathematical model of rock medium has been proposed and the strategy to solve this model is discussed. Then, the drilling unloading, fluid flow, and thermal effect of wellbore are analyzed in overbalanced drilling and under balanced drilling, and the influence of drilling fluid and thermal effect on pore pressure and stress redistribution around borehole with time are all discussed.

Coupled THM Modelling Framework for Wellbore Stability
. . e Basic Idea of Excavation Unloading. Before the formation is drilled, the stress state of rock mass is mainly controlled by in-situ stress, and the stress state in the formation is in equilibrium. During the drilling process, the initial stress equilibrium in the formation is disturbed, causing the stress redistribution around the surrounding rock until a new balance is reached (as shown in Figure 1). The change of stress and displacement field in the formation caused by drilling is essentially an unloading process. After the formation is drilled, the excavated rock in the hole is replaced by drilling fluid, and the stress state near the wellbore is redistributed and the mechanical characteristics of the rock around the wellbore will inevitably be affected.
For unloading process modelling, external boundary loading method and excavation unloading method are usually used [6,7,15]. The displacement obtained by the external boundary loading method is actually a comprehensive result of the initial stress field and the excavation boundary, which has no obvious practical significance, while the excavation unloading method truly reflects the displacement changes caused by excavation. Correct modelling the unloading process is the key to analyze the wellbore stability in drilling engineering. The element removing and reactivating technique is an ideal method to realize wellbore unloading process, which is realized by modifying the stiffness of element set. For element removing process, the element is not really removed, but instead the stiffness of the element is multiplied by a very small coefficient (usually 1×10 −6 ). The mass of the removed element and other parameters are set to 0; thus the load and strain of the element is equal to 0.
Before the elements are removed, the mechanical equilibrium equation is defined as [15]

{P} = [K] {U}
(1) where [K] is stiffness matrix of element, {U} is node displacement array, and {P} is the load array. Introducing the parameter to control the element removing and reactivating, the modified equilibrium equation is given where . . ermal Transfer Model of Rock. Assuming the porosity of a rock element is , the proportion of solid skeleton is 1 − . According to the principle of conservation of energy and Fourier's law [16,17], the energy conservation equation of the solid skeleton is as follows: Mathematical Problems in Engineering   3 where s , s , and s are the thermal conductivity, specific heat, and density of the solid skeleton, respectively, s is the energy generated per unit volume of skeleton per unit time, is the temperature, 0 is the initial temperature, is linear expansion coefficient, is the volume modulus of rock medium, and V is volume strain. The energy conservation equation of the fluid phase can be established in a similar way. Since the fluid flow velocity in the formation is small, the kinetic energy of the fluid can be ignored. The energy equation considering convection and conduction is given as where , , and are thermal conductivity, specific heat, and density of the fluid respectively, is the energy generated per unit volume of fluid per unit time, is Darcy flow velocity, and ( ∇) is convection item, which indicates the change rate of temperature caused by the movement of fluid particles.
According to theory of mixtures [18], the total energy equation of the rock medium is superimposed by the proportion of material component.
. . Coupled Flow Model. Combined with Darcy's law, the flow equation considering rock deformation is defined as [16,17] where f is fluid viscosity, is permeability of rock, is pore pressure, = f + (1 − ) s is comprehensive compressibility, s is the compressibility of solid skeleton, f is the compressibility of fluid, = f + (1 − ) s is comprehensive thermal expansion coefficient, s = 3 l is volumetric thermal expansion coefficient of solid skeleton, f is volumetric thermal expansion coefficient of fluid, is Biot coefficient, and f is internal and external fluid sources.
In the process of fluid flow, the permeability is related to the change of porosity, and the change of pore volume is related to pore pressure and stress, so the permeability is dynamic change with fluid flow [15]. The relations between permeability, porosity, and volumetric strain are given: where 0 and 0 are initial permeability and porosity, respectively.
The relationship between permeability and volume strain is shown in Figure 2. The smaller the porosity, the more obvious is the effect of deformation on permeability.
. . Coupled Mechanical Model. According to Biot's consolidation theory, the effective stress can be expressed as where is the effective stress, is total stress, and is Kronecker symbol. In this study, tensile stress is positive and compressive stress is negative.
To describe the thermal expansion and plastic deformation of a rock medium, the total stress can be expressed in the form of increment: Mathematical Problems in Engineering  Equation (11) can be rearranged as where , is displacement components of rock, and ep is elastic-plastic stiffness matrix.
Drilling excavation causes stress redistribution around wellbore, and the zone where the stress state changes is usually called the excavation disturbed zone (as shown in Figure 3). That drilling changes the original equilibrium state of the formation would cause stress concentration around the wellbore and may lead to wellbore instability. Wellbore instability is generally divided into two types: one is collapse caused by shear failure and another is the fracturing caused by tensile stress.
Due to the simplicity and practicability of the Mohr-Coulomb criterion, it is the most commonly used for evaluating wellbore stability (as shown in Figure 4). The shear Mohr-Coulomb criterion is generally written in terms of stress invariants [19]: where m , , , denote the average stress, equivalent stress, cohesion, and friction angle, respectively. is a function of the Lode angle and friction angle .
In (13), m , and are defined by where 1 , 2 , 3 are the three principal stresses and 2 is the second invariant of the deviatoric stress. The tensile Mohr-Coulomb criterion is defined in terms of stress invariants [15]: where t is the tensile strength.
Considering both tensile strength and shear strength of a formation, the yield surface is a combination of two different failure criteria: shear failure and tensile failure. The modified Mohr-Coulomb (MMC) function is defined as follows (as shown in Figure 5): where 0 ≤ ≤ 1 is a parameter reflecting the tensile strength of formation. When = 0 indicating a higher tensile strength, the MMC yield criterion becomes the shear Mohr-Coulomb criterion; = 1 denotes that there is no tensile strength. In addition, the parameter can smooth the vertex of the yield surface and avoid the numerical divergence and slow convergence.
The plastic potential function can be defined from the yield function as follows: where is the dilatancy angle. Similarly, is a function of the Lode angle and the dilatancy angle . Associated flow occurs if = , while nonassociated flow occurs if 0 ≤ < .
According to the MMC criterion, the quantitative evaluation indices for wellbore collapse and rupture are given as [14] = where , are normal stress and shear stress on the shear surface. The wellbore would collapse for ≤ 1 and fracture for ≤ 1.

Solution Strategy
The decoupled solution method is applied for the solution of this THM coupling model. The HM coupling module and thermal module in ABAQUS software is used to solve the stress field, flow field, and temperature field. The decoupled solution between two modules is performed in MATLAB. Then, the THM coupling software for rock medium can be developed using MATLAB language as the platform and ABAQUS software as the solver, in which different parts such as thermal, hydraulic, and mechanical modules can be called. The software framework for THM is shown in Figure 6. The solution method has been verified through a typical example of one-dimensional thermal consolidation [20,21]. The decoupled method is also efficient and easy to be developed, which can be applied to the analysis of thermalhydraulic-mechanical-chemical (THMC) coupling for rock. The programming process is given as follows.
(1) The main program of MATLAB is compiled, and the calculation parameter setting and control precision of cycle iteration are defined.
(2) Different physical fields are analyzed by ABAQUS solver in sequence: write command file ( * . inp) for different physical field based on ABAUQS and USDFLD subroutine file ( * . for) for parameter evolution such as permeability change with strain based on FORTRAN. Call ABAUQS solver to solve different field modules by the command code (system) in MATLAB.
(3) Compile ABQMAIN subroutine for data conversion: the field variable node data file ( * . fil) is formed for other physical field modules to call as field variable values.
(4) The coupled parameters in the command file ( * .inp) of ABAQUS is modified and updated in MATLAB.
(5) Compile the coupled cycle processing subroutine in MATLAB: call the updated command file ( * .inp) and ABAQUS solver to calculate. Then, read the current result file ( * .dat) of ABAQUS and compare it with the previous result file until the convergence tolerance error is satisfied and the cyclic iteration is completed.

Parametric Study
During drilling, fluid flow in the formation and the temperature difference between formation and drilling fluid will cause changes of pore pressure and temperature around the   wellbore, which will inevitably affect wellbore stability. In this section, a conceptual numerical model is applied to study the wellbore stability influenced by flow conditions and temperature.
. . Modelling Approach. A plane strain model was studied based on the proposed method (as shown in Figure 7). The radius of the wellbore is 0.108 m. Due to the symmetry of the model, a 1/4 portion with the length and the width of 5 m was established.
The overburden strata stress V , horizontal maximum total stress , the minimum total stress ℎ , and the formation pressure are 69MPa, 75MPa, 54MPa, and 45Mpa, respectively. The fluid column pressure of drilling fluid is 50MPa for overbalanced drilling and 40 for underbalanced drilling. The properties of the rock and the fluid are listed in  The numerical calculation process of wellbore stability includes three steps, which are defined as follows. (1) In first step, in-situ stress, initial temperature and pore pressure are applied to the rock mass to simulate the undisturbed state of formation (Figure 8(a)). (2) To simulate the drilling unloading, the element removing technique is used to deal with the excavated part of wellbore (Figure 8 The drilling fluid is injected and fluid temperature, pore pressure, and fluid pressure surface load are applied on the wall of borehole (Figure 9).
. . Effect of Filtrate Cake Quality on Wellbore Stability. For overbalanced drilling, the drilling fluid column pressure is greater than the formation pressure, and the drilling fluid can form a mud filtrate cake on the wellbore (as shown in Figure 10). The fluid flow depends greatly on the differential pressure and filtrate cake quality, which are the direct causes of fluid flow and pore pressure changes around wellbore. Therefore, the wellbore stability is influenced by filtrate cake quality.
Two kinds of extreme conditions for filtrate cake quality were discussed. One is impermeable wall or closed wellbore which indicates that the filtrate cake is impermeable and thus there is no fluid communication between the wellbore and the formation. Another scenario is permeable wall condition which means drilling fluid and formation is connected. To study the influence of the filtrate cake quality on wellbore  stability, the temperature of drilling fluid is assumed the same as that of the formation in this section.
Compared with the initial stress state (Figure 11), the effective radial stress near the wellbore decreased rapidly after drilling unloading, and the stress continues to decrease gradually due to pore pressure change with time, as shown in Figure 12. For a permeable wall condition, the pressure support on the wellbore by drilling fluid will be reduced by seepage, and the effective radial stress at the wellbore decreases gradually until to zero ultimately. As for an impermeable wall condition, the hydraulic connection between the drilling fluid and the formation water is eliminated by the closure of good filtrate cake, and then the liquid column 8 Mathematical Problems in Engineering S, S11 (Avg: 75%) -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 -3.000e+07 Y X (a) S, S22 (Avg: 75%) -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.000e+06 -9.  of drilling fluid can provide an effective supporting for the wellbore. The effective radial stress on the inner wall is -6.5 MPa after 24 hours for impermeable wall, whereas it is 0 MPa for permeable wall. Figure 13 shows the effective hoop stress distribution near wellbore region. The maximum hoop stress is at the location 0.5 from the wall after drilling unloading. For a permeable wall, the hoop stress decreases gradually with time and shows a trend from compression to tension after 24 hours. For an impermeable wall, the hoop stress is gradually reduced and stays in a compressive stress state.
Due to the impact of drilling unloading, the pore pressure near wellbore decreases rapidly after drilling excavation (Figure 14). For a permeable wall, the minimum pore pressure is 28MPa at the location 0.4 from the wall. Since the drilling fluid in the wellbore is connected with formation fluid, the pore pressure around wellbore increases to 50MPa gradually after 24 hours. For an impermeable wall, the pore pressure around wellbore goes to initial pore pressure 45MPa gradually ( Figure 11). Figure 15 shows the distribution of quantitative evaluation indices and near wellbore region, respectively. It can be found that the wellbore is easiest to collapse when is equal to 90 ∘ or 270 ∘ , and the wellbore is easiest to fracture when is equal to 0 ∘ or 180 ∘ (the direction of the maximum principal stress is horizontal).
For an impermeable wall, drilling fluid and formation water is separated by the filtrate cake. The pore pressure near wellbore changes obviously by drilling unloading. After that, the pore pressure near wellbore gradually tends to initial pore pressure. In the condition of overbalanced drilling, the collapse resistance of the wellbore is improved with the supporting effect by liquid column pressure on the borehole wall.
For a permeable wall, the drilling fluid is connected with the formation fluid. Due to the drilling liquid column as the inner boundary of seepage field, the supporting effect by liquid column pressure on the formation is gradually reduced to zero, while the effective load on far field boundary is the difference between initial total stress and pore pressure. The pore pressure gradient is gradually decreased along the direction of far field which can make a certain inhibitory effect on the wellbore deformation, but the reduction of wall support load promotes the wellbore deformation. In the condition of overbalanced drilling, the pore pressure of borehole wall gradually tends to the pressure of liquid column, which is higher than that of far field. With the decrease of effective stress, the capacity of the wellbore to resist tension failure decreased. Due to the effect of pore fluid pressure gradient and reduction of wellbore supporting load, both radial stress and hoop stress of permeable wellbore are lower than those of impermeable wellbore. Thus, wellbore stability of impermeable wall is higher than that of permeable wellbore, so it is particularly important to improve the filter cake quality of drilling fluid under the condition of overbalanced drilling.

. . Effect of Fluid Flow on Wellbore Stability for Underbalanced Drilling.
In the process of underbalanced drilling, formation fluid can flow into the borehole, and the stress around the wellbore will be redistributed with the seepage of formation fluid, thus affecting the wellbore stability. Two kinds of wellbore conditions were analyzed in this section: one is impermeable wall of wellbore another scenario is permeable condition. The radial stress near wellbore decreased rapidly after drilling unloading, as shown in Figure 16. After that, the seepage effect makes the stress decrease continuously. For a permeable wall, the supporting load of wellbore wall is decreased by the seepage, and the effective radial stress of rock around wellbore decreased gradually and tends to 0 MPa. For an impermeable wall, the pressure of liquid column is lower than the initial pore pressure of formation, and the deformation of wellbore makes the stress state of rock into tension. Figure 17 shows the hoop stress distribution around wellbore. The maximum hoop stress is at the location 0.4 from the wall. The hoop stress keeps in compressive stress state with time. Figure 18 shows pore pressure distribution of near wellbore region. For a permeable wall, the pore pressure around wellbore tends to 40MPa after 24 hours. For an impermeable wall, the drilling fluid in wellbore cannot connect with the far pore pressure field and the pore pressure goes to initial pore pressure gradually with time.
For an impermeable wall in underbalanced drilling condition, the wellbore shrinkage and a large change of pore pressure near the wellbore occurred. The pore pressure around wellbore has the tendency to initial formation pressure and the radial stress is in tension gradually, which make wellbore stability worse. For a permeable wall, the pore pressure near wellbore wall tends to the pressure of liquid column, which is lower than the far formation pressure, and the stability of wellbore increased. Therefore, the wellbore stability with permeable wall is higher than that of impermeable wall in underbalanced drilling condition (Figure 19).
. . Effect of Temperature on Wellbore Stability. To discuss the effect of temperature on wellbore stability, different coupled relationships are illustrated through the following three cases, which are shown in Figure 20. Case 1 and case 2 are two types of incomplete coupling forms, while case 3 is fully coupled THM form. For case 1, the coupling effect between temperature field and stress field is not considered, but the influence of temperature field and stress field on seepage field is considered. For case 2, the coupling effect between temperature field and seepage field is not considered, but the influence of seepage field and temperature field on stress field is considered. Taking overbalanced drilling with good filtrate cake quality as an example, the importance of temperature coupling on wellbore stability is analyzed, in which the temperature of drilling fluid and formation are 120 degrees Celsius and 70 degrees Celsius, respectively. Figures 21 and 22 show the comparison of temperature distribution after 24 hours under different cases. Due to lower thermal conductivity and higher specific heat of fluid, the temperature field is least affected in case 1. The temperature field is most affected in case 2 for higher thermal conductivity of solid skeleton. Because the theory of mixtures is applied in this fully coupled THM model [18], the temperature field of case 3 has shown a moderate degree. Figure 23 shows the comparison of pore pressure after 24 hours under different cases. Due to the fact that the thermal-hydraulic coupling (TH) is not considered in case 2, the pore pressure difference between wellbore wall and the formation is the smallest. For case 3, the THM coupling is fully considered, so the pore pressure difference between wellbore and the formation is the largest. For stress field (Figure 24), the radial and hoop stress are basically in a compressive state, and the distributions around wellbore are similar for different cases. The temperature on stress field is affected most in case 3 and least in case 2, which is similar to pore pressure. Figure 25 shows the distribution of collapse index and fracture index under different cases. It can be found that the Mathematical Problems in Engineering For collapse failure, the collapse index under case 2 is the smallest and the wellbore is most easy to collapse. For fracture failure, the fracture index under case 2 is the smallest and the wellbore is most easy to break. There are obvious differences in stress, pore pressure, and temperature of wellbore by fullcoupled model and partial-coupled models. During drilling, the temperature difference between the drilling fluid and the formation will lead to temperature change near wellbore. The disturbed range of temperature is related to the thermal diffusivity of the formation and the fluid inside, which leads to the great change of pore pressure and effective stress in the near-wellbore zone, and inevitably affects the wellbore stability. Therefore, it is necessary to take full-coupled model into account when dealing with the evaluation of wellbore stability [22]. The influence of temperature fluctuation on wellbore stability is studied by changing drilling fluid temperature using THM fully coupled model. The results are shown in Figures 26 and 27. When the formation is heated by the drilling fluid, the larger the temperature difference, the smaller the collapse index and fracture index of wellbore, which indicates that the instability possibility of wellbore is greater. When the drilling fluid temperature is lower than the formation temperature, the stress around wellbore tends to decrease with the increase of temperature difference. The larger the temperature difference is, the smaller the instability possibility of wellbore occurs.    Normalized radial distance r/r w Figure 23: Comparison of pore pressure distribution along BC. . . Discussion. A well located in the Jidong oilfield of China is selected as an example to investigate the wellbore stability in a mudstone formation. According to laboratory tests, field data, logs, and related geological data, the calculation parameters at the depth of 4100m are defined as overburden strata stress 91.6MPa, the maximum horizontal stress 81.2MPa, the minimum horizontal stress 68.7MPa, and the formation pressure 42.5MPa; the elastic modulus of rock 20.2GPa, Poisson's ratio 0.16, cohesion 24.1MPa, and the internal friction angle 21.7 ∘ ; permeability of formation 1.01mD; and the porosity 9%. The thermal expansion of rock media is 5.4×10 −5 , specific heat 886 J/kg• ∘ C, and thermal conductivity 3.25 J/(m• ∘ C). In drilling stage, the temperature of drilling fluid is 10 degrees Celsius lower than that of the formation.

Mathematical Problems in Engineering
The failure process of wellbore is simulated with the drilling fluid density of 1.1, 1.2, and 1.3 g/cm 3 , respectively. The distribution of damaged zone caused by drilling unloading is shown in Figures 28 and 29. In the drilling unloading, the drilling time is short, and thus the borehole stability is mainly controlled by the liquid column pressure of drilling fluid. When the drilling fluid density is 1.1, the failure 8 depth of wellbore is about 0.07m in the horizontal minimum principal stress (EF line) with a maximum equivalent plastic strain of 4.574×10 −3 after drilling unloading. With the increasing of drilling fluid density, the failure depth decreases gradually. When the drilling fluid density is 1.3, the failure depth of wellbore is about 0.02m. When the drilling fluid density is greater than 1.35, the rock mass around the wellbore is basically in an elastic state and has reached a stable state, which is consistent with that of field testing.
According to the failure depths in horizontal and vertical direction, the enlargement rate of wellbore can be defined as ( Figure 3)   Mathematical Problems in Engineering where 0 is the radius of the bit and = ( a + b )/2 is the average radius of wellbore after drilling. The traditional analytical model for collapse pressure prediction is given as [12] where is the formation pressure of = (45 − /2), = (1 − 2 )/2(1 − ), and is the seepage ability coefficient of filtrate cake. = 1 for a permeable wellbore wall, and = 0 for an impermeable wall.
Although the influence of temperature and seepage ability can be considered in the analytical models (Eq. (21)), it only provide a static value that cannot consider the pore pressure and effective stress change with time as a transient response. In the traditional model, the wellbore is thought instability once shear failure occurs around wellbore. In practical engineering, the wellbore can be allowed certain local failure and the enlargement rate of wellbore can be controlled within 15% if the rock debris can be carried from the bottom of the well in time [23]. In this numerical model, the initial enlargement rate of wellbore is 10.8% when the drilling fluid density is 1.3. As for the drilling fluid density of 1.2 and 1.1, the wellbore enlargement rate is 18.1% and 32.5%, respectively. The actual density of drilling fluid used in this well is 1.3, there is only slight sloughing in the construction process, and the stability condition of borehole meets the requirements of drilling construction, which is consistent with the numerical results.
According to field data, the studied wellbore is stable in the drilling stages, but local instability of wellbore appears after 2 days. The reason is that the open-hole section is too long (500m) in mudstone formation and thus the immersion time of wellbore is increased by drilling fluid. Thus, the hydration effect cannot be neglected for long time drilling. The hydration effect causes the decrease of rock strength and the increase of collapse pressure. The drilling fluid density of 1.3 at this early stage can support the wellbore to make wellbore stable, but the drilling fluid density must be increased to 1.59 after 2 days for keeping the wellbore stable in engineering field. Therefore, the wellbore instability is a dynamic process, and the collapse pressure is also a dynamic value, which cannot be predicted by the traditional analytical model.
According to the laboratory tests [24,25], the hydration effect causes the surrounding rock to be deteriorated. The initial water content of this mudstone formation is 2%, and the saturated water content is 10% after long time immersion by drilling fluid. The evolution of strength parameters is given: where 0 and 0 are the cohesion and friction angle at initial water content 0 , respectively, and is current water content.
The evolution of strength parameters of rock is coded by USDFLD (user subroutine to redefine field variables at a material point) subroutine. The enlargement rate curves of wellbore considering the effect hydration are shown in Figure 30. It can be found that the enlargement rate increases with the increase of drilling fluid immersion time, and the collapse pressure increases accordingly.
Wellbore stability is affected by in-situ stress, drilling fluid, temperature change, wellbore seepage, hydration effect, creep, etc. According to the above analysis, if the immersion time of wellbore is too long, the hydration and creep effects of wellbore need to be considered, which should be coupled with three fields of THM model. The complete coupling between hydration, creep, and THM needs to be further studied theoretically and experimentally.

Conclusions
Based on continuum mechanics and the theory of mixtures, a THM coupling model of rock medium is established, and then it is coded with MATLAB language as platform and ABAQUS software as the solver. Considering the drilling unloading, the numerical model for wellbore analysis is established. The variation laws of temperature, pore pressure, and stress around wellbore under different physical field coupling are compared. The results show that the stress, pore pressure, and temperature around wellbore obtained by full-coupling and partial-coupling are quite different, which prove that the fully coupled model is reasonable for wellbore stability analysis.
The change of pore pressure near wellbore is large at the moment of drilling unloading. For overbalanced drilling, wellbore stability with permeable wall is lower than that of impermeable wall. Under the condition of underbalanced drilling, the wellbore shrinkage will occur in impermeable wall which makes the wellbore in tension state and reduces the wellbore stability, but the wellbore stability will be improved by the connectivity of permeable wall. Under the condition of overbalanced drilling, something should pay attention to the filtrate cake by reducing the disturbance of drill string and improving the protected capacity of drilling fluid as much as possible. For underbalanced drilling, the drilling fluid should be selected to delay the formation of filtrate cake.
As for designing drilling fluid for deep wells and wells with large temperature gradient, the variation of pore pressure and thermal stress caused by temperature change should be taken into account. When the formation temperature is lower than the drilling fluid, the possibility of wellbore failure increases with the increase of temperature difference. If the drilling fluid temperature is lower than that of the formation, the larger the temperature difference is, the smaller the instability possibility of wellbore occurs. In deep well drilling, if equipped with cooling drilling fluid equipment, it can increase wellbore stability, but also conducive to the normal operation of logging and downhole instruments. The method proposed in this study provides an effective way to analyze the THM coupling process for wellbore and lays a foundation for further study of the THMC coupling problem on wellbore stability.

Data Availability
The data used to support the findings of this study are included within the article.

Disclosure
Caoxuan Wen is a co-first author.

Conflicts of Interest
There is no conflict of interests regarding this paper.