Fully Coupled Thermo-hydro-mechanical Model for Wellbore Stability Analysis in Deep Gas-Bearing Unsaturated Formations Based on Thermodynamics

Abstract

A fully rigorous unsaturated thermo-hydro-mechanical (THM) coupling model based on non-equilibrium thermodynamics and continuum mechanics is developed to investigate the instability mechanism of the borehole drilled in deep gas-bearing formations. It incorporates the gas–water two-phase seepage phenomena, geomechanical loading, and temperature-difference disturbance during deep drilling process. Comparing the predictions of the previous saturated THM model and the presented model, the maximum relative errors of pore pressure, effective radial stress, effective tangential stress, and temperature are 30.2%, 34.6%, 3.4%, and 0.5%, respectively. The maximum prediction deviation in the width and depth of the failure zone occurs at the early time (t = 10⁻⁵ d). The unsaturated thermal coupling effect, although less influential on pore pressure and effective radial stress, makes a great contribution to the effective tangential stress with a maximum relative error of 10.1% for the temperature-difference of 30 K, and enhances the time dependence of the failure zone. The cooling condition is favorable to avoid shear collapse failure of the borehole, while the heating condition is not. The results also show that the gas-bearing formations with lower initial water saturation, higher rock permeability, or higher water phase relative permeability correspond to the higher risk of wellbore instability. The thermal conductivity and thermal expansion coefficient of the solid phase have a much greater effect on wellbore stability than those of the fluid phases. The effect of unsaturated THM coupling on wellbore stability may be more significant at lower rock permeabilities. The research findings provide insight into the wellbore stability analysis during drilling in deep gas-bearing formations.

Highlights

• A fully rigorous unsaturated thermo-hydro-mechanical (THM) coupling model based on thermodynamics is developed.

• The unsaturated THM responses of pore pressure, stress components, and failure zone of a borehole were clarified.

• The prediction differences between the presented model and previous saturated THM and unsaturated HM models were quantitatively analyzed.

• The influence of two-phase seepage and thermal transport on wellbore stability were discussed parametrically.

Appendices

Appendix 1

1.1 Balance Equations

1.1.1 Mass Balance Equations

The balance equations based on the mixture coupling theory for a thermodynamic open system have recently been widely recognized, discussed, and extended, and the details can be referred to found in previous literature (Heidug and Wong 1996; Chen et al. 2013; Coussy 2004), will not be repeated here. For the porous medium mixture of interest containing solids and fluids (including liquid and gas phases), assuming no mass exchange between components, the mass balance equation is generally expressed as

$$\frac{{\rm D}}{{\rm D}t}\left( {\int_{V} {\rho^{\xi } {\rm d}V} } \right) = - \int_{\Gamma } {{\varvec{I}}^{\xi } \cdot {\varvec{n}}{\rm d}\Gamma } ,$$

(50)

where the superscript ξ = s, l and g represent solid, liquid, and gas, respectively, ρ^ξ denotes the mass density of each component, n is the outward normal vector, and the component mass flux I^ξ can be given by

$${\varvec{I}}^{\xi } = \rho^{\xi } \left( {{\mathbf{v}}_{\xi } - {\mathbf{v}}_{{\rm s}} } \right),$$

(51)

where v_ξ represents the velocity of liquid (ξ = l) and gas (ξ = g), v_s denotes the solid velocity, and the material time derivative following the motion of the solid is expressed as

$$\frac{{\rm D}\left( \cdot \right)}{{{\rm D}t}} = \frac{\partial \left( \cdot \right)}{{\partial t}} + {\mathbf{v}}_{{\rm s}} \cdot \nabla \left( \cdot \right).$$

(52)

Considering that there is no solid mass flux into V, while it is open to fluid mass flow (liquid and gas), by substituting Eqs. (51) and (52) into Eq. (50), the localized version of the mass balance equations for solid and fluid components can be derived as

For solid phase (ξ = s):

$$\dot{\rho }^{{\rm s}} + \rho^{{\rm s}} \nabla \cdot {\mathbf{v}}_{{\rm s}} = 0.$$

(53)

For liquid and gas phase (ξ = l and g):

$$\dot{\rho }^{\xi } + \rho^{\xi } \nabla \cdot {\mathbf{v}}_{{\rm s}} + \nabla \cdot {\varvec{I}}^{\xi } = 0,$$

(54)

where the time derivative ∂(·) is replaced by the over dot ‘·’.

The mass density of each component ρ^ξ is defined with respect to the total volume of the fluid–solid mixture system, and can be expressed in terms of respective volume fraction φ_ξ and true state density ρ_ξ a as

$$\rho^{\xi } = \varphi_{\xi } \rho_{\xi } ,$$

(55)

where the component volume fraction φ_ξ is related to the total porosity φ of the porous medium through

$$\varphi_{{\rm s}} = 1 - \varphi , \, \varphi_{{\rm l}} = \varphi S_{{\rm l}} , \, \varphi_{{\rm g}} = \varphi S_{{\rm g}} ,$$

(56)

where S_ξ is the saturation of liquid (ξ = l) and gas (ξ = g).

The following density relationship is then obtained:

$$\rho^{{\rm s}} = \left( {1 - \varphi } \right)\rho_{{\rm s}} , \, \rho^{{\rm l}} = \varphi S_{{\rm l}} \rho_{{\rm l}} , \, \rho^{{\rm g}} = \varphi S_{{\rm g}} \rho_{{\rm g}} ,$$

(57)

in which S_l + S_g = 1.

1.1.2 Heat Balance Equation

In the REV, the thermal density, is influenced only by the heat flow flux across the boundary Γ, and can be defined in a unified form for different components as

$$q_{\xi } = \rho^{\xi } C_{\xi } T.$$

(58)

Substituting Eq. (57) into Eq. (58), one can obtain

$$q_{{\rm s}} = \left( {1 - \varphi } \right)\rho_{{\rm s}} C_{{\rm s}} T, \quad q_{{\rm l}} = \varphi S_{{\rm l}} \rho_{{\rm l}} C_{{\rm l}} T, \quad q_{{\rm g}} = \varphi S_{{\rm g}} \rho_{{\rm g}} C_{{\rm g}} T,$$

(59)

where q_s, q_l, and q_g are the heat density of solid, liquid, and gas, respectively, C_s, C_l, and C_g represent the specific heat capacity of solid, liquid, and gas, respectively, and T is the temperature.

The total heat flow q includes two aspects (Ma et al. 2022c; Siddiqui and Roshan 2022):

• (1) The heat flow contained in the liquid and gas flow, which are defined as h^wI^w and h^gI^g, where h^w and h^g represent the enthalpy of the liquid and gas, respectively.

• (2) The reduced heat flow q′ = q-h^wI^w − h^gI^g, which means the difference between the total heat flow and the heat flow carried by the fluids (Katchalsky and Curran 1965).

Therefore, the heat balance equation can be obtained based on the balance law for the thermodynamic open system as

$$\frac{{\rm D}}{{\rm D}t}\int_{V} {\left( {q_{{\rm s}} + q_{{\rm l}} + q_{{\rm g}} } \right)} {\rm d}V = - \int_{\Gamma } {\left( {{{\varvec{q}}^{\prime}} + h^{{\rm l}} {\varvec{I}}^{{\rm l}} + h^{{\rm g}} {\varvec{I}}^{{\rm g}} } \right)} \cdot {\varvec{n}}{\rm d}\Gamma .$$

(60)

Using Eq. (52) and neglecting any other heat source or heat change, Eq. (60) can be rewritten as

$$\left( {q_{{\rm s}} + q_{{\rm l}} + q_{{\rm g}} } \right)^{ \cdot } + \left( {q_{{\rm s}} + q_{{\rm l}} + q_{{\rm g}} } \right)\nabla \cdot {\mathbf{v}}_{{\rm s}} + \nabla \cdot {{\varvec{q}}^{\prime}} + \nabla \cdot h^{{\rm l}} {\varvec{I}}^{{\rm l}} + \nabla \cdot h^{{\rm g}} {\varvec{I}}^{{\rm g}} = 0.$$

(61)

1.1.3 Helmholtz Free Energy Balance Equation

First, involving in the concepts of the internal energy density ε and entropy density η, the Helmholtz free energy density ψ is conveniently defined as (Chen et al. 2013)

$$\psi = \varepsilon - T\eta .$$

(62)

The local form of the internal energy balance equation can be written as

$$\dot{\varepsilon } + \varepsilon \nabla \cdot {\mathbf{v}}_{{\rm s}} - \nabla \cdot \left( {{\varvec{\sigma}}{\mathbf{v}}_{{\rm s}} } \right) + \nabla \cdot {{\varvec{q}}^{\prime} + }\nabla \cdot h^{{\rm w}} {\varvec{I}}^{{\rm w}} + \nabla \cdot h^{{\rm g}} {\varvec{I}}^{{\rm g}} = 0,$$

(63)

where σ denotes the Cauchy stress tensor.

Likewise, the local form of the entropy balance equation can be obtained as

$$\left( {\eta^{{\rm s}} + \eta^{{\rm l}} + \eta^{{\rm g}} } \right)^{ \cdot } + \left( {\eta^{{\rm s}} + \eta^{{\rm l}} + \eta^{{\rm g}} } \right)\nabla \cdot {\mathbf{v}}_{{\rm s}} + \nabla \cdot {\varvec{I}}_{\eta } - \gamma = 0,$$

(64)

where η^s, η^l, and g^g represent the entropy density of the solid, liquid, and gas phases, respectively, γ denotes the entropy production per unit mixture volume, and I_η means the entropy exchange with surroundings (Katchalsky and Curran 1965), and can be expressed as

$${\varvec{I}}_{\eta } = \frac{{{\varvec{q}} - \omega^{{\rm l}} {\varvec{I}}^{{\rm l}} - \omega^{{\rm l}} {\varvec{I}}^{{\rm l}} }}{T} = \frac{{{{\varvec{q}}^{\prime}}}}{T} + \eta^{{\rm l}} {\varvec{I}}^{{\rm l}} + \eta^{{\rm g}} {\varvec{I}}^{{\rm g}} ,$$

(65)

in which ω^ξ = h^ξ− Tη^ξ (ξ = l and g) represent the chemical potential of the liquid and gas.

Therefore, according to Eq. (62), using Reynold’s transport theorem and combining Eqs. (63) and (64), the balance equation for the Helmholtz free energy density can ultimately be derived as

$$\dot{\psi } + \psi \nabla \cdot {\mathbf{v}}_{{\rm s}} - \nabla \cdot \left( {{\varvec{\sigma}}{\mathbf{v}}_{{\rm s}} } \right) + \nabla \cdot {{\varvec{q}}^{\prime} + }\nabla \cdot h^{{\rm w}} {\varvec{I}}^{{\rm w}} + \nabla \cdot h^{{\rm g}} {\varvec{I}}^{{\rm g}} + \left( {\dot{T} + T\nabla \cdot {\mathbf{v}}_{{\rm s}} } \right)\eta^{{\rm mix}} - T\nabla \cdot {\varvec{I}}_{\eta } = - T\gamma \le 0,$$

(66)

where η^mix = η^s + η^l + η^g represents the entropy density of the mixture.

1.2 Entropy Production and Transport Law

It is essential to interpret the entropy production to quantify the dissipation mechanism within the rock. The entropy production is assumed to be due to the frictional resistance at the solid–fluid interface as the fluid flow in the porous medium here. Thus, using standard arguments of non-equilibrium thermodynamics, a macroscopic expression for dissipation under non-isothermal conditions is given by (Ma et al. 2022c; Siddiqui and Roshan 2022)

$$T\gamma = - {\varvec{I}}_{\eta } \cdot \nabla T - {\varvec{I}}^{{\rm l}} \cdot \nabla \omega^{{\rm l}} - {\varvec{I}}^{{\rm g}} \cdot \nabla \omega^{{\rm g}} .$$

(67)

Moreover, according to the Gibbs–Duhem equation (Moran et al. 2010), the relationship between the fluid chemical potential ω^ξ and the respective pressure p_ξ can be established based on the assumption of the local thermal equilibrium condition as

$$\rho_{{\rm l}} \nabla \omega^{{\rm l}} = \nabla p_{{\rm l}}$$

(68)

$$\rho_{{\rm g}} \nabla \omega^{{\rm g}} = \nabla p_{{\rm g}} ,$$

(69)

where p_l and p_g are the pore liquid pressure and gas pressure, respectively.

In addition, the Darcy velocity of fluids are given by

$${\varvec{u}}_{{\rm l}} = \varphi_{{\rm l}} \left( {{\mathbf{v}}_{{\rm l}} - {\mathbf{v}}_{{\rm s}} } \right)$$

(70)

$${\varvec{u}}_{{\rm g}} = \varphi_{{\rm g}} \left( {{\mathbf{v}}_{{\rm g}} - {\mathbf{v}}_{{\rm s}} } \right).$$

(71)

Therefore, neglecting the chemical potential of the liquid and gas phases, Eqs. (66) and (67) are available to rewrite the dissipation function by making the necessary rearrangements for driving forces and fluxes as

$$0 \le T\gamma = - {{\varvec{q}}^{\prime}} \cdot \frac{\nabla T}{T} - {\varvec{u}}_{{\rm w}} \nabla p_{{\rm w}} - {\varvec{u}}_{{\rm g}} \nabla p_{{\rm g}} = - {{\varvec{q}}^{\prime}} \cdot \frac{\nabla T}{T} - {\varvec{u}}_{{\rm w}} \nabla p_{{\rm w}} - {\varvec{u}}_{{\rm g}} \nabla p_{{\rm g}} .$$

(72)

Furthermore, the linear dependence of fluid fluxes and heat flux on the corresponding driving forces can then be expressed with the help of phenomenological equations as (Ma et al. 2022c; Siddiqui and Roshan 2022)

$$\rho_{{\rm l}} {\varvec{u}}_{{\rm l}} = - L_{11} \left( {\frac{{\nabla p_{{\rm l}} }}{{\rho_{{\rm l}} }}} \right) - L_{12} \left( {\frac{{\nabla p_{{\rm g}} }}{{\rho_{{\rm l}} }}} \right) - L_{13} \left( {\frac{\nabla T}{T}} \right)$$

(73)

$$\rho_{{\rm g}} {\varvec{u}}_{{\rm g}} = - L_{21} \left( {\frac{{\nabla p_{{\rm w}} }}{{\rho_{{\rm w}} }}} \right) - L_{22} \left( {\frac{{\nabla p_{{\rm g}} }}{{\rho_{{\rm g}} }}} \right) - L_{23} \left( {\frac{\nabla T}{T}} \right)$$

(74)

$${{\varvec{q}}^{\prime}} = - L_{31} \left( {\frac{{\nabla p_{{\rm w}} }}{{\rho_{{\rm w}} }}} \right) - L_{32} \left( {\frac{{\nabla p_{{\rm g}} }}{{\rho_{{\rm g}} }}} \right) - L_{33} \left( {\frac{\nabla T}{T}} \right).$$

(75)

1.3 Equations of State

Following the above Helmholtz free energy balance and dissipation equations, the solid/fluid coupling interaction response to stress, strain, pressure, and temperature can be further determined.

1.3.1 Basic Equation of State

Considering that the rock is in mechanical equilibrium, there is

$$\nabla \cdot {\varvec{\sigma}} = 0.$$

(76)

Substituting Eq. (67) into Eq. (66), thus the Helmholtz free energy balance equation is then rewritten as

$$\dot{\psi } + \psi \nabla \cdot {\mathbf{v}}_{{\rm s}} - {\rm tr}\left( {{\varvec{\sigma}}\nabla {\mathbf{v}}_{{\rm s}} } \right) + \omega^{{\rm l}} \nabla \cdot {\varvec{I}}^{{\rm l}} + \omega^{{\rm g}} \nabla \cdot {\varvec{I}}^{{\rm g}} + \left( {\dot{T} + T\nabla \cdot {\mathbf{v}}_{{\rm s}} } \right)\eta^{{\rm mix}} = 0.$$

(77)

For the small strain ε problem of interest, according to the relative concepts of continuum mechanics, invoking the fluid mass balance (Eq. 54) and making the necessary simplifications, one can obtain the referential equivalent expression of Eq. (77) as

$$\dot{\Psi } = {\rm tr}\left( {{{\varvec{\sigma}} \dot{\varvec{\varepsilon }}}} \right) + \omega^{{\rm l}} \dot{m}^{{\rm l}} + \omega^{{\rm g}} \dot{m}^{{\rm g}} - J\eta^{{\rm mix}} \dot{T},$$

(78)

where Ψ = Jψ represents the free energy in the reference configuration, m^l = JφS_lρ_l and m^g = JφS_gρ_g are the masses of the liquid and gas per unit referential volume, respectively, ε is the infinitesimal strain tensor, and J denotes the Jacobian of the deformation, the time derivative of which is given by

$$\dot{J} = J\nabla \cdot {\mathbf{v}}_{{\rm s}} .$$

(79)

1.3.2 Helmholtz Free Energy Density of the Pore Fluids

According to the classical thermodynamics, the free energy density ψ_pore of the fluids (liquid and gas) contained in the free pore space can be defined as

$$\psi_{{\rm pore}} = - p_{{\rm pore}} + \rho_{{\rm l}} S_{{\rm l}} \omega^{{\rm l}} + \rho_{{\rm g}} S_{{\rm g}} \omega^{{\rm g}} ,$$

(80)

where p_pore = p_lS_l + p_gS_g represents the average pore pressure.

Under the non-isothermal condition, where the presence of both liquid and gas is considered, based on the Gibbs–Duhem equation, there is

$$\dot{p}_{{\rm pore}} = S_{{\rm l}} \eta_{{\rm l}} \dot{T} + S_{{\rm g}} \eta_{{\rm g}} \dot{T} + S_{{\rm l}} \dot{\omega }^{{\rm l}} \rho_{{\rm l}} + S_{{\rm g}} \dot{\omega }^{{\rm g}} \rho_{{\rm g}} + \dot{S}_{{\rm l}} p_{{\rm l}} + \dot{S}_{{\rm g}} p_{{\rm g}} ,$$

(81)

where η_l and η_g represent the liquid and gas entropies per liquid and gas volume, respectively, and η^l = φS_lη_l and η^g = φS_gη_g. Combining Eq. (81) with the time derivative of Eq. (80) can finally yield

$$\dot{\psi }_{{\rm pore}} = - S_{{\rm l}} \eta_{{\rm l}} \dot{T} - S_{{\rm g}} \eta_{{\rm g}} \dot{T} + \omega^{{\rm l}} \left( {S_{{\rm l}} \rho_{{\rm l}} } \right)^{ \cdot } + \omega^{{\rm g}} \left( {S_{{\rm g}} \rho_{{\rm g}} } \right)^{ \cdot } - \dot{S}_{{\rm l}} p_{{\rm l}} - \dot{S}_{{\rm g}} p_{{\rm g}} .$$

(82)

1.3.3 Helmholtz Free Energy Density of the Solid Matrix

By subtracting the contribution of the pore fluid free energy Jφψ_pore from the total free energy Ψ of the solid/fluid mixture system, the free energy of the solid matrix can be expressed as

$$\left( {\Psi - J\varphi \psi_{{\rm pore}} } \right)^{ \cdot } = \dot{\Psi } - J\varphi \dot{\psi }_{{\rm pore}} - \psi_{{\rm pore}} \left( {J\varphi } \right)^{ \cdot } .$$

(83)

Substituting Eqs. (78), (80) and (82) into (83) yields

$$\begin{aligned} \left( {\Psi - J\varphi \psi_{{\rm pore}} } \right)^{ \cdot } & = {\rm tr}\left( {{{\varvec{\sigma}} \dot{\varvec{\varepsilon }}}} \right) + \omega^{{\rm l}} \dot{m}^{{\rm l}} + \omega^{{\rm g}} \dot{m}^{{\rm g}} - J\eta^{{\rm mix}} \dot{T} \\ & \quad - J\varphi \left( { - S_{{\rm l}} \eta_{{\rm l}} \dot{T} - S_{{\rm g}} \eta_{{\rm g}} \dot{T} + \omega^{{\rm l}} \left( {S_{{\rm l}} \rho_{{\rm l}} } \right)^{ \cdot } + \omega^{{\rm g}} \left( {S_{{\rm g}} \rho_{{\rm g}} } \right)^{ \cdot } - \dot{S}_{{\rm l}} p_{{\rm l}} - \dot{S}_{{\rm g}} p_{{\rm g}} } \right) \\ & \quad - \left( { - p_{{\rm l}} S_{{\rm l}} - p_{{\rm g}} S_{{\rm g}} + \rho_{{\rm l}} S_{{\rm l}} \omega^{{\rm l}} + \rho_{{\rm g}} S_{{\rm g}} \omega^{{\rm g}} } \right)\left( {J\varphi } \right).^{ \cdot } \\ \end{aligned}$$

(84)

And Eq. (84) can be further rearranged as

$$\begin{aligned} \left( {\Psi - J\varphi \psi_{{\rm pore}} } \right)^{ \cdot } & = {\rm tr}\left( {{{\varvec{\sigma}} \dot{\varvec{\varepsilon }}}} \right) + \omega^{{\rm l}} \dot{m}^{{\rm l}} + \omega^{{\rm g}} \dot{m}^{{\rm g}} - J\eta^{{\rm mix}} \dot{T} \, \\ & \quad + \left( {J\varphi } \right)^{ \cdot } \left( {p_{{\rm l}} S_{{\rm l}} + p_{{\rm g}} S_{{\rm g}} } \right) + J\varphi \left( {S_{{\rm l}} \eta_{{\rm l}} + S_{{\rm g}} \eta_{{\rm l}} } \right)\dot{T} + J\varphi \left( {\dot{S}_{{\rm l}} p_{{\rm l}} + \dot{S}_{{\rm g}} p_{{\rm g}} } \right) \\ & \quad - \omega^{{\rm l}} \left( {J\varphi S_{{\rm l}} \rho_{{\rm l}} } \right)^{ \cdot } - \omega^{{\rm g}} \left( {J\varphi S_{{\rm g}} \rho_{{\rm g}} } \right)^{ \cdot } . \\ \end{aligned}$$

(85)

Recalling the previous definition of Jη^mix, m^l and m^g, Eq. (85) can finally be simplified as

$$\left( {\Psi - J\varphi \psi_{{\rm pore}} } \right)^{ \cdot } = {\rm tr}\left( {{{\varvec{\sigma}} \dot{\varvec{\varepsilon }}}} \right) + p_{{\rm l}} \dot{\upsilon }_{{\rm l}} + p_{{\rm g}} \dot{\upsilon }_{{\rm g}} - H_{{\rm s}} \dot{T},$$

(86)

where υ_l = JφS_l and υ_g = JφS_g represent the pore liquid and gas volumes per unit referential volume, respectively, H_s = Jη^s denotes the entropy density of the solid matrix in the reference configuration.

To simplify, the following potential can be easily defined based on Eq. (86) as:

$$W = \left( {\Psi - J\varphi \psi_{{\rm pore}} } \right) - p_{{\rm l}} \upsilon_{{\rm l}} - p_{{\rm g}} \upsilon_{{\rm g}} .$$

(87)

Next, substituting Eq. (86) into the time derivative of W, we obtain

$$\dot{W}\left( {{\varvec{\varepsilon}},p_{{\rm l}} ,p_{{\rm g}} ,T} \right) = {\rm tr}\left( {{{\varvec{\sigma}} \dot{\varvec{\varepsilon }}}} \right) - \upsilon_{{\rm l}} \dot{p}_{{\rm l}} - \upsilon_{{\rm g}} \dot{p}_{{\rm g}} - H_{{\rm s}} \dot{T}.$$

(88)

Equation (88) indicates that W is a function of state variables ε, p_l, p_g, and T, and thus, the expressions for their conjugate thermodynamics state variables σ, υ_l, υ_g and H_s can be obtained.

Appendix 2

To solve system variables u, p_w, S_w, and T, the governing equations (Eqs. (34), (40), (41) and (45)) needs to be integrated and simplified as follows:

$$G\nabla^{2} {\varvec{u}} + \left( {\frac{G}{1 - 2v}} \right)\nabla \left( {\nabla \cdot {\varvec{u}}} \right) - \alpha \nabla p_{{\rm w}} - \alpha \left( {1 - S_{{\rm w}} } \right)\nabla p_{{\rm c}} - K\beta_{{\rm s}} \nabla T = 0$$

(89)

$$\begin{aligned} & \left( {\frac{{\varphi S_{{\rm w}} }}{{K_{{\rm w}} }} + \frac{{\varphi S_{{\rm g}} }}{{K_{{\rm g}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}} \right)\frac{{\partial p_{{\rm w}} }}{\partial t} + \left( {\frac{{\varphi S_{{\rm g}} }}{{K_{{\rm g}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm g}} } \right)\frac{{\partial p_{{\rm c}} }}{\partial t} + \alpha \frac{{\partial \varepsilon_{v} }}{\partial t} + \left[ { - \left( {\alpha - \varphi } \right)\beta_{{\rm s}} - \varphi S_{{\rm w}} \beta_{{\rm w}} - \varphi S_{{\rm g}} \beta_{{\rm g}} } \right]\frac{\partial T}{{\partial t}} \\ & \quad + \nabla \cdot \left[ { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - \frac{{kk_{{\rm rg}} }}{{\mu_{{\rm g}} }}\nabla \left( {p_{{\rm w}} + p_{{\rm c}} } \right) - \left( {k_{{\rm Tw}} + k_{{\rm Tg}} } \right)\nabla T} \right] = 0 \\ \end{aligned}$$

(90)

$$\begin{aligned} & \left( {\frac{{\varphi S_{{\rm w}} }}{{K_{{\rm w}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm w}} } \right)\frac{{\partial p_{{\rm w}} }}{\partial t} + \left( {\frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm w}} S_{{\rm g}} } \right)\frac{{\partial p_{{\rm c}} }}{\partial t} + \varphi \frac{{\partial S_{{\rm w}} }}{\partial t} + \alpha S_{{\rm w}} \frac{{\partial \varepsilon_{v} }}{\partial t} \\ & \quad + \left[ { - \left( {\alpha - \varphi } \right)S_{{\rm w}} \beta_{{\rm s}} - \varphi S_{{\rm w}} \beta_{{\rm w}} } \right]\frac{\partial T}{{\partial t}} + \nabla \cdot \left( { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - k_{{\rm Tw}} \nabla T} \right) = 0 \\ \end{aligned}$$

(91)

$$\begin{aligned} & \left[ {\left( {1 - \varphi } \right)\rho_{{\rm s}} C_{{\rm s}} + \varphi S_{{\rm w}} \rho_{{\rm w}} C_{{\rm w}} + \varphi S_{{\rm g}} \rho_{{\rm g}} C_{{\rm g}} } \right]\frac{\partial T}{{\partial t}} \\ & \quad + \nabla \cdot \left[ { - \lambda_{{\rm eff}} \nabla T + \rho_{{\rm w}} C_{{\rm w}} T\left( { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - k_{{\rm Tw}} \nabla T} \right) + \rho_{{\rm g}} C_{{\rm g}} T\left( { - \frac{{kk_{{\rm rg}} }}{{\mu_{{\rm g}} }}\nabla p_{{\rm g}} - k_{{\rm Tg}} \nabla T} \right)} \right] = 0. \\ \end{aligned}$$

(92)

Thus, the weak form of above governing equations then can be derived by multiplying them by correspond test functions based on divergence and Gauss’s theorem as follows:

$$\begin{aligned} & \int_{\Omega } {\nabla \tilde{u} \cdot \left[ {G\nabla {{\varvec{u}} + }\left( {\frac{G}{1 - 2v}} \right)\nabla \cdot {\varvec{u}} - \alpha p_{{\rm w}} - \alpha \left( {1 - S_{{\rm w}} } \right)p_{{\rm c}} - K\beta_{{\rm s}} T} \right]} {\rm d}\Omega \\ & \quad - \int_{\Gamma } {\tilde{u} \cdot \left[ {G\nabla {{\varvec{u}} + }\left( {\frac{G}{1 - 2v}} \right)\nabla \cdot {\varvec{u}} - \alpha p_{{\rm w}} - \alpha \left( {1 - S_{{\rm w}} } \right)p_{{\rm c}} - K\beta_{{\rm s}} T} \right]} \cdot {\mathbf{n }}{\rm d}\Gamma = 0 \\ \end{aligned}$$

(93)

$$\begin{aligned} & \int_{\Omega } {\tilde{p}_{{\rm w}} \cdot \left( {\frac{{\varphi S_{{\rm w}} }}{{K_{{\rm w}} }} + \frac{{\varphi S_{{\rm g}} }}{{K_{{\rm g}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}} \right)} \frac{{\partial p_{{\rm w}} }}{\partial t}{\rm d}\Omega + \int_{\Omega } {\tilde{p}_{{\rm w}} \cdot } \left[ {\left( {\frac{{\varphi S_{{\rm g}} }}{{K_{{\rm g}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm g}} } \right)\frac{{\partial p_{{\rm c}} }}{{\partial S_{{\rm w}} }}} \right]\frac{{\partial S_{{\rm w}} }}{\partial t}{\rm d}\Omega \\ & \quad + \int_{\Omega } {\tilde{p}_{{\rm w}} \cdot \left( \alpha \right)} \frac{{\partial \varepsilon_{v} }}{\partial t}{\rm d}\Omega + \int_{\Omega } {\tilde{p}_{{\rm w}} \cdot } \left( { - \left( {\alpha - \varphi } \right)\beta_{{\rm s}} - \varphi S_{{\rm w}} \beta_{{\rm w}} - \varphi S_{{\rm g}} \beta_{{\rm g}} } \right)\frac{\partial T}{{\partial t}} {\rm d}\Omega \\ & \quad - \int_{\Omega } {\nabla \tilde{p}_{{\rm w}} \cdot \left[ { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - \frac{{kk_{{\rm rg}} }}{{\mu_{{\rm g}} }}\nabla \left( {p_{{\rm w}} + p_{{\rm c}} } \right) - \left( {k_{{\rm Tw}} + k_{{\rm Tg}} } \right)\nabla T} \right]} {\rm d}\Omega \\ & \quad + \int_{\Gamma } {\tilde{p}_{{\rm w}} \cdot \left[ { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - \frac{{kk_{{\rm rg}} }}{{\mu_{{\rm g}} }}\nabla \left( {p_{{\rm w}} + p_{{\rm c}} } \right) - \left( {k_{{\rm Tw}} + k_{{\rm Tg}} } \right)\nabla T} \right]} \cdot {\mathbf{n}} {\rm d}\Gamma = 0 \\ \end{aligned}$$

(94)

$$\begin{aligned} & \int_{\Omega } {\tilde{S}_{{\rm w}} \cdot \left( {\frac{{\varphi S_{{\rm w}} }}{{K_{{\rm w}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm w}} } \right)} \frac{{\partial p_{{\rm w}} }}{\partial t}{\rm d}\Omega + \int_{\Omega } {\tilde{S}_{{\rm w}} \cdot } \left[ {\left( {\frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm w}} S_{{\rm g}} } \right)\frac{{\partial p_{{\rm w}} }}{{\partial S_{{\rm w}} }} + \varphi } \right]\frac{{\partial S_{{\rm w}} }}{\partial t}{\rm d}\Omega \\ & \quad + \int_{\Omega } {\tilde{S}_{{\rm w}} \cdot \left( {\alpha S_{{\rm w}} } \right)} \frac{{\partial \varepsilon_{v} }}{\partial t}{\rm d}\Omega + \int_{\Omega } {\tilde{S}_{{\rm w}} \cdot } \left( { - \left( {\alpha - \varphi } \right)S_{{\rm w}} \beta_{{\rm s}} - \varphi S_{{\rm w}} \beta_{{\rm w}} } \right)\frac{\partial T}{{\partial t}} {\rm d}\Omega \\ & \quad - \int_{\Omega } {\nabla \tilde{S}_{{\rm w}} \cdot \left( { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - k_{{\rm Tw}} \nabla T} \right)} {\rm d}\Omega + \int_{\Gamma } {\tilde{S}_{{\rm w}} \cdot \left( { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - k_{{\rm Tw}} \nabla T} \right)} \cdot {\mathbf{n}} {\rm d}\Gamma = 0 \\ \end{aligned}$$

(95)

$$\begin{aligned} & \int_{\Omega } {\tilde{T} \cdot } \left[ {\left( {1 - \varphi } \right)\rho_{{\rm s}} C_{{\rm s}} + \varphi S_{{\rm w}} \rho_{{\rm w}} C_{{\rm w}} + \varphi S_{{\rm g}} \rho_{{\rm g}} C_{{\rm g}} } \right]\frac{\partial T}{{\partial t}}{\rm d}\Omega \\ & \quad - \int_{\Omega } {\nabla \tilde{T} \cdot \left[ { - \lambda_{{\rm eff}} \nabla T + \rho_{{\rm w}} C_{{\rm w}} T\left( { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - k_{{\rm Tw}} \nabla T} \right) + \rho_{{\rm g}} C_{{\rm g}} T\left( { - \frac{{kk_{{\rm rg}} }}{{\mu_{{\rm g}} }}\nabla p_{{\rm g}} - k_{{\rm Tg}} \nabla T} \right)} \right]} {\rm d}\Omega \\ & \quad + \int_{\Gamma } {\tilde{T} \cdot \left[ { - \lambda_{{\rm eff}} \nabla T + \rho_{{\rm w}} C_{{\rm w}} T\left( { - \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}\nabla p_{{\rm w}} - k_{{\rm Tw}} \nabla T} \right) + \rho_{{\rm g}} C_{{\rm g}} T\left( { - \frac{{kk_{{\rm rg}} }}{{\mu_{{\rm g}} }}\nabla p_{{\rm g}} - k_{{\rm Tg}} \nabla T} \right)} \right]} \cdot {\mathbf{n}} {\rm d}\Gamma = 0, \\ \end{aligned}$$

(96)

where $\widetilde{\text{u}}$, ${\widetilde{\text{p}}}_{\text{w}}$, ${\widetilde{\text{S}}}_{\text{w}}$ and $\widetilde{\text{T}}$ are the test functions.

To simplify Eqs. (94)–(96), the following intermediate variables are defined, including:

$$\zeta_{1} = \frac{{\varphi S_{{\rm w}} }}{{K_{{\rm w}} }} + \frac{{\varphi S_{{\rm g}} }}{{K_{{\rm g}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}$$

(97)

$$\zeta_{2} = \left( {\frac{{\varphi S_{{\rm g}} }}{{K_{{\rm g}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm g}} } \right)\frac{{\partial p_{{\rm w}} }}{{\partial S_{{\rm w}} }}$$

(98)

$$\zeta_{3} = - \left( {\alpha - \varphi } \right)\beta_{{\rm s}} - \varphi S_{{\rm w}} \beta_{{\rm w}} - \varphi S_{{\rm g}} \beta_{{\rm g}}$$

(99)

$$\chi_{1} = \frac{{kk_{rw} }}{{\mu_{{\rm w}} }}$$

(100)

$$\chi_{2} = \frac{{kk_{{\rm rg}} }}{{\mu_{{\rm g}} }}$$

(101)

$$\chi_{3} = k_{{\rm Tw}}$$

(102)

$$\chi_{4} = k_{{\rm Tg}}$$

(103)

$$\gamma_{1} = \frac{{\varphi S_{{\rm w}} }}{{K_{{\rm w}} }} + \frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm w}}$$

(104)

$$\gamma_{2} = \left( {\frac{\alpha - \varphi }{{K_{{\rm s}} }}S_{{\rm w}} S_{{\rm g}} } \right)\frac{{\partial p_{{\rm w}} }}{{\partial S_{{\rm w}} }} + \varphi$$

(105)

$$\gamma_{3} = - \left( {\alpha - \varphi } \right)S_{{\rm w}} \beta_{{\rm s}} - \varphi S_{{\rm w}} \beta_{{\rm w}}$$

(106)

$$\varsigma_{1} = \left( {1 - \varphi } \right)\rho_{{\rm s}} C_{{\rm s}} + \varphi S_{{\rm w}} \rho_{{\rm w}} C_{{\rm w}} + \varphi S_{{\rm g}} \rho_{{\rm g}} C_{{\rm g}}$$

(107)

$$\varsigma_{2} = \lambda_{{\rm eff}}$$

(108)

$$\varsigma_{3} = \rho_{{\rm w}} C_{{\rm w}} T$$

(109)

$$\varsigma_{4} = \rho_{{\rm g}} C_{{\rm g}} T$$

(110)

$$\omega = \frac{{\partial p_{{\rm w}} }}{{\partial S_{{\rm w}} }}.$$

(111)

Moreover, by employing the finite-element interpolation function, the field variables to be solved, including displacement (u), water phase pressure (p_w), water saturation (S_w), and temperature (T) can be rewritten in FEM formulation as

$${\varvec{u}} = {\mathbf{N}}_{{\rm u}} \cdot {\mathbf{U}}$$

(112)

$$p_{{\rm w}} = {\mathbf{N}}_{{\rm p}} \cdot {\mathbf{P}}$$

(113)

$$S_{{\rm w}} = {\mathbf{N}}_{{\rm s}} \cdot {\mathbf{S}}$$

(114)

$$T = {\mathbf{N}}_{{\rm T}} \cdot {\mathbf{T}}$$

(115)

$${\varvec{\varepsilon}} = {\mathbf{B}}_{{\rm u}} \cdot {\mathbf{U}}$$

(116)

$$\nabla p_{{\rm w}} = {\mathbf{B}}_{{\rm p}} \cdot {\mathbf{P}}$$

(117)

$$\nabla S_{{\rm w}} = {\mathbf{B}}_{{\rm s}} \cdot {\mathbf{S}}$$

(118)

$$\nabla T = {\mathbf{B}}_{{\rm T}} \cdot {\mathbf{T,}}$$

(119)

where N_u, N_p, N_s, and N_T are the shape functions of these field variables; U, P, S, and T are the vectors of variables at each node.

Thus, Eqs. (93)–(96) can be rewritten as

$$\begin{aligned} & \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm u}}^{{\rm T}} {\mathbf{DB}}_{{\rm u}} {\rm d}\Omega } \right)} {\mathbf{ U}} - \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm u}}^{{\rm T}} \alpha {\mathbf{N}}_{{\rm p}} {\rm d}\Omega } \right)} {\mathbf{ P}} - \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm u}}^{{\rm T}} \alpha \left( {1 - S_{{\rm w}} } \right)f\left( {S_{{\rm w}} } \right){\mathbf{N}}_{{\rm s}} {\rm d}\Omega } \right]} {\mathbf{ S}} \\ & \quad - \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm u}}^{{\rm T}} K\beta_{{\rm s}} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)} {\mathbf{ T}} = \int_{\Gamma } {{\mathbf{N}}_{{\rm u}}^{{\rm T}} \overline{t}_{n} \cdot {\mathbf{n}} {\rm d}\Gamma } \\ \end{aligned}$$

(120)

$$\begin{aligned} & \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \zeta_{1} {\mathbf{N}}_{{\rm p}} {\rm d}\Omega } \right)} {\dot{\mathbf{P}}} + \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \zeta_{2} {\mathbf{N}}_{{\rm s}} {\rm d}\Omega } \right)} {\dot{\mathbf{S}}} + \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \alpha {\mathbf{B}}_{{\rm u}} {\rm d}\Omega } \right)} {\dot{\mathbf{U}}} + \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \zeta_{3} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)} {\dot{\mathbf{T}}} \\ & \quad + \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm p}}^{{\rm T}} \left( {\chi_{1} + \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\rm d}\Omega } \right]} {\mathbf{ P}} + \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm p}}^{{\rm T}} \chi_{2} \omega {\mathbf{B}}_{{\rm s}} {\rm d}\Omega } \right)} {\mathbf{ S}} + \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm p}}^{{\rm T}} \left( {\chi_{3} + \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\rm d}\Omega } \right]} {\mathbf{ T}} \\ & \quad - \int_{\Gamma } {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \left( {\chi_{1} + \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\mathbf{P}} \cdot {\mathbf{n}} {\rm d}\Gamma } - \int_{\Gamma } {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \chi_{2} \omega {\mathbf{B}}_{{\rm s}} {\mathbf{S}}} \cdot {\mathbf{n}} {\rm d}\Gamma - \int_{\Gamma } {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \left( {\chi_{3} + \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\mathbf{T}}} \cdot {\mathbf{n}} {\rm d}\Gamma = 0 \\ \end{aligned}$$

(121)

$$\begin{aligned} & \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \gamma_{1} {\mathbf{N}}_{{\rm p}} {\rm d}\Omega } \right)} {\dot{\mathbf{P}}} + \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \gamma_{2} {\mathbf{N}}_{{\rm s}} {\rm d}\Omega } \right)} {\dot{\mathbf{S}}} + \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \alpha {\mathbf{N}}_{{\rm s}} {\mathbf{SB}}_{{\rm u}} {\rm d}\Omega } \right)} {\dot{\mathbf{U}}} + \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \gamma_{3} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)} {\dot{\mathbf{T}}} \\ & \quad + \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm s}}^{{\rm T}} \chi_{1} {\mathbf{B}}_{{\rm p}} {\rm d}\Omega } \right)} {\mathbf{ P}} + \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm s}}^{{\rm T}} \chi_{3} {\mathbf{B}}_{{\rm T}} {\rm d}\Omega } \right)} {\mathbf{ T}} - \int_{\Gamma } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \chi_{1} {\mathbf{B}}_{{\rm p}} \cdot {\mathbf{P}}} \right)} \cdot {\mathbf{n}} {\rm d}\Gamma - \int_{\Gamma } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \chi_{3} {\mathbf{B}}_{{\rm T}} \cdot T} \right)} \cdot {\mathbf{n}} {\rm d}\Gamma = 0 \\ \end{aligned}$$

(122)

$$\begin{aligned} & \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \varsigma_{1} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)} {\dot{\mathbf{T}}} + \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{2} + \varsigma_{3} \chi_{3} + \varsigma_{4} \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\rm d}\Omega } \right]} {\mathbf{ T}} + \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{3} \chi_{1} + \varsigma_{4} \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\rm d}\Omega } \right]} {\mathbf{ P}} \\ & \quad + \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm T}}^{{\rm T}} \varsigma_{4} \chi_{2} \omega {\mathbf{B}}_{S} {\rm d}\Omega } \right)} {\mathbf{ S}} - \int_{\Gamma } {\left[ {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{2} + \varsigma_{3} \chi_{3} + \varsigma_{4} \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\mathbf{T}}} \right]} \cdot {\mathbf{n}} {\rm d}\Gamma \\ & \quad - \int_{\Gamma } {\left[ {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{3} \chi_{1} + \varsigma_{4} \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\mathbf{P}}} \right]} \cdot {\mathbf{n}} {\rm d}\Gamma - \int_{\Gamma } {\left( {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \varsigma_{4} \chi_{2} \omega {\mathbf{B}}_{S} {\mathbf{S}}} \right)} \cdot {\mathbf{n}} {\rm d}\Gamma = 0. \\ \end{aligned}$$

(123)

Then, the following variables are further introduced for simplification:

$${\mathbf{F}}_{{\rm u}} = \int_{\Gamma } {{\mathbf{N}}_{{\rm u}}^{{\rm T}} \overline{t}_{n} \cdot {\mathbf{n}} {\rm d}\Gamma }$$

(124)

$${\mathbf{F}}_{{\rm p}} = \int_{\Gamma } {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \left( {\chi_{1} + \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\mathbf{P}} \cdot {\mathbf{n}} {\rm d}\Gamma } + \int_{\Gamma } {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \chi_{2} \omega {\mathbf{B}}_{{\rm s}} {\mathbf{S}}} \cdot {\mathbf{n}} {\rm d}\Gamma + \int_{\Gamma } {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \left( {\chi_{3} + \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\mathbf{T}}} \cdot {\mathbf{n}} {\rm d}\Gamma$$

(125)

$${\mathbf{F}}_{{\rm s}} = \int_{\Gamma } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \chi_{1} {\mathbf{B}}_{{\rm p}} \cdot {\mathbf{P}}} \right)} \cdot {\mathbf{n}} {\rm d}\Gamma + \int_{\Gamma } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \chi_{3} {\mathbf{B}}_{{\rm T}} \cdot T} \right)} \cdot {\mathbf{n}} {\rm d}\Gamma$$

(126)

$$\begin{aligned} {\mathbf{F}}_{{\rm T}} & = \int_{\Gamma } {\left[ {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{2} + \varsigma_{3} \chi_{3} + \varsigma_{4} \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\mathbf{T}}} \right]} \cdot {\mathbf{n}} {\rm d}\Gamma \\ & \quad + \int_{\Gamma } {\left[ {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{3} \chi_{1} + \varsigma_{4} \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\mathbf{P}}} \right]} \cdot {\mathbf{n}} {\rm d}\Gamma + \int_{\Gamma } {\left( {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \varsigma_{4} \chi_{2} \omega {\mathbf{B}}_{S} {\mathbf{S}}} \right)} \cdot {\mathbf{n}} {\text{d}}\Gamma \\ \end{aligned}$$

(127)

$${\mathbf{M}}_{1} = \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm u}}^{{\rm T}} {\mathbf{DB}}_{{\rm u}} {\rm d}\Omega } \right)}$$

(128)

$${\mathbf{Q}}_{1} = \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm u}}^{{\rm T}} \alpha {\mathbf{N}}_{{\rm p}} {\rm d}\Omega } \right)}$$

(129)

$${\mathbf{W}}_{1} = \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm u}}^{{\rm T}} \alpha \left( {1 - S_{{\rm w}} } \right)f\left( {S_{{\rm w}} } \right){\mathbf{N}}_{{\rm s}} {\rm d}\Omega } \right]}$$

(130)

$${\mathbf{H}}_{1} = \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm u}}^{{\rm T}} K\beta_{{\rm s}} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)}$$

(131)

$${\mathbf{M}}_{2t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \alpha {\mathbf{B}}_{{\rm u}} {\rm d}\Omega } \right)}$$

(132)

$${\mathbf{Q}}_{2t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \zeta_{1} {\mathbf{N}}_{{\rm p}} {\rm d}\Omega } \right)}$$

(133)

$${\mathbf{W}}_{2t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \zeta_{2} {\mathbf{N}}_{{\rm s}} {\rm d}\Omega } \right)}$$

(134)

$${\mathbf{H}}_{2t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm p}}^{{\rm T}} \zeta_{3} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)}$$

(135)

$${\mathbf{Q}}_{2} = \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm p}}^{{\rm T}} \left( {\chi_{1} + \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\rm d}\Omega } \right]}$$

(136)

$${\mathbf{W}}_{2} = \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm p}}^{{\rm T}} \chi_{2} \omega {\mathbf{B}}_{{\rm s}} {\rm d}\Omega } \right)}$$

(137)

$${\mathbf{H}}_{2} = \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm p}}^{{\rm T}} \left( {\chi_{3} + \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\rm d}\Omega } \right]}$$

(138)

$${\mathbf{M}}_{3t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \alpha {\mathbf{N}}_{{\rm s}} {\mathbf{SB}}_{{\rm u}} {\rm d}\Omega } \right)}$$

(139)

$${\mathbf{Q}}_{3t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \gamma_{1} {\mathbf{N}}_{{\rm p}} {\rm d}\Omega } \right)}$$

(140)

$${\mathbf{W}}_{3t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \gamma_{2} {\mathbf{N}}_{{\rm s}} {\rm d}\Omega } \right)}$$

(141)

$${\mathbf{H}}_{3t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm s}}^{{\rm T}} \gamma_{3} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)}$$

(142)

$${\mathbf{Q}}_{3} = \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm s}}^{{\rm T}} \chi_{1} {\mathbf{B}}_{{\rm p}} {\rm d}\Omega } \right)}$$

(143)

$${\mathbf{H}}_{3} = \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm s}}^{{\rm T}} \chi_{3} {\mathbf{B}}_{{\rm T}} {\rm d}\Omega } \right)}$$

(144)

$${\mathbf{H}}_{3t} = \int_{\Omega } {\left( {{\mathbf{N}}_{{\rm T}}^{{\rm T}} \varsigma_{1} {\mathbf{N}}_{{\rm T}} {\rm d}\Omega } \right)}$$

(145)

$${\mathbf{Q}}_{3} = \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{3} \chi_{1} + \varsigma_{4} \chi_{2} } \right){\mathbf{B}}_{{\rm p}} {\rm d}\Omega } \right]}$$

(146)

$${\mathbf{W}}_{3} = \int_{\Omega } {\left( {{\mathbf{B}}_{{\rm T}}^{{\rm T}} \varsigma_{4} \chi_{2} \omega {\mathbf{B}}_{S} {\rm d}\Omega } \right)}$$

(147)

$${\mathbf{H}}_{3} = \int_{\Omega } {\left[ {{\mathbf{B}}_{{\rm T}}^{{\rm T}} \left( {\varsigma_{2} + \varsigma_{3} \chi_{3} + \varsigma_{4} \chi_{4} } \right){\mathbf{B}}_{{\rm T}} {\rm d}\Omega } \right]} .$$

(148)

Combining Eqs. (120)–(148), the governing equations (Eqs. (34), (40), (41) and (45)) can then be spatially discretized by FEM, so that a set of nonlinear equations can be written in FEM form as

$$\left\{ \begin{array}{ll} {\mathbf{M}}_{1} \cdot {\mathbf{U}} - {\mathbf{Q}}_{1} \cdot {\mathbf{P}} - {\mathbf{W}}_{1} \cdot {\mathbf{S}} - {\mathbf{H}}_{1} \cdot {\mathbf{T}} = {\mathbf{F}}_{{\rm u}} \hfill \\ {\mathbf{M}}_{2t} \cdot {\dot{\mathbf{U}}} + {\mathbf{Q}}_{2t} \cdot {\dot{\mathbf{P}}} + {\mathbf{W}}_{2t} \cdot {\dot{\mathbf{S}}} + {\mathbf{H}}_{2t} \cdot {\dot{\mathbf{T}}} + {\mathbf{Q}}_{2} \cdot {\mathbf{P}} + {\mathbf{W}}_{2} \cdot {\mathbf{S}} + {\mathbf{H}}_{2} \cdot {\mathbf{T}} = {\mathbf{F}}_{{\rm p}} \hfill \\ {\mathbf{M}}_{3t} \cdot {\dot{\mathbf{U}}} + {\mathbf{Q}}_{3t} \cdot {\dot{\mathbf{P}}} + {\mathbf{W}}_{3t} \cdot {\dot{\mathbf{S}}} + {\mathbf{H}}_{3t} \cdot {\dot{\mathbf{T}}} + {\mathbf{Q}}_{3} \cdot {\mathbf{P}} + {\mathbf{H}}_{3} \cdot {\mathbf{T}} = {\mathbf{F}}_{{\rm s}} \hfill \\ {\mathbf{H}}_{4t} \cdot {\dot{\mathbf{T}}} + {\mathbf{Q}}_{4} \cdot {\mathbf{P}} + {\mathbf{W}}_{4} \cdot {\mathbf{S}} + {\mathbf{H}}_{4} \cdot {\mathbf{T}} = {\mathbf{F}}_{{\rm T}} \hfill \\ \end{array} \right..$$

(149)

Then, using the backward difference method for time-domain discretization, the partial derivative with time of the variables in Eqs. (112)–(115) can be expressed as

$${\dot{\mathbf{U}}} = \frac{{{\mathbf{U}} - {\mathbf{U}}_{{\left( {t_{n} - 1} \right)}} }}{\Delta t}$$

(150)

$${\dot{\mathbf{P}}} = \frac{{{\mathbf{P}} - {\mathbf{P}}_{{\left( {t_{n} - 1} \right)}} }}{\Delta t}$$

(151)

$${\dot{\mathbf{S}}} = \frac{{{\mathbf{S}} - {\mathbf{S}}_{{\left( {t_{n} - 1} \right)}} }}{\Delta t}$$

(152)

$${\dot{\mathbf{T}}} = \frac{{{\mathbf{T}} - {\mathbf{T}}_{{\left( {t_{n} - 1} \right)}} }}{\Delta t},$$

(153)

where t_n-1 is the previous time step.

Finally, substituting Eqs. (150)–(153) into Eq. (149), the FEM formulation of the unsaturated THM coupling model can be written in a compact form as

$$\left[ {\begin{array}{*{20}c} {{\mathbf{M}}_{1} } & { - {\mathbf{Q}}_{1} } & { - {\mathbf{W}}_{1} } & { - {\mathbf{H}}_{1} } \\ {{\mathbf{M}}_{2t} } & {{\mathbf{Q}}_{2t} + \Delta t{\mathbf{Q}}_{2} } & {{\mathbf{W}}_{2t} + \Delta t{\mathbf{W}}_{2} } & {{\mathbf{H}}_{2t} + \Delta t{\mathbf{H}}_{2} } \\ {{\mathbf{M}}_{3t} } & {{\mathbf{Q}}_{3t} + \Delta t{\mathbf{Q}}_{3} } & {{\mathbf{W}}_{3t} } & {{\mathbf{H}}_{3t} + \Delta t{\mathbf{H}}_{3} } \\ 0 & {\Delta t{\mathbf{Q}}_{4} } & {\Delta t{\mathbf{W}}_{4} } & {{\mathbf{H}}_{4t} + \Delta t{\mathbf{H}}_{4} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\mathbf{U}} \\ {\mathbf{P}} \\ {\mathbf{S}} \\ {\mathbf{T}} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {{\mathbf{F}}_{{\rm u}} } \\ {\Delta t{\mathbf{F}}_{{\rm p}} + {\mathbf{M}}_{2t} {\mathbf{U}}_{{\left( {t_{n} - 1} \right)}} + {\mathbf{Q}}_{2t} {\mathbf{P}}_{{\left( {t_{n} - 1} \right)}} + {\mathbf{W}}_{2t} {\mathbf{S}}_{{\left( {t_{n} - 1} \right)}} + {\mathbf{H}}_{2t} {\mathbf{T}}_{{\left( {t_{n} - 1} \right)}} } \\ {\Delta t{\mathbf{F}}_{{\rm s}} + {\mathbf{M}}_{3t} {\mathbf{U}}_{{\left( {t_{n} - 1} \right)}} + {\mathbf{Q}}_{3t} {\mathbf{P}}_{{\left( {t_{n} - 1} \right)}} + {\mathbf{W}}_{3t} {\mathbf{S}}_{{\left( {t_{n} - 1} \right)}} + {\mathbf{H}}_{3t} {\mathbf{T}}_{{\left( {t_{n} - 1} \right)}} } \\ {\Delta t{\mathbf{F}}_{{\rm T}} + + {\mathbf{H}}_{4t} {\mathbf{T}}_{{\left( {t_{n} - 1} \right)}} } \\ \end{array} } \right].$$

(154)