## 2 Mathematical Formulations
### Governing Equations  
The volume fraction (\(\alpha\)) represents the composition of the control cell. Therein, \(\alpha_{i}=1\) suggests the cell is fully occupied by the \(i\) phase, \(\alpha_{i}=0\) indicates there are no \(i\) phase in the cell and \(0<\alpha_{i}<1\) means that only part of the cell contains \(i\) phase. The heat transfer process during cryogenic fracturing includes three phases which are the water, ice, and cryogenic fluid, respectively. For the tracking of interface (water-cryogenic fluid, water-ice), the VOF method in Euler-Euler multiphase system is adopted, in which the interface is tracked by solving volume fraction equations. In the present work, the volume fraction of the \(i^{\text{th}}\) phase can be calculated as [43; 44]:  
\[\frac{\partial\alpha_{i}}{\partial t}+\
abla\cdot(\alpha_{i}\overset{\rightarrow }{u}_{i})=S_{\alpha_{i}} \tag{1}\]  
where \(S_{\alpha i}\) is the source term of the ice phase, which is induced by the mass transfer during ice formation, and the details would be discussed in the following section. In the computational domain with \(n\) phases, the volume fractions of these phases can be referred to the following constraint:  
\[\sum\
olimits_{i=1}^{n}\alpha_{i}=1\;\;. \tag{2}\]  
Therein, the fluid density in the cell can be expressed according to the volume fractions of the \(n\) phases:  
\[\rho=\sum\alpha_{i}\rho_{i}\;. \tag{3}\]  
The mass conservation equations of each phase are solved, and the continuity equation of the \(i^{\text{th}}\) phase could be written as:  
\[\frac{\partial\alpha_{i}\rho_{i}}{\partial t}+\
abla\cdot(\alpha_{i}\rho_{i} \overset{\rightarrow}{u}_{i})=S_{m_{i}} \tag{4}\]  
where \(\rho_{i}\) and \(\overset{\rightarrow}{u}_{i}\) represent the fluid density and the fluid velocity of the \(i^{\text{th}}\) phase, respectively. Furthermore, \(S_{mi}\) represents the mass source which equals to zero at the flow interface (water-cryogenic fluid) and equals to mass transfer rate at the water-ice interface. The formulation of the mass source is presented in Section 2.2. In this case, the momentum equations are applied to model the flow behavior of various phases, which can be written as following for the \(i^{\text{th}}\) phase. \[\frac{\partial\rho_{i}\alpha_{i}\overset{\rightarrow}{u}_{i}}{\partial t}+ \
abla\cdot(\rho_{i}\alpha_{i}\overset{\rightarrow}{u}_{i}\overset{ \rightarrow}{u}_{i})=-\alpha_{i}\
abla p+\
abla\cdot\overset{\rightarrow}{ \tau}_{i}+K_{ji}(\overset{\rightarrow}{u}_{j}-\overset{\rightarrow}{u}_{i})+S _{u_{i}} \tag{5}\]  
where \(p\) is pressure, \(K_{ji}\) is interfacial drag coefficient, \(S_{ui}\) is the source term of the momentum induced by the mass transfer which would be given in the next section in details. Furthermore, \(\overset{=}{\tau}\) is the tensor of the stress-strain and it is described as:  
\[\overset{=}{\tau}_{i}=\mu_{i}\alpha_{i}(\
abla\overset{\rightarrow}{u}_{i}+ \
abla\overset{\rightarrow}{u}_{i}^{T})+\alpha_{i}\bigg{(}\omega_{i}-\frac{2}{ 3}\mu_{i}\bigg{)}\
abla\cdot\overset{\rightarrow}{u}_{i}^{\overset{=}{T}} \tag{6}\]where \(\overset{=}{I}\) is unit tensor. In all the cases described in this work, Reynolds number is not allowed to exceed the critical value as about 2000, suggesting that a laminar flow model is used for the simulation. The energy conservation equation of the \(i^{\text{th}}\) phase is expressed as:  
\[\frac{\partial\rho_{i}\alpha_{i}H_{i}}{\partial t}+\
abla\cdot(\rho_{i}\alpha_{ i}\overset{\rightarrow}{u_{i}}H_{i})=\
abla\cdot(k_{i}\
abla T_{i})+S_{H_{i}} \tag{7}\]  
where \(H\) represents specific enthalpy, \(k\) represents the thermal conductivity, \(T\) represents the fluid temperature, and \(S_{Hi}\) represents the heat source due to mass transfer which would be described in the next section. As a result, the governing equations mentioned above constitute a time-dependent and space-dependent model for the physical problems of the temperature distribution, the fluid flow and the phase change of solid-liquid.