## 2 Governing Equations for Fractured Unsaturated Porous Media under the Non-Isothermal Condition
### Crack-Porous Media Interaction Equations  
Since geological materials such as soil and rock do not exhibit the singularity of stress predicted by the linear elastic fracture mechanics (LEFM) at the crack tip, the cohesive crack model, which develops when the tip stress exceeds the tensile strength, could be a powerful alternative for the quasi-brittle materials. In addition to cohesive traction, the suction developed in the pores of the porous media acts on the crack faces as follows:  
\[t_{d}=t_{d}(\llbracket u\rrbracket_{m})+S_{w}P_{w} \tag{11}\]  
where \(\llbracket u\rrbracket_{m}\) is the modified crack opening due to the temperature difference between crack faces (\(\llbracket u\rrbracket_{m}=\llbracket u\rrbracket(1-\alpha_{s}\llbracket T\rrbracket)\), and \(t_{d}\) is the dependant cohesive force, which is considered as a linear function of \(\llbracket u\rrbracket_{m}\), as in Figure 2a. In this study, the linear function is adopted to reduce the complexity associated with cracking such that we focus more on the calculation of coupled suction and thermal impacts on porous media. Moreover, this treatment is aligned with the nature of cracks in fluid infiltrating porous media, which shows in-between the brittle and ductile fracture behaviors.  
As shown in Figure 2b, the crack in porous media leads to the exchange of moisture (\(\overline{q}_{\text{{\tiny{ind}}}}\)) and heat (\(\overline{q}_{\text{{\tiny{T}}}d}\)) between the crack and the porous media. These exchanged flows are calculated by rewriting the moisture mass and energy conservation equations around the crack as follows:  
\[\begin{array}{c}\overline{q}_{\text{{\tiny{ind}}}}(x,t)=-C_{wud}\Big{(}2h \Big{\langle}\frac{\partial v_{s}}{\partial X^{\prime}}\Big{\rangle}+ \llbracket w\rrbracket_{s}\Big{)}-2h\Big{(}C_{wud}\dot{P}_{w}+C_{wtd}\dot{T} \Big{)}\\ \qquad\qquad\qquad\qquad-2h\frac{\partial}{\partial X^{\prime}}(\rho_{w}(U+V) )\\ \overline{q}_{\text{{\tiny{T}}}d}(x,t)=-\chi_{\text{{\tiny{cd}}}}\Big{(}2h \Big{\langle}\frac{\partial v_{s}}{\partial X^{\prime}}\Big{\rangle}+ \llbracket w\rrbracket_{s}\Big{)}-2h\Big{(}\chi_{\text{{\tiny{rad}}}}\dot{P}_{w} +\chi_{\text{{\tiny{T}}}d}\dot{T}\Big{)}\\ \qquad\qquad\qquad-2h\Big{(}\chi_{\text{{\tiny{rad}}}}\dot{P}_{w}+\chi_{\text{ {\tiny{T}}}d}\dot{T}\Big{)}-2h\frac{\partial}{\partial X^{\prime}}(Q)\end{array} \tag{12}\]  
Here, the crack opening is expressed by \(2h\), and \(v_{s}\) and \(w_{s}\) are the tangential and normal velocities of solid with respect to the crack. We note that \(\langle*\rangle\) and \([[*]]\) demonstrate the average and difference of variables at two opposite sides of the crack. The subscript \(d\) defines the coefficients at the damaged area and there is a difference between damaged coefficients and continuum coefficients for the intrinsic permeability and porosity. The intrinsic permeability is defined using the cubic law, and the porosity is obtained as \(n_{d}=[[u]]/[[u]]_{c}\). We refer more details about other model parameters including \(C_{wud}\), \(C_{wud}\), \(C_{wud}\), \(C_{wud}\), \(\chi_{\text{{\tiny{cd}}}d}\), \(\chi_{\text{{\tiny{T}}}d}\), and \(\chi_{\text{{\tiny{T}}}d}\) to [23].