## 2 Numerical Methods  
### Governing Equations and Details for Simulation  
The present study solves a set of conservation equations for mass, momentum, and species that are expressed as follows. \[\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u_{i})}{\partial x_{i}}=0, \tag{1}\]  
\[\frac{\partial}{\partial t}(\rho u_{i})+\rho u_{j}\frac{\partial u_{i}}{ \partial x_{j}}=-\frac{\partial p}{\partial x_{i}}+\frac{\partial\tau_{ji}}{ \partial x_{j}}+\rho g, \tag{2}\]  
\[\frac{\partial}{\partial t}(\rho Y_{m})+\frac{\partial(\rho u_{i}Y_{m})}{ \partial x_{i}}=-\frac{\partial J_{jm}}{\partial x_{j}}, \tag{3}\]  
where \(\rho\), \(u\), \(x\), and \(t\) are the density, velocity, position, and time, respectively. \(p\) is the pressure and \(\tau_{ji}\) is the stress, \(J_{jm}\) is the diffusion flux of species \(m\), and \(Y_{m}\) is the mass fraction of species \(m\). The standard \(k\)-\(\epsilon\) turbulence model was used, and the standard wall function was employed to describe the turbulent dissipation energy near the wall. The equations for the turbulent kinetic energy \(k\) and the rate of energy dissipation \(\varepsilon\) are as follows. \[\frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_{i}}(\rho u_{i }k)=\frac{\partial}{\partial x_{j}}\Bigg{[}\left(\mu+\frac{\mu_{t}}{\sigma_{k }}\right)\frac{\partial k}{\partial x_{j}}\Bigg{]}+C_{k}+G_{b}-\rho\varepsilon, \tag{4}\]  
\[\frac{\partial}{\partial t}(\rho\varepsilon)+\frac{\partial}{\partial x_{i}}( \rho u_{i}\varepsilon)=\frac{\partial}{\partial x_{j}}\Bigg{[}\left(\mu+\frac {\mu_{t}}{\sigma_{\varepsilon}}\right)\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}\]  
where \(G_{k}\) is the production of turbulent kinetic energy due to mean velocity gradients and \(G_{b}\) is the production of turbulence kinetic energy due to buoyancy. \(C_{1}\) and \(C_{2}\) are model constants and \(\sigma_{k}\) and \(\sigma_{\varepsilon}\) are the turbulent Prandtl numbers of \(k\) and \(\varepsilon\), respectively. The model constants are as follows. \[ \begin{array}{ll}C_{1}=1.44,&C_{2}=1.92,&C_{\mu}=0.09\\ \sigma_{k}=1.0,&\sigma_{\varepsilon}=1.3\end{array}  \tag{6}\]  
In addition, \(\mu_{t}\) denotes the turbulent viscosity obtained by  
\[\mu_{t}=C_{\mu}\rho\frac{k^{2}}{\varepsilon}. \tag{7}\]  
Figure 1 shows the schematic of the two-dimensional computational domain of the HNS accident ship. The shape of the ship was provided by Korea Research Institute of Ship and Ocean Engineering (KRISO). It is assumed that the ship moves parallel to the direction of the current. The current scenario considers the specific crack position and leakage flow rate, as illustrated in Figure 1, and assumes that the spilled HNS is mainly affected by the sea current. In spite of the two-dimensional analysis, the present study is expected to provide a useful analysis of HNS propagation affected by the current velocity and substance density. The domain size was 1000 m in length and 100 m in depth (\(y\)). The length of the accident ship (\(I_{ship}\)) was 80 m and its height (\(H_{draft}\)) was 8 m. The crack area was 0.1 m\({}^{2}\), located 30 m from the top of the ship. The number of grids in the computation domain was taken as 355,496, determined from grid-independent tests. In general, the sea current consists of three different components: tide-induced, density-driven, and wind-driven currents. Among them, the density-driven current and tide-induced current have almost uniform distributions in the depth direction [14]. Furthermore, the wind-driven current velocity is smaller than the other two. Thus, the ocean current is assumed to be a simplified uniform velocity profile for calculation. The sea-current velocity was varied in the range from 0. 5 to 1.5 m/s for the present simulation. The Neumann condition was applied to the outlet boundary surface and the slip condition was used for the top surface of the computational domain.