Numerical Simulation of Propagation Characteristics of Hazardous Noxious Substances Spilled from Transport Ships

: This study numerically investigates the propagation characteristics of hazardous noxious substances (HNSs) spilled from transport ships and suggests the metal model for predicting the HNS propagation velocity varied with the current velocity and HNS density. The commercial computational ﬂuid dynamics (CFD) code ANSYS FLUENT (V. 17.2) was used for two-dimensional simulation based on the Reynolds-averaged Navier–Stokes (RANS) equation together with the standard k – ε model. The scalar transport equation was also solved to estimate the spatial and transient behaviors of HNS. The main parameters to analyze the near-ﬁeld propagation characteristics of HNSs spilled from the ship were layer thickness, HNS concentration, and propagation velocity. It was found that advection becomes more dominant in propagating an HNS layer that becomes thinner as the current velocity increases. When the current velocity increased beyond a certain level (~0.75 m/s), the mixing effect made the HNS layer less dense but thicker. Consequently, lower-density HNS causes increased HNS concentrations at sea level. As the current velocity increased, the concentration distribution became homogeneous regardless of HNS density. In particular, the second-order response surface model provided for three variables on the basis of the numerical results for 15 cases with the use of the general least-squares regression method, showing a good ﬁt. This model would be useful in estimating the propagation velocity of HNS spilled from a ship.


Introduction
Hazardous and noxious substances (HNSs) are widely known as liquid and mixed-liquid substances that can have harmful effects on the marine environment [1]. Large quantities of HNSs are transported by sea. If an HNS-related accident occurs in a port or on a ship, harmful substances are diffused in the vicinity of the ships in different ways depending on the type of spill material (e.g., sinking, floating, evaporating, dissolving) [2]. Indeed, it is difficult to cope with HNS accidents because of the high risk associated with various HNSs with bioaccumulative, biodegradable, toxic, and explosive characteristics [3][4][5][6][7][8]. To date, many government authorities and environmental organizations have concentrated their effort on establishing emergency planning against HNS disasters. There is still a lack of information on HNS characteristics under HNS accident scenarios.
Indeed, it is difficult to perform experiments directly on HNS propagation in real circumstances because of the toxicity of HNSs with various chemical properties. Nevertheless, Fuhrer et al. [2]

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.
∂ρ ∂t where ρ, u, x, and t are the density, velocity, position, and time, respectively. p is the pressure and τ 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-ε 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 ε are as follows.
∂ ∂t ∂ ∂t 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 σ k and σ ε are the turbulent Prandtl numbers of k and ε, respectively. The model constants are as follows.
C 1 = 1.44, C 2 = 1.92, C µ = 0.09 σ k = 1.0, σ ε = 1.3 (6) In addition, µ t denotes the turbulent viscosity obtained by 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 (L 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. The floating HNS refers to a substance whose density is less than seawater density (1024 kg/m 3 ), and the HNS density was conventionally in the range from 700 to 900 kg/m 3 [2]. The present simulation varied the HNS density in this range, resulting in floating behavior. The viscosity and diffusion coefficient of the HNS were set to 0.00076 kg/m·s and 8 × 10 −8 m 2 /s, based on styrene. It was assumed that the HNS spilled with a mass flow rate of 45 kg/s·m 2 , estimated from a possible accident scenario. A first-order upwind scheme was used, and the density of the mixture was calculated by using the volume-weighted mixing law, as follows.
The time step was fixed as 0.25 s, and its averaged Courant number during the calculation was 0.183. The total flow time was taken as 1000 s for all calculations.

Response Surface Model
The meta-model can be suggested for providing a surrogate mathematical model of the original CFD simulation [15][16][17]. The present study used the response surface model (RSM) [17][18][19]. First, the surface function of interest was set as where y(x) is the surface function of interest, f(x) is a polynomial function of x, and e is the error or noise observed in y(x). The polynomial function f(x) used to approximate y(x) is assumed to be of typical quadratic form, expressed by The regression coefficients, β0, βi, βii, and βij, were determined by the least-squares method that minimizes the sum of the squares of the residual between ' y and y(x). The regression coefficients can be determined by where y is an n × 1 column vector of the response values, X is a k × n matrix of the sample points, and X′ is its transpose. 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. The floating HNS refers to a substance whose density is less than seawater density (1024 kg/m 3 ), and the HNS density was conventionally in the range from 700 to 900 kg/m 3 [2]. The present simulation varied the HNS density in this range, resulting in floating behavior. The viscosity and diffusion coefficient of the HNS were set to 0.00076 kg/m·s and 8 × 10 −8 m 2 /s, based on styrene. It was assumed that the HNS spilled with a mass flow rate of 45 kg/s·m 2 , estimated from a possible accident scenario. A first-order upwind scheme was used, and the density of the mixture was calculated by using the volume-weighted mixing law, as follows.
The time step was fixed as 0.25 s, and its averaged Courant number during the calculation was 0.183. The total flow time was taken as 1000 s for all calculations.

Response Surface Model
The meta-model can be suggested for providing a surrogate mathematical model of the original CFD simulation [15][16][17]. The present study used the response surface model (RSM) [17][18][19]. First, the surface function of interest was set as where y(x) is the surface function of interest, f (x) is a polynomial function of x, and e is the error or noise observed in y(x). The polynomial function f (x) used to approximate y(x) is assumed to be of typical quadratic form, expressed by The regression coefficients, β 0 , β i , β ii , and β ij , were determined by the least-squares method that minimizes the sum of the squares of the residual between y and y(x). The regression coefficients can be determined by where y is an n × 1 column vector of the response values, X is a k × n matrix of the sample points, and X is its transpose.

Validation of CFD Simulation
Unfortunately, it was not possible to find any experimental or numerical results for the abovementioned scenario. As an alternative, additional CFD simulations were conducted to evaluate the flow-rate estimation made in this study by comparing the results from a previous report [20]. This is because the flow rate model becomes an important boundary condition that strongly affects the HNS propagation. Indeed, Tavakoli et al.
[20] conducted a lab-scale experiment for the oil-spill behavior with the damaged tank below the waterline. According to Reference [20], the present study constructed the computational domain, as shown in Figure 2. The size of the room was 12 m × 5 m × 5 m and the depth of water was 3 m. An oil tank with a size of 1 m × 0.5 m × 1 was placed at the center of the water bath. The draft of the tank refers to the depth of the submerged part of the tank, and its value was taken as 0.6 m in this study. A puncture with a diameter of 2.2 cm was drilled 0.1 m from the bottom of the oil tank. The number of grids in the computational domain was taken to be approximately 370,000, determined from grid-independent tests. The oil tank was fully filled with oil, initially. The top surface was set as an atmospheric-pressure outlet condition.
The oil in the tank began to flow out through the puncture initially, and the volume flow rate was estimated by the change in the oil head in the tank. After the spill began, the oil level and volume flow rate decreased owing to the change in pressure difference. The predicted volume flow rate of spilled oil was compared with the experimental result, as shown in Figure 3. It shows reasonable agreement between the prediction and measurement, showing a difference of less than 8%.

Validation of CFD Simulation
Unfortunately, it was not possible to find any experimental or numerical results for the abovementioned scenario. As an alternative, additional CFD simulations were conducted to evaluate the flow-rate estimation made in this study by comparing the results from a previous report [20]. This is because the flow rate model becomes an important boundary condition that strongly affects the HNS propagation. Indeed, Tavakoli et al.
[20] conducted a lab-scale experiment for the oil-spill behavior with the damaged tank below the waterline. According to Reference [20], the present study constructed the computational domain, as shown in Figure 2. The size of the room was 12 m × 5 m × 5 m and the depth of water was 3 m. An oil tank with a size of 1 m × 0.5 m × 1 was placed at the center of the water bath. The draft of the tank refers to the depth of the submerged part of the tank, and its value was taken as 0.6 m in this study. A puncture with a diameter of 2.2 cm was drilled 0.1 m from the bottom of the oil tank. The number of grids in the computational domain was taken to be approximately 370,000, determined from grid-independent tests. The oil tank was fully filled with oil, initially. The top surface was set as an atmospheric-pressure outlet condition.
The oil in the tank began to flow out through the puncture initially, and the volume flow rate was estimated by the change in the oil head in the tank. After the spill began, the oil level and volume flow rate decreased owing to the change in pressure difference. The predicted volume flow rate of spilled oil was compared with the experimental result, as shown in Figure 3. It shows reasonable agreement between the prediction and measurement, showing a difference of less than 8%.

Validation of CFD Simulation
Unfortunately, it was not possible to find any experimental or numerical results for the abovementioned scenario. As an alternative, additional CFD simulations were conducted to evaluate the flow-rate estimation made in this study by comparing the results from a previous report [20]. This is because the flow rate model becomes an important boundary condition that strongly affects the HNS propagation. Indeed, Tavakoli et al.
[20] conducted a lab-scale experiment for the oil-spill behavior with the damaged tank below the waterline. According to Reference [20], the present study constructed the computational domain, as shown in Figure 2. The size of the room was 12 m × 5 m × 5 m and the depth of water was 3 m. An oil tank with a size of 1 m × 0.5 m × 1 was placed at the center of the water bath. The draft of the tank refers to the depth of the submerged part of the tank, and its value was taken as 0.6 m in this study. A puncture with a diameter of 2.2 cm was drilled 0.1 m from the bottom of the oil tank. The number of grids in the computational domain was taken to be approximately 370,000, determined from grid-independent tests. The oil tank was fully filled with oil, initially. The top surface was set as an atmospheric-pressure outlet condition.
The oil in the tank began to flow out through the puncture initially, and the volume flow rate was estimated by the change in the oil head in the tank. After the spill began, the oil level and volume flow rate decreased owing to the change in pressure difference. The predicted volume flow rate of spilled oil was compared with the experimental result, as shown in Figure 3. It shows reasonable agreement between the prediction and measurement, showing a difference of less than 8%.    The mass fraction of HNS in the wake region was measured in the range from 0.02 to 0.1. The dissolved HNS floated near the sea surface because it had a lower density than seawater. Here, the HNS layer formed behind the wake was mainly analyzed. The maximum mass fraction appeared at sea level, and it showed decreasing patterns in the depth direction. For quantitative analysis, three main parameters, namely, the HNS layer thickness (δ layer ), HNS propagation velocity (v p ), and averaged HNS mass fraction at sea level (φ s−l ) were introduced in the present study.

HNS Propagation Characteristics
Appl. Sci. 2018, 8, 2409 6 of 12 Figure 4 shows the distributions of HNS mass fractions at t = 1000 s after leakage. After 1000 s, the flow field reached the steady state, and the shape of the HNS layer was almost constant over time. The mass fraction of HNS in the wake region was measured in the range from 0.02 to 0.1. The dissolved HNS floated near the sea surface because it had a lower density than seawater. Here, the HNS layer formed behind the wake was mainly analyzed. The maximum mass fraction appeared at sea level, and it showed decreasing patterns in the depth direction. For quantitative analysis, three main parameters, namely, the HNS layer thickness (δlayer), HNS propagation velocity (vp), and averaged HNS mass fraction at sea level ( φ s−l) were introduced in the present study.    Figure 5a represents the change in the HNS layer thickness for different current velocities and HNS densities at t = 1000 s. The HNS layer thickness is defined as the depth at which the integrated HNS mass reaches 90% of the total amount of spilled HNS by integrating the mass fraction of HNS (φ hns ) from sea level at t = 1000 s, as follows.

HNS Propagation Characteristics
where x max and y max represent the length of the computational domain in the x and y directions, respectively. The denominator represents the total amount of spilled HNS. The HNS layer thickness indicates the upper limit for the y-axis of the integral in the numerator, as shown in Figure 5b. From the results, the lower the HNS density, the thinner the layer that was formed under the same current velocity conditions, because of the buoyancy effect. When the density of HNS decreases, the dissolved HNS floats rapidly to the sea surface and forms a thinner layer. The change in HNS layer thickness is presented in Figure 5a.
where xmax and ymax represent the length of the computational domain in the x and y directions, respectively. The denominator represents the total amount of spilled HNS. The HNS layer thickness indicates the upper limit for the y-axis of the integral in the numerator, as shown in Figure 5b. From the results, the lower the HNS density, the thinner the layer that was formed under the same current velocity conditions, because of the buoyancy effect. When the density of HNS decreases, the dissolved HNS floats rapidly to the sea surface and forms a thinner layer. The change in HNS layer thickness is presented in Figure 5a. The HNS layer thickness decreases with the increase in current velocity, but it increases after a certain value of current velocity because of the competitive role between advection and diffusion. The inflection point appears at approximately vc = 0.75 m/s. When vc is 0.5 m/s, the corresponding δlayer values are 4.56, 4.21, and 4.01 m for ρhns values of 700, 800, and 900 kg/m 3 , respectively. As vc increases, the advection becomes dominant to overcome the diffusion effect and acts as the main driving force in HNS propagation. However, when vc increases above a certain value (~0.75 m/s) at which an inflection point appears, the mixing of HNS becomes significant and makes the HNS layer less dense but thicker.
where t and v x are the time and velocity of the HNS interface in the downstream direction. T denotes the period of integration, which is adopted as the whole computation time in the present study. Figure 6b shows the maximum propagation length of the HNS interface in the downstream direction. The maximum propagation length linearly increases with time and is significantly affected by the current velocity. As shown in Figure 6c, however, very similar tendencies are found for different densities because the advection and diffusion effects become more dominant than the buoyancy effect. The propagation velocity v p is slightly greater than v c . For example, in the case of v c = 1.5 m/s, v p is approximately 1.6 m/s. In addition, v p varies with HNS density. When v c is 0.5 m/s, v p is 0.56, 0.60 and 0.65 m/s with respect to the HNS density of 700, 800, and 900 kg/m 3 , respectively. The reason is as follows. The dissolved HNS floats with vertical-directional momentum because of the buoyancy-driven force. The vertical-directional momentum of dissolved HNS will change to the streamwise-directional momentum at the sea surface region. Therefore, floating HNS with a density of 700 kg/m 3 propagates faster than that with a density of 900 kg/m 3 . This effect is valid even in the case of v c = 1.5 m/s but is not significant at higher-velocity conditions. This is because the vertical direction momentum component of the dissolved HNS decreases as v c increases.
where t and vx are the time and velocity of the HNS interface in the downstream direction. T denotes the period of integration, which is adopted as the whole computation time in the present study. Figure 6b shows the maximum propagation length of the HNS interface in the downstream direction. The maximum propagation length linearly increases with time and is significantly affected by the current velocity. As shown in Figure 6c, however, very similar tendencies are found for different densities because the advection and diffusion effects become more dominant than the buoyancy effect. The propagation velocity vp is slightly greater than vc. For example, in the case of vc = 1.5 m/s, vp is approximately 1.6 m/s. In addition, vp varies with HNS density. When vc is 0.5 m/s, vp is 0.56, 0.60, and 0.65 m/s with respect to the HNS density of 700, 800, and 900 kg/m 3 , respectively. The reason is as follows. The dissolved HNS floats with vertical-directional momentum because of the buoyancydriven force. The vertical-directional momentum of dissolved HNS will change to the streamwisedirectional momentum at the sea surface region. Therefore, floating HNS with a density of 700 kg/m 3 propagates faster than that with a density of 900 kg/m 3 . This effect is valid even in the case of vc = 1.5 m/s but is not significant at higher-velocity conditions. This is because the vertical direction momentum component of the dissolved HNS decreases as vc increases.  Figure 7 represents the change in the averaged HNS mass fraction at sea level for different current velocity and HNS density. The averaged HNS mass fraction at the sea level, φ s−l , is defined as the area-averaged mass fraction of the HNS at sea level (y = 0) as follows: where φ(x) represents the HNS mass fraction at sea-level and A is the area at the sea surface. As v c increases, φ s−l decreases proportionally, since spilled HNS is not stagnant but is actively propagated by advection. From 1 to 1.5 m/s v c , the averaged HNS mass fraction at sea level increased for less dense HNS because buoyancy is a significant driving force for floating HNS at sea level. However, when v c is 0.5 m/s, φ s−l is estimated to be about 0.021 regardless of HNS density. This is because the spilled HNS with a density of 700 kg/m 3 propagates faster than that with a density 900 kg/m 3 due to the buoyancy force, and it makes the HNS layer less dense. Therefore, the φ s−l values become similar regardless of HNS density at a v c of 0.5 m/s.  Figure 7 represents the change in the averaged HNS mass fraction at sea level for different current velocity and HNS density. The averaged HNS mass fraction at the sea level, φ s−l, is defined as the area-averaged mass fraction of the HNS at sea level (y = 0) as follows: where φ (x) represents the HNS mass fraction at sea-level and A is the area at the sea surface. As vc increases, φ s−l decreases proportionally, since spilled HNS is not stagnant but is actively propagated by advection. From 1 to 1.5 m/s vc, the averaged HNS mass fraction at sea level increased for less dense HNS because buoyancy is a significant driving force for floating HNS at sea level. However, when vc is 0.5 m/s, φ s−l is estimated to be about 0.021 regardless of HNS density. This is because the spilled HNS with a density of 700 kg/m 3 propagates faster than that with a density 900 kg/m 3 due to the buoyancy force, and it makes the HNS layer less dense. Therefore, the φ s−l values become similar regardless of HNS density at a vc of 0.5 m/s.

Response Surface Model for HNS Spill Predictions
The CFD simulation provides deterministic solutions that are bounded by the boundary and operating conditions, which vary with different scenarios. Thus, obvious limitations exist in taking the CFD approach whenever accident scenarios are changed. To reduce the computational cost, the statistical approach has been widely used in many engineering applications and is also used in the present study. As in the meta-modeling, the present study adopted the RSM to provide the surrogate mathematical model from the CFD simulation results. To our knowledge, this is the first study of the near-field propagation characteristics. Of course, it is not easy to create a meta-model with high accuracy because there are insufficient simulation results. However, a more accurate model can be made by accumulating the simulation data and updating the model continuously. The main objective of this work is, therefore, to construct the numerical procedure required for making the meta-model of HNS propagation. The three main parameters, HNS layer thickness (δlayer), HNS propagation velocity (vp), and averaged HNS mass fraction at sea level (φ s−l) are the dependent variables used in

Response Surface Model for HNS Spill Predictions
The CFD simulation provides deterministic solutions that are bounded by the boundary and operating conditions, which vary with different scenarios. Thus, obvious limitations exist in taking the CFD approach whenever accident scenarios are changed. To reduce the computational cost, the statistical approach has been widely used in many engineering applications and is also used in the present study. As in the meta-modeling, the present study adopted the RSM to provide the surrogate mathematical model from the CFD simulation results. To our knowledge, this is the first study of the near-field propagation characteristics. Of course, it is not easy to create a meta-model with high accuracy because there are insufficient simulation results. However, a more accurate model can be made by accumulating the simulation data and updating the model continuously. The main objective of this work is, therefore, to construct the numerical procedure required for making the meta-model of HNS propagation. The three main parameters, HNS layer thickness (δ layer ), HNS propagation velocity (v p ), and averaged HNS mass fraction at sea level (φ s−l ) are the dependent variables used in the mathematical model. The second-order RSM for the three main parameters is fitted to 15 cases using a general least-squares regression method. Each input/dependent variable is obtained from the CFD simulation, as discussed above. The general form of second-order RSM is expressed by where, y is the main parameter, x 1 = v c , and x 2 = ρ hns ; a 0 to a 5 are model constants. The obtained model constants are listed in Table 1. Three-dimensional response surface contours of each dependent variable are shown in Figure 8. The predicted values are compared with the corresponding CFD simulation results (black dots), and the goodness of fit is measured. The variables and scatter of the data points around the fitted regression curve are always evaluated between 0 and 1. As listed in Table 2, both values are greater than 0.9, which means that the prediction model provides a good fit.
where, y′ is the main parameter, x1 = vc, and x2 = ρhns; a0 to a5 are model constants. The obtained model constants are listed in Table 1. Three-dimensional response surface contours of each dependent variable are shown in Figure 8. The predicted values are compared with the corresponding CFD simulation results (black dots), and the goodness of fit is measured. The variables and scatter of the data points around the fitted regression curve are always evaluated between 0 and 1. As listed in Table 2, both values are greater than 0.9, which means that the prediction model provides a good fit.

Conclusions
The present study conducted a numerical simulation of an HNS spill and developed a prediction model. The HNS spill and propagation behaviors were evaluated by numerical analysis according to the water-current velocity conditions and HNS density. The following three main parameters were defined: HNS layer thickness, HNS propagation velocity, and averaged HNS mass fraction at sea level. In addition, corresponding mathematical prediction models were developed using the response surface-model method. The results are as follows, and such data would be useful in HNS prediction technology and understanding the detailed physics behind the near-field characteristics of HNS spills and propagation.
(1) Lower-density HNS forms thinner layers in all water-current velocity ranges. As current velocity increases, the advection effect becomes a significant driving force of HNS propagation. In this case, the spilled HNS was not stagnant, but was actively propagated by advection, and the layer

Conclusions
The present study conducted a numerical simulation of an HNS spill and developed a prediction model. The HNS spill and propagation behaviors were evaluated by numerical analysis according to the water-current velocity conditions and HNS density. The following three main parameters were defined: HNS layer thickness, HNS propagation velocity, and averaged HNS mass fraction at sea level. In addition, corresponding mathematical prediction models were developed using the response surface-model method. The results are as follows, and such data would be useful in HNS prediction technology and understanding the detailed physics behind the near-field characteristics of HNS spills and propagation.
(1) Lower-density HNS forms thinner layers in all water-current velocity ranges. As current velocity increases, the advection effect becomes a significant driving force of HNS propagation. In this case, the spilled HNS was not stagnant, but was actively propagated by advection, and the layer thickness decreased. However, when the current velocity increased above a certain level (~0.75 m/s), the mixing effect became significant and made the HNS layer less dense, but the layer thickness increased.
(2) The current velocity strongly affected the propagation velocity, which linearly increased in proportion to current velocity. The averaged HNS mass fraction at sea level increased for less dense HNS because buoyancy was a significant driving force for floating HNS at sea level. Moreover, the HNS mass fraction decreased as the current velocity increased.