Development of a Metamodel for Predicting Near-Field Propagation of Hazardous and Noxious Substances Spilled From a Ship

: This study aims to numerically analyze the near-ﬁeld propagation behavior of hazardous and noxious substances (HNSs) and to develop a new metamodel for HNS propagation. Extensive computational ﬂuid dynamics (CFD) simulations were conducted using the ANSYS FLUENT (V. 17.2) code for various HNS spill scenarios. We newly introduced several key parameters, including the streamwise propagation velocity, transverse propagation velocity, and averaged HNS mass fraction. From the results, the advection e ﬀ ect is more dominant with an increase in the current velocity and streamwise propagation velocity, and with a decrease in the transverse propagation velocity. Also, the HNS mass fraction decreases as the current velocity increases with the change of concentration and propagation area. Particularly, a new metamodel of HNS propagation based on the current CFD results was validated by the hidden point test, showing very good ﬁt. We believe this model would make useful predictions under various scenarios without CFD simulations.


Introduction
Hazardous and noxious substances (HNSs) have harmful effects on human health, biological resources, and the marine environment when they are spilled into the marine environment [1]. The volume of maritime trade of HNSs is increasing every year, thus increasing the possibility of HNS spill at sea. HNSs are distinguished by one or more of the following five characteristics according to their inherent properties: Flammable, explosive, toxic, corrosive, and reactive. Owing to these hazardous characteristics, a response strategy to, and rapid removal of, a HNS spill must be determined. However, it is difficult to cope with a HNS spill quickly and accurately because of its complicated characteristics [2][3][4][5][6][7]. Furthermore, there is a lack of information on HNS propagation near ships over time under various accident scenarios. Thus, a more elaborate prediction method to measure HNS propagation characteristics is required that accounts for various effects of the HNS properties and environmental conditions. An experiment may be the best way to obtain the detailed propagation characteristics when an HNS spills in the ocean. Fuhrer et al. [8] conducted experiments under the limited condition, and made floating cells to monitor the residual concentration over time. In their experiment, they released a small amount of styrene, then captured the HNS behavior for several hours and a few days after. This result was taken under very limited conditions but it can be used for creating an emergency response guide for spill accidents. Despite this, due to the hazardous chemical properties, the study did not consider realistic scenarios, as releasing a HNS into the ocean may cause danger to the aquatic ecosystem as well as human safety.
A computational study could be a promising alternative to the experiment for providing useful information on the propagation behavior of HNSs in more detail. Janeiro et al. [9] numerically studied the management of oil pollution accidents in the Tuscany Archipelago region using a set of nested models with downscaled conditions. As a result, the model showed good agreement with the trajectories, captured by satellite. Elhakeem et al. [10] established three-dimensional numerical models by combining the hydrodynamic and oil spill models. They evaluated the predictions by comparison with the observed results and used the model to build a risk assessment map for an oil spill accident. Berry et al. [11] developed the semi-empirical model, called OILTRANS, with theoretical equations and an empirical database to predict the spatial propagation. For oil spills, this model has been known to be accurate when compared with experimental observation. To the best of our knowledge, there are few computational fluid dynamics (CFD) studies for HNS propagation, especially regarding near-field characteristics, from the earlier published literature. Thus, more detailed information on the spatial change of HNS concentration over time is required to build a response strategy to an HNS spill, particularly in the early stage of the accident.
As emphasized above, the experimental approach found it difficult to obtain useful and sufficient data related to HNS propagation, whereas the CFD approach could be useful for predicting the transient evolution of HNS concentration around a ship. If an accident occurs, all possible responses, e.g., setting up the risk zone, separation, evacuation, and chemical treatment, would be decided as quickly as possible. However, considerable additional time is required after the accident to provide HNS distributions over time. To overcome this limitation, it is necessary to develop a metamodel that predicts all information on the HNS propagation more quickly without CFD simulation. Furthermore, this model should be based on the CFD results for all possible scenarios, taken from literature and accident reports.
Meanwhile, most simulations have been conducted for a space of a few hundred kilometers. Hence, information within 1 km of the ship cannot be yielded, and, accordingly, cannot contribute to the early response stage of the HNS spill in detail. Previously, our group conducted a two-dimensional CFD simulation for predicting the near-field propagation behavior in the vicinity of the ship, within a few hundred meters [12], and presented the spatial HNS distribution over time.
Therefore, the current study aims to examine HNS propagation characteristics by conducting three-dimensional and transient CFD simulations with the key parameters of current velocity, HNS density, and crack location. Based on the Kriging model and the CFD results for 18 test cases, a new metamodel is suggested for predicting the near-field propagation of HNS spilled from a ship. Moreover, we newly define the key parameters, such as the propagation velocity in a streamwise/transverse direction and the averaged HNS mass fraction at sea level, to aid the analysis of the change in HNS concentration over time, in a quantitative manner.

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-ε 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 Tavakoli et al. [13] who provided the 3 of 14 measured data for oil propagation in the oil tank under different scenarios. The governing equations are expressed as follows: where ρ denotes the density, µ 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, ε, are expressed as follows: where µ 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 σ k and σ ε are the turbulent Prandtl numbers of k and ε, respectively. The model constants were taken as 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]. Hence, the present study was set under the following specific conditions: First, the whole computational domain is 1000 m in length (x), 1000 m in width (y), and 100 m in depth (z). The ship is 160 m long, 30 m wide, and has a 14.5 m draft. In the cases of a side and bottom crack, the area was set to 0.25π m 2 located 80 m from the bow of the ship. In addition, the current scenario assumed that the accidental crack only occurred at the side and bottom surfaces of the ship. The side crack was located 10 m below sea level, while the bottom crack was located at the center of the ship. Besides this, we conducted grid-independent tests to ensure grid convergence. As a result, the selected number of grids was 293,162, which were applied for all computational cases. It is necessary to use the appropriate boundary and initial conditions for CFD simulations. In fact, the sea-current (or current) velocity is affected by tidal, wind-driven, and density-driven currents [16]. Due to the lack of information regarding this variation, it may be difficult to consider all factors accurately in simulation. Practically, the current velocity is assumed to be uniformly distributed with a mean value. The current velocity (vc) in the range of 0.514-1.542 m/s was used as the inlet velocity. Accordingly, the Reynolds number was estimated to be larger than 0.51 × 1010 with a scale length of 100 km, which has usually been used for large-scale ocean currents [17]. For simplicity, the free surface flow was not considered and the HNS propagation at sea level was assumed to not be affected by the interfacial motion. The slip condition was used at sea level, assuming the normal gradient of the tangential velocity component becomes zero. Similarly, a zerogradient pressure condition was set to all the outlet surfaces, as the type of Neumann condition. As mentioned above, there is a lack of sufficient information on HNS properties. To examine the density effect on the propagation behavior, we set such artificial cases to show the influence of the HNS density on the propagation behavior. In fact, there are a lot of factors that can affect the HNS spill behavior. Thus, it is not possible to consider all possible scenarios on CFD simulations. Here, we took two important parameters such as current flow velocity and HNS density in CFD simulations for developing the metamodel. Styrene is one of the top 20 imports and exports of chemicals in Korea. Therefore, this study used styrene, which has a high possibility of spillage in case of spill accident. Styrene shows floating behavior and was used as the base material for the simulation. Its density change was considered to be in the range of 700-900 kg/m 3 without changes to other properties, such as the viscosity, and the diffusion coefficients were 0.00076 kg/m·s and 8·10-8 m 2 /s. It was also assumed that the HNS initially spills with a mass flow rate of 900 kg/s at the crack of the ship. For the simulation, the timestep was fixed to 0.05 s and the averaged Courant number was estimated to be 0.1231 for a total flow time of 1000 s. Table 1 summarizes the different cases for this simulation. It is necessary to use the appropriate boundary and initial conditions for CFD simulations. In fact, the sea-current (or current) velocity is affected by tidal, wind-driven, and density-driven currents [16]. Due to the lack of information regarding this variation, it may be difficult to consider all factors accurately in simulation. Practically, the current velocity is assumed to be uniformly distributed with a mean value. The current velocity (v c ) in the range of 0.514-1.542 m/s was used as the inlet velocity. Accordingly, the Reynolds number was estimated to be larger than 0.51 × 1010 with a scale length of 100 km, which has usually been used for large-scale ocean currents [17]. For simplicity, the free surface flow was not considered and the HNS propagation at sea level was assumed to not be affected by the interfacial motion. The slip condition was used at sea level, assuming the normal gradient of the tangential velocity component becomes zero. Similarly, a zero-gradient pressure condition was set to all the outlet surfaces, as the type of Neumann condition. As mentioned above, there is a lack of sufficient information on HNS properties. To examine the density effect on the propagation behavior, we set such artificial cases to show the influence of the HNS density on the propagation behavior. In fact, there are a lot of factors that can affect the HNS spill behavior. Thus, it is not possible to consider all possible scenarios on CFD simulations. Here, we took two important parameters such as current flow velocity and HNS density in CFD simulations for developing the metamodel. Styrene is one of the top 20 imports and exports of chemicals in Korea. Therefore, this study used styrene, which has a high possibility of spillage in case of spill accident. Styrene shows floating behavior and was used as the base material for the simulation. Its density change was considered to be in the range of 700-900 kg/m 3 without changes to other properties, such as the viscosity, and the diffusion coefficients were 0.00076 kg/m·s and 8·10-8 m 2 /s. It was also assumed that the HNS initially spills with a mass flow rate of 900 kg/s at the crack of the ship. For the simulation, the timestep was fixed to 0.05 s and the averaged Courant number was estimated to be 0.1231 for a total flow time of 1000 s. Table 1 summarizes the different cases for this simulation.

Metamodeling of HNS Propagation
As mentioned above, the present study developed a metamodel for predicting HNS propagation based on the CFD results. The well-known Kriging model was adopted for to make the metamodel [18][19][20]. Basically, the unknown function y(x) includes the regression model and the stochastic process and has two terms, as follows: where f (x) represents the approximation function, and Z(x) indicates the localized deviation interpolating the n sampled values. Here, the mean value of the variance was assumed to be zero. The correlation between Z(x i ) and Z(x j ) is given by the covariance matrix: where R is the correlation matrix and σ 2 is the process variance. R also indicates Gaussian correlation function expressed by Consequently, the Kriging predictorŷ(x) can be expressed at an arbitrary point x as follows: where y represents the CFD simulation result with a column vector of length n, f is a unit column vector of length n s , and r T refers to a column vector of length n between each point {x 1 , ..., x n } and arbitrary points x. The column vector is expressed as follows: In Equation (11),β is obtained by ( f T R −1 f ) −1 f T R −1 y and the correlation parameter θ k in Equation (10) is determined by using the Maximum likelihood estimates (MLE), expressed by Indeed, θ k was iteratively determined for each variable based on the CFD simulation, requiring a large computational cost. Figure 2 shows the time-dependent propagation behavior of the HNS mass fraction field. This study determined the interfaces between the HNS and seawater with a mass fraction of 100 ppm taken as the criterion with reference to the exposure limit for styrene, as selected by the Occupational Safety and Health Administration [21]. During propagation over time, there was little difference between the side and bottom crack cases, noticing that the crack location has practically negligible effect. This indicates that the current velocity and HNS property may be important parameters in HNS propagation.

HNS Propagation Behaviors
Indeed, θk was iteratively determined for each variable based on the CFD simulation, requiring a large computational cost. Figure 2 shows the time-dependent propagation behavior of the HNS mass fraction field. This study determined the interfaces between the HNS and seawater with a mass fraction of 100 ppm taken as the criterion with reference to the exposure limit for styrene, as selected by the Occupational Safety and Health Administration [21]. During propagation over time, there was little difference between the side and bottom crack cases, noticing that the crack location has practically negligible effect. This indicates that the current velocity and HNS property may be important parameters in HNS propagation.  Figure 3a compares the HNS propagation characteristics with respect to the current velocity for different crack locations. As the current velocity increases, the propagation becomes faster owing to higher convective motions, resulting in relatively lower diffusive behavior and rapid propagation in the streamwise direction. As shown in Figure 3b, it is interesting that even for a heavier HNS, a big difference is not observed, indicating that for the floating material, such as styrene, convective motion becomes dominant. Under the present scenario, the spilled HNS interacted with the vertical directional (z direction) momentum only by the buoyancy force. When the HNS moves towards sea level because of floating behavior, the vertical motion rapidly changes owing to the streamwise motion, which is dominant. In addition, from Figure 3, it can be seen that there is no clear difference  Figure 3a compares the HNS propagation characteristics with respect to the current velocity for different crack locations. As the current velocity increases, the propagation becomes faster owing to higher convective motions, resulting in relatively lower diffusive behavior and rapid propagation in the streamwise direction. As shown in Figure 3b, it is interesting that even for a heavier HNS, a big difference is not observed, indicating that for the floating material, such as styrene, convective motion becomes dominant. Under the present scenario, the spilled HNS interacted with the vertical directional (z direction) momentum only by the buoyancy force. When the HNS moves towards sea level because of floating behavior, the vertical motion rapidly changes owing to the streamwise motion, which is dominant. In addition, from Figure 3, it can be seen that there is no clear difference in the propagation behaviors according to the crack location, except for the asymmetric propagation in the side crack cases.

HNS Propagation Behaviors
in the propagation behaviors according to the crack location, except for the asymmetric propagation in the side crack cases.

Quantitative Analysis
For quantitative analysis, the present study introduces three key parameters: The streamwise propagation velocity ( ̅ ), transverse propagation velocity ( ̅ ), and averaged HNS mass fraction

Quantitative Analysis
For quantitative analysis, the present study introduces three key parameters: The streamwise propagation velocity (v x ), transverse propagation velocity (v y ), and averaged HNS mass fraction (w s−l ). First, v x and v y are the time-averaged propagation velocities at the HNS interface in the x direction and the y direction, respectively. (15) where T is the total time at which HNS interface reaches the outlet region in the computational domain. We also suggest a new parameter, the area-averaged mass fraction of the HNS at sea level (z = 0), given by: where w and A hns represent the local mass fraction of the HNS and the area covered by the HNS at sea level. The parameters defined above would be useful for understanding the HNS propagation behavior in a quantitative manner. As far as we know, this is the first effort for a quantitative analysis of HNS propagation. Figure 4a shows the behavior of the HNS interface with respect to the current velocity for 200 s. The convective effect becomes stronger as the current velocity increases. For instance, at v c = 1.542 m/s, the local interface velocity may be approximately 1.6 m/s, which is slightly faster than the current velocity because of the flow change induced by the shape of the transport ship, as well as the buoyancy effect. Figure 4b,c provides the estimated propagation distance (D p ) of the HNS interface in the streamwise and transverse directions according to the current velocity, respectively. It is seen that the streamwise propagation substantially increases, while the propagation in the transverse direction rather decreases with the increase in the current velocity because the total flow momentum is redistributed in both directions. ). First, ̅ and ̅ are the time-averaged propagation velocities at the HNS interface in the x direction and the y direction, respectively.
where T is the total time at which HNS interface reaches the outlet region in the computational domain. We also suggest a new parameter, the area-averaged mass fraction of the HNS at sea level (z = 0), given by: where w and Ahns represent the local mass fraction of the HNS and the area covered by the HNS at sea level. The parameters defined above would be useful for understanding the HNS propagation behavior in a quantitative manner. As far as we know, this is the first effort for a quantitative analysis of HNS propagation. Figure 4a shows the behavior of the HNS interface with respect to the current velocity for 200 s. The convective effect becomes stronger as the current velocity increases. For instance, at vc = 1.542 m/s, the local interface velocity may be approximately 1.6 m/s, which is slightly faster than the current velocity because of the flow change induced by the shape of the transport ship, as well as the buoyancy effect. Figure 4b,c provides the estimated propagation distance (Dp) of the HNS interface in the streamwise and transverse directions according to the current velocity, respectively. It is seen that the streamwise propagation substantially increases, while the propagation in the transverse direction rather decreases with the increase in the current velocity because the total flow momentum is redistributed in both directions.   Figure 5a represents the effects of the current velocity on ̅ for different HNS densities. The results show that ̅ is proportional to the current velocity, and no difference can be seen according to the HNS density. This is also seen in Figure 3b. The reason for this is because the advection effect becomes more dominant in the streamwise direction compared to the buoyancy effect. The maximum propagation distance is shown in Figure 5b. From the results, there is no difference in ̅ between the cases of the side crack and bottom crack. This indicates that the crack position is no longer involved, and its effect is negligible.    Figure 5a represents the effects of the current velocity on v x for different HNS densities. The results show that v x is proportional to the current velocity, and no difference can be seen according to the HNS density. This is also seen in Figure 3b. The reason for this is because the advection effect becomes more dominant in the streamwise direction compared to the buoyancy effect. The maximum propagation distance is shown in Figure 5b. From the results, there is no difference in v x between the cases of the side crack and bottom crack. This indicates that the crack position is no longer involved, and its effect is negligible.   Figure 5a represents the effects of the current velocity on ̅ for different HNS densities. The results show that ̅ is proportional to the current velocity, and no difference can be seen according to the HNS density. This is also seen in Figure 3b. The reason for this is because the advection effect becomes more dominant in the streamwise direction compared to the buoyancy effect. The maximum propagation distance is shown in Figure 5b. From the results, there is no difference in ̅ between the cases of the side crack and bottom crack. This indicates that the crack position is no longer involved, and its effect is negligible.  Figure 6a shows the effects of the current velocity on ̅ for different HNS densities. As the current velocity increases, ̅ decreases, because the advection effect becomes dominant with high current velocity. Similar to Figure 4c, the propagation in the transverse direction decreases owing to the increase in HNS density. Meanwhile, Figure 6b shows the maximum propagation distance according to the crack position. For the side crack, the averaged transverse propagation distance was estimated to be approximately 30 m larger than that for the bottom crack. This deviation could be made by the geometrical effect induced by the ship shape. For the bottom crack, the spilled HNS flows along the bottom surface of the ship. Then, it moves upward and reaches sea level.   Figure 7a shows the effect of the current velocity on the averaged HNS mass fraction at sea level, , according to the HNS density. As expected, decreases as the current velocity increases, showing faster propagation. The mass fraction at sea level increases with HNS density. For example, when vc is 1.028 m/s, increases by approximately 30%. This is because a lower HNS density causes faster HNS propagation velocity in the transverse direction, resulting in wider area occupied by the HNS. Thus, the heavier the HNS, the wider the area it occupies at sea level. In addition, the influence of the crack location on the mass fraction becomes negligible as the current velocity increases.    Figure 7a shows the effect of the current velocity on the averaged HNS mass fraction at sea level, , according to the HNS density. As expected, decreases as the current velocity increases, showing faster propagation. The mass fraction at sea level increases with HNS density. For example, when vc is 1.028 m/s, increases by approximately 30%. This is because a lower HNS density causes faster HNS propagation velocity in the transverse direction, resulting in wider area occupied by the HNS. Thus, the heavier the HNS, the wider the area it occupies at sea level. In addition, the influence of the crack location on the mass fraction becomes negligible as the current velocity increases.

Metamodeling Using Kriging Model
After the HNS accident occurs, the response must be made as fast as possible. At this time, it would be meaningless to perform a CFD simulation practically. Hence, another focus of this study is on the development of a metamodel to predict the HNS propagation. This model would be useful in making the maneuver, especially at the early stage of the HNS-spill accident. Of course, the current study has the limitation of not including all possible scenarios. In this regard, the current study is at the first step of establishing the metamodel platform based on the well-known Kriging model and the CFD results.
Using Equations (11)-(14), we considered 18 numerical cases, as listed in Table 1, to obtain the parameters, such as θ k andβ, for the input variables. Here, subscript k indicates the input variable, such as current velocity, current density, and crack location. Finally, the present study determined the coefficients involved in the prediction model, as summarized in Table 2. Here, we alternatively performed hidden point tests for evaluation because there was no experimental and numerical data for comparison. In doing so, the CFD simulations for eight cases, listed in Table 3, different from the cases listed in Table 1, were carried out. Their results were compared with the estimation of the metamodel proposed in this work. Figure 8 presents the comparison between the CFD results and the hidden point test results. The black dots represent the CFD results of 18 cases (Table 1), and the red dots indicate the additional CFD results for eight cases. Moreover, the surface was formed by the solutions of the metamodel. For evaluation, the hidden point error δ was checked by: where V is the result of the prediction model, and V is the result of additional cases. As a result of the comparison using Equation (17), the average error of the prediction model for each parameter was estimated to be 8.48% in the streamwise propagation velocity, 5.38% in the transverse propagation velocity, and 6.71% in the averaged HNS mass fraction, respectively. All errors are relatively small, showing that the current model is well fitted. In addition, to improve the accuracy of the prediction model, the eight additional cases were simulated and the metamodel was updated based on the results for better correction. Finally, we obtained Figure 9, which shows the results of the metamodel for the key parameters. Putting the possible inputs, such as HNS density, current velocity, and crack location, into the metamodel proposed in the present study would make the estimated outputs, such as the averaged HNS propagation velocity in the streamwise and transverse directions and the HNS mass fraction at sea level. For evaluation, the hidden point error δ was checked by: where V' is the result of the prediction model, and V is the result of additional cases. As a result of the comparison using Equation (17), the average error of the prediction model for each parameter was estimated to be 8.48% in the streamwise propagation velocity, 5.38% in the transverse propagation velocity, and 6.71% in the averaged HNS mass fraction, respectively. All errors are relatively small, showing that the current model is well fitted. In addition, to improve the accuracy of the prediction model, the eight additional cases were simulated and the metamodel was updated based on the results for better correction. Finally, we obtained Figure 9, which shows the results of the metamodel for the key parameters. Putting the possible inputs, such as HNS density, current velocity, and crack location, into the metamodel proposed in the present study would make the estimated outputs, such as the averaged HNS propagation velocity in the streamwise and transverse directions and the HNS mass fraction at sea level.

Conclusions
The present study performed a numerical analysis of the HNS propagation characteristics based on CFD simulations under limited scenarios. Using the CFD results, it involves developing a new metamodel for predicting the key parameters, which are newly introduced, together with the Kriging model. The following conclusions are drawn:

Conclusions
The present study performed a numerical analysis of the HNS propagation characteristics based on CFD simulations under limited scenarios. Using the CFD results, it involves developing a new metamodel for predicting the key parameters, which are newly introduced, together with the Kriging model. The following conclusions are drawn: (1) From the CFD results, it was found that the streamwise propagation velocity increases linearly with the current velocity because the advection effect is dominant, whereas the transverse propagation velocity decreases. Furthermore, the crack position did not cause a large difference in the streamwise propagation velocity.
(2) As the current velocity increases, the averaged HNS mass fraction at sea level decreases because the HNS is actively transported by advection. In the case of the low HNS density, the buoyancy effect becomes stronger and the propagation in the radial direction is increased. Thus, the area occupied by spilled HNS increases, and the averaged HNS mass fraction decreases. For side cracks, the averaged HNS mass fraction decreases because the spilled HNS has lower area. (3) The metamodel was developed using the Kriging model for the key parameters, which were newly defined for the HNS propagation analysis. The hidden point tests for evaluation show that the mean error of the streamwise propagation velocity is 8.48%, the transverse propagation velocity is 5.38%, and the averaged HNS mass fraction is 6.71%, approximately. The current model is the first version developed under the restricted scenario. Thus, it could be updated continuously by adding further CFD results under various scenario.