## 2 Mathematical Modeling  
The governing equations of the melting process of the ice in a shell and coil tube ice storage are solved for laminar, unsteady and three-dimensional flow configurations using the enthalpy-porosity method [21]. The viscous dissipation term is considered as negligible [21]. The enthalpy-porosity method has proven to be an efficient method to simulate the melting and solidification of PCMs (like ice) [3; 13; 20]. In order to include the effects of buoyancy on the natural convection in the storage system, the Boussinesq approximation was applied to the problem [22]  
\[\rho=\rho_{0}\Big{[}\beta\Big{(}T-T_{liquid}\Big{)}+1\Big{]}^{-1} \tag{1}\]  
Continuity:  
\[\
abla.\stackrel{{\rightarrow}}{{V}}=0 \tag{2}\]  
Momentum:  
\[\frac{\partial V}{\partial t}+\stackrel{{\rightarrow}}{{V}}. \
abla\stackrel{{\rightarrow}}{{V}}=\frac{1}{\rho}(-\
abla\rho+ \mu\
abla^{2}\stackrel{{\rightarrow}}{{V}}+\rho\beta(T-T_{ref}) )+\stackrel{{\rightarrow}}{{S}} \tag{3}\]  
Energy:  
\[\frac{\partial h_{sens}}{\partial t}+\frac{\partial h_{lat}}{\partial t}+ \
abla.(\stackrel{{\rightarrow}}{{V}}h_{sens})=\
abla.(\frac{k }{\rho c_{p}}\
abla h_{sens}). \tag{4}\]  
The material total enthalpy, which can be computed by the sum of sensible and latent heat is  
\[h_{tot}=h_{sens}+h_{lat} \tag{5}\]  
where  
\[h_{sens}=h_{ref}+\int\limits_{T_{ref}}^{T}c_{p}dT=h_{ref}+c_{p}\int\limits_{T_ {ref}}^{T}dT \tag{6}\]  
\[h_{lat}=\sum\limits_{i=1}^{n}\lambda_{i}L. \tag{7}\]  
in which \(h_{lat}\) can vary from zero (for solid phase) to \(L\) (for the liquid phase) and \(\lambda\) is the liquid fraction. According to the correlation, the total latent heat can be computed by the summation of the latent heat from every cell in the computational domain at each time step. After that, the sensible heat can be obtained from  
\[h_{sens}=h_{tot}-h_{lat}. \tag{8}\]  
and the liquid fraction, \(\lambda\), can be defined as [23; 24]  
\[\lambda=\left\{\begin{array}{l l}\frac{h_{lat}}{h_{lat}}=0&if&T<T_{s}\\ \frac{h_{lat}}{h_{lat}}=1&if&T<T_{liq}\\ \frac{h_{lat}}{L}=\frac{T-T_{s}}{T_{liq}-T_{s}}&if&T_{s}<T<T_{liq}\end{array} \right\}. \tag{9}\]  
[PERSON]'s law damping term is usually added to the momentum equation to consider the influence of the phase change on the convective heat transfer. Generally, the solidus temperature is the highest temperature at which a solid material starts melting. The liquidus is the temperature at which a material is completely melted. The mushy zone corresponds to intermediate temperatures between the solidus and liquidus temperatures. According to [24; 25], the temperature span in the mushy layer remains around 0.1 K. Also, [PERSON] et al. [26] noted that for pure material (like water), the difference between the solidus and liquidus temperatures is below 0.5 K and can be set to zero without introducing a noticeable error when simulating phase change processes.  
The term S in Equation (3) is the source term representing the [PERSON]'s law damping term. It can be defined as  
\[\overrightarrow{S}=-\frac{\left(1-\lambda\right)^{2}}{\lambda^{3}}A_{mush} \overrightarrow{V} \tag{10}\]  
where \(A_{mush}\) is the mushy zone constant. This constant is usually fixed to values within the range \(10^{4}\)-\(10^{7}\). In this study, \(A_{mush}\) is assumed to be constant and its value was set to \(10^{5}\).