## 2 Mathematical Representation  
### CFD Simulation  
The present study used the commercial code ANSYS FLUENT (V.17.2), which has been widely used and evaluated, for simulations under the chosen scenarios. We simultaneously solved the Reynolds-averaged Navier-Stokes (RANS) equations and the scalar transport equation of HNS. Considering the large scales of the current velocity and representative length, the Reynolds number becomes much larger than the critical value, representing fully turbulent flows. In addition, we used the well-known standard k-\(\varepsilon\) model. In our previous study, we conducted simulation to evaluate the CFD methods and compared the predictions with measured data obtained in with the lab-scale oil spill experiment [12; 13]. This experiment was conducted by [PERSON] et al. [13] who provided the measured data for oil propagation in the oil tank under different scenarios. The governing equations are expressed as follows:  
\[\frac{\partial\rho}{\partial t}+\frac{\partial(\rho\overline{u}_{i})}{\partial x _{i}}=0 \tag{1}\]  
\[\frac{\partial}{\partial t}(\rho\overline{u}_{i})+\frac{\partial}{\partial x_{i }}\Big{(}\rho\overline{u}_{i}\overline{u}_{j}\Big{)}=-\frac{\partial p}{ \partial x_{i}}+\frac{\partial}{\partial x_{i}}\bigg{[}\mu\bigg{(}\frac{ \partial\overline{u}_{i}}{\partial x_{j}}+\frac{\partial\overline{u}_{j}}{ \partial x_{i}}-\frac{2}{3}\delta_{ij}\frac{\partial\overline{u}_{i}}{ \partial x_{l}}\bigg{)}\bigg{]}+\frac{\partial}{\partial x_{j}}\Big{(}-\rho \overline{u^{\prime}iu^{\prime}}_{i}\Big{)} \tag{2}\]  
\[\frac{\partial}{\partial t}(\rho Y_{k})+\frac{\partial(\rho\overline{u}_{i}Y_ {k})}{\partial x_{i}}=-\frac{\partial f_{jk}}{\partial x_{j}} \tag{3}\]  
where \(\rho\) denotes the density, \(\mu\) is the dynamic viscosity, and \(J_{jk}\) and \(Y_{k}\) are the diffusion flux and mass fraction of a species \(k\), respectively. The equations for the turbulent kinetic energy, \(k\), and rate of energy dissipation, \(\varepsilon\), are expressed as follows:  
\[\frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_{i}}(\rho \overline{u}_{i}k)=\frac{\partial}{\partial x_{j}}\bigg{[}\bigg{(}\mu+\frac{ \mu_{t}}{\sigma_{k}}\bigg{)}\frac{\partial k}{\partial x_{j}}\bigg{]}+G_{k}+G_ {b}-\rho\varepsilon \tag{4}\]  
\[\frac{\partial}{\partial t}(\rho\varepsilon)+\frac{\partial}{\partial x_{i}} (\rho\overline{u}_{i}\varepsilon)=\frac{\partial}{\partial x_{j}}\bigg{[} \bigg{(}\mu+\frac{\mu_{t}}{\sigma_{\varepsilon}}\bigg{)}\frac{\partial \varepsilon}{\partial x_{j}}\bigg{]}+C_{1}\frac{\varepsilon}{k}(G_{k}+G_{b}) -C_{2}\rho\frac{\varepsilon^{2}}{k} \tag{5}\]  
\[\mu_{t}=C_{\mu}\rho\frac{k^{2}}{\varepsilon} \tag{6}\]  
where \(\mu_{t}\) is the turbulent viscosity, \(G_{k}\) is the production of the turbulent kinetic energy due to the mean velocity gradient, \(G_{b}\) is the production of the turbulence kinetic energy due to buoyancy, \(C_{1}\) and \(C_{2}\) are the model constants, and \(\sigma_{k}\) and \(\sigma_{\varepsilon}\) are the turbulent Prandtl numbers of \(k\) and \(\varepsilon\), respectively. The model constants were taken as  
\[C_{1}=1.44,\quad C_{2}=1.92,\quad C_{\mu}=0.09,\quad\sigma_{k}=1.0,\quad \sigma_{\varepsilon}=1.3 \tag{7}\]  
In simulation, we used the first-order upwind scheme for diffusion and transient terms in the governing equations because the Peclet number is relatively high for all cases. We first constructed the computational grid systems shown in Figure 1. The geometry of the ship was basically provided by the Korea Research Institute of Ships and Ocean Engineering (KRISO). Regarding the scenarios of HNS spill accidents, KRISO suggested the key parameters, which should consider and determine specific conditions based on real circumstances. It is assumed that the HNS spills directly through a crack formed on a ship before, inevitably, spreading away from the ship. Empirically, grounding and collision are the main cause of HNS spills when accidents occur [14; 15].