## 2 Governing Equations for Fractured Unsaturated Porous Media under the Non-Isothermal Condition
### Description of Strong Discontinuities and Solution Strategy  
Spatial discretization is performed using the extended finite element method (XFEM) that enriches ordinary shape functions with additional terms for describing strong discontinuities. Temporal discretization is performed using the finite difference method. In the following parts of this section, each type of discretization is described briefly. To begin with, spatial discretization should be performed to satisfy the physical condition of the proposed problem. For instance, the liquid phase should be entered into the crack, and the type of function for spatial discretization plays an important role in this action. These requirements include (1) the strong discontinuity of displacement field for fracture opening and (2) the weak discontinuity of water pressure and temperature fields for moisture and heat exchanges. The strong discontinuity for crack opening means the displacement field and the gradient of displacement, strain, which has to be discontinuous between two faces of the fracture. The weak discontinuity of water pressure and temperature results in discontinuous water pressure gradients and temperature normal to the fracture. To have such fields, the classic finite element should be enriched to ensure that aforementioned assumptions are fulfilled. We note that only the nodes bisected by cracks are enriched. Since the fracture opening is zero around the crack tip and the moisture and heat flows are continuous, the nodes on the side where the crack tip lay on are not enriched. In the XFEM, all the unknown variables are decomposed into standard and enriched parts. To model the desired discontinuity, the enrichment function is embedded in the enriched part. The extended finite element approximation for the displacement field can be given as follows:  
Figure 2: (**a**) Schematics of crack-opening with traction-separation law, (**b**) Demonstration of the boundary condition on the crack and model surface. \[u(x,t)=\sum_{i\in\Omega}N_{i}^{u}(x)u_{i}(t)+\sum_{i\in\Omega_{curr}}N_{i}^{u}(x)( H(x)-H(x_{i}))\bar{u}_{i}(i) \tag{14}\]  
where \(N_{i}^{u}\) is the shape function and \(u_{i}(t)\) and \(\bar{u}_{i}(t)\) are the standard and enriched degrees of displacement, respectively. \(\Omega\) and \(\Omega_{curr}\) are the set of all the nodes and the enriched nodes in the domain. \(H(x)\) denotes the Heaviside function, taken into account for the consideration of the strong discontinuity of displacement field that can be given as follows:  
\[H(x)=\left\{ \begin{array}{rl}+1,&if(x-x_{\Gamma}).n_{\Gamma_{d}}\geq 0\\ &-1,\;otherwise\end{array} \right\} \tag{15}\]  
where \(x_{\Gamma_{d}}\) is the nearest projection of any point on the crack path. To obtain the value of the field variable \(u\) on any enriched point, the Heaviside function is shifted by \((H(x)-H(x_{i}))\) around any node. Similarly, the pressure and temperature unknowns are approximated as follows:  
\[P_{w}(x,t)=\sum_{i\in\Omega}N_{i}^{P_{w}}(x)P_{wi}(t)+\sum_{i\in\Omega_{curr} }N_{i}^{P_{w}}(x)(D(x)-D(x_{i}))R(x)\bar{P}_{wi}(t) \tag{16}\]  
\[T(x,t)=\sum_{i\in\Omega}N_{i}^{T}(x)T_{i}(t)+\sum_{i\in\Omega_{curr}}N_{i}^{T }(x)(D(x)-D(x_{i}))R(x)\bar{T}_{i}(t) \tag{17}\]  
where \(N_{i}^{P_{w}}\) and \(N_{i}^{T}\) are standard shape functions for water pressure and temperature, respectively. \(P_{wi}(t)\) and \(T_{i}(t)\) are standard degrees of freedom and \(\bar{P}_{wi}(t)\) and \(\bar{T}_{i}(t)\) are enriched degrees of freedom for water pressure and temperature, respectively. In order to implement the weak discontinuity in the formulation, the distance function (\(D(x)\)) is used. The weak form of governing equations should also be discretized over time domain [24].