Modeling of Accidental Oil Spills at Different Phases of LNG Terminal Construction

: Accidental oil spills not only deteriorate biodiversity but also cause immediate threats to coastal environments. This study quantitatively investigates the initial dispersion of spilled oil using the environmental ﬂuid dynamics code (EFDC) model, loosely coupled with an endorsed oil spill model (MEDSLIK-II) accounting for time-dependent advection, diffusion, and physiochemical weathering of the surface oil slick. Focusing on local contributing factors (i.e., construction activities) to oil dispersion, the current model is applied to likely oil spills occurring at three different phases of the Songdo LNG terminal construction on a reclaimed site in South Korea. Applied phases pose detailed ship collision scenarios generated based on a proposed construction plan of the terminal. The effects of permeable revetments, required for reclamation, on the currents were also investigated and applied in subsequent oil spill modeling. For each scenario, the simulated results showed distinct patterns in the advection, dispersion, and transformation of the oil slick. Oil absorption into the coast, which causes immense damage to the coastal communities, is found to be highly dependent on the tidal currents, volume of oil spilled, and nearby construction activities.


Introduction
The demand for transport facilities has dramatically increased with the rapid economic development of coastal areas, and so has the frequency of oil spill accidents [1]. Moreover, unprecedented weather conditions due to global climate change may increase the possibility of an oil spill in the ocean. These accidental oil spills often have immense impact on marine environments and coastal communities. To mitigate damage, detecting oil slick and predicting its fate (i.e., slick displacement and dispersion) often require a specific and timely intervention at sea by governmental response agencies [2]. The most renowned historical oil spills due to ship collisions include the Exxon Valdez (1989), Atlantic Empress (1979), Amoco Cadiz (1978), and Torrey Canyon (1967). The potential for maritime accidents generally increases nearshore construction sites and harbors where oil transports are frequent, and the various routes of oil tanker and barges are concurrent.
In South Korea, a peninsula where northward land transportations are limited, marine transport allows for efficient transcontinental transportation with low freight rates by enabling simultaneous transportation of large amounts of cargo. Therefore, more than 90% of the trading volume relies on marine transport, with oil traffic accounting for about 30% of all marine transportation [3]. Continuous oil-related accidents have occurred mainly due to frequent traffic and related issues. For example, an oil tanker was stranded on the reef, resulting in an oil spill of about 5000 tons near Yeosu, South Korea, in 1995. In 2007, the Hebei Spirit oil spill accident was occurred near Manripo beach, Taean, in the advection-diffusion and transformation processes [21]. Wind velocity fields, computed based on extreme value analysis utilizing the probability density function of a 5-year easterly wind, were also fed into MEDSLIK-II, accounting for atmospheric conditions ( Figure 1). To investigate initial dispersion, 12-h current velocity fields computed from EFDC during maximum phases of ebb to flood at spring tide were inputted to MEDSLIK-II. For efficient data transfer, a spatial grid resolution of 10 m and a time step of 6 min were set for both EFDC and the MEDSLIK-II. Since the computations in EFDC are based on the Eulerian coordinate system, compared with the Eulerian/Lagrangian coordinate system used in MEDSLIK-II, a coordinate conversion was required and applied in the current study to collocate the nodes. Specifically, MEDSLIK-II firstly computes transport, diffusion, and transformation of surface oil in the Eulerian coordinate system based on the Eulerianbased current fields from EFDC, which are then utilized to compute oil spill model state variables [17], including oil concentration and oil slick state variables. As a result, forecasts of oil spill transport and transformation are enabled through the Lagrangian representation of advection-diffusion processes and oil slick variables for weathering processes. based on extreme value analysis utilizing the probability density function of a 5-year easterly wind, were also fed into MEDSLIK-II, accounting for atmospheric conditions ( Figure  1). To investigate initial dispersion, 12-h current velocity fields computed from EFDC during maximum phases of ebb to flood at spring tide were inputted to MEDSLIK-II. For efficient data transfer, a spatial grid resolution of 10 m and a time step of 6 min were set for both EFDC and the MEDSLIK-II. Since the computations in EFDC are based on the Eulerian coordinate system, compared with the Eulerian/Lagrangian coordinate system used in MEDSLIK-II, a coordinate conversion was required and applied in the current study to collocate the nodes. Specifically, MEDSLIK-II firstly computes transport, diffusion, and transformation of surface oil in the Eulerian coordinate system based on the Eulerian-based current fields from EFDC, which are then utilized to compute oil spill model state variables [17], including oil concentration and oil slick state variables. As a result, forecasts of oil spill transport and transformation are enabled through the Lagrangian representation of advection-diffusion processes and oil slick variables for weathering processes.

Numerical Domain and Selected Scenarios
To provide hydrodynamic conditions, EFDC calculated the tidal and wind-induced currents at each construction stage (Daewoo E&C, 2019). Figure 2a and 2b shows the bathymetry and corresponding orthogonal curvilinear grid system in an expanded domain encompassing approximately 130 km in the east-west and 200 km the south-north directions applied in the current study. The grid size varies from 10 ~ 500 m with finer resolutions set close to the LNG terminal construction site (Figure 2c,d). The total number of valid computational grids is 34,408. Five vertical layers (1: surface; 2~4: middle; and 5: bottom) were constructed through the model. The current EFDC model was run for 30 days with a temporal resolution of 3 s. The modeled tide was decomposed to four major tidal constituents including M2, S2, K1, and O1.
Three ship collision scenarios were generated based on a proposed construction plan of the Songdo LNG terminal that is expected to be built on reclaimed land, as depicted in Figure 3. The three scenarios are essentially representative of the three distinguished construction phases (phase 1, 2, and 3) depending on the time frames of construction. Thus, configurations of the construction vessels and probable hazards were reconstructed accordingly. Similar to the Hebei Spirit oil spill accident, an oil spillage was assumed to have originated from a fuel tanker for each scenario (Figure 3). In phase 1, an oil spill of 29,400 L was initiated near the western revetment due to the collision of the fuel tanker. In phase 2 and phase 3, 23,100 and 10,500 L of oil were spilled, both near the eastern revetment, due to the collision of the fuel tanker. A list of the selected properties of the fuel tanker for

Numerical Domain and Selected Scenarios
To provide hydrodynamic conditions, EFDC calculated the tidal and wind-induced currents at each construction stage (Daewoo E&C, 2019). Figure 2a,b shows the bathymetry and corresponding orthogonal curvilinear grid system in an expanded domain encompassing approximately 130 km in the east-west and 200 km the south-north directions applied in the current study. The grid size varies from 10~500 m with finer resolutions set close to the LNG terminal construction site (Figure 2c,d). The total number of valid computational grids is 34,408. Five vertical layers (1: surface; 2~4: middle; and 5: bottom) were constructed through the model. The current EFDC model was run for 30 days with a temporal resolution of 3 s. The modeled tide was decomposed to four major tidal constituents including M 2 , S 2 , K 1 , and O 1 .
Three ship collision scenarios were generated based on a proposed construction plan of the Songdo LNG terminal that is expected to be built on reclaimed land, as depicted in Figure 3. The three scenarios are essentially representative of the three distinguished construction phases (phase 1, 2, and 3) depending on the time frames of construction. Thus, configurations of the construction vessels and probable hazards were reconstructed accordingly. Similar to the Hebei Spirit oil spill accident, an oil spillage was assumed to have originated from a fuel tanker for each scenario (Figure 3). In phase 1, an oil spill of 29,400 L was initiated near the western revetment due to the collision of the fuel tanker. In phase 2 and phase 3, 23,100 and 10,500 L of oil were spilled, both near the eastern revetment, due to the collision of the fuel tanker. A list of the selected properties of the fuel tanker for each scenario is provided in Table 1. A gradual spill of Low-Sulfur Fuel Oil (LSFO) with a prescribed constant spill rate for a 1-h time span was set, derived from the total volume of the spilled oil being 70% of the tanker capacity for all scenarios. Sample properties of LSFO applied in the current model are a density of 0.91 tonne/m 3 , a viscosity of 95.50 mm 2 /s at 15 • C, and a vapor pressure of 0.184 bar. Note that this gradual spill condition improves the relatively unrealistic oil spill scenarios with an instantaneous total spill condition applied in Lee et al. [3]. To evaluate the model predictability, the predicted tidal elevation and current speeds were compared to the observed data by the Korea Hydrographic and Oceanographic Agency. The validation was executed in the vicinity of the construction site at two points (T1 and T2) for the tidal elevation and three points (PC1, PC2, and PC3) for the current speed, as shown in Figure 2e.

Hydrodynamic Model
The EFDC model solves the three-dimensional primitive equations of motion for turbulent flow in a horizontal curvilinear-orthogonal and vertical σ-stretched coordinate system [23]. Transforming the vertically hydrostatic boundary layer form of the turbulent equations of motion and utilizing the Boussinesq approximation for variable density provides the following forms of continuity equation and momentum equations [22]. The continuity equation is where m x and m y are the horizontal scale factors; H is the sum of the mean water depth (h) and free surface elevation (ζ) H = h + ζ; U and V are the horizontal velocity components in the curvilinear, orthogonal coordinates x and y, respectively; and W is the vertical velocity in the σ-stretched dimensionless vertical coordinate z [22]. The momentum equations are ∂ t m x m y HV + ∂ x m y HUV + ∂ y (m x HVV) + ∂ z m x m y WV + f e m x m y HU where f e is the Coriolis parameter; g is the gravitational acceleration; P is the pressure; A v is the vertical eddy viscosity, determined based on the second moment turbulence closure model [24,25]; Q u and Q v are the source and sink terms, respectively. Further details of the above formulations and numerical techniques employed in the EFDC model are provided in Hamrick [22].

Oil Spill Model
In the current study, MEDSLIK-II [17] is employed to simulate the processes of transportation, diffusion, and transformation of surface oil in seawater. In MEDSLIK-II, a tracer concentration, with units of mass over volume, mixed in the marine environment is formulated as ∂C ∂t where ∂/∂t is the local time-rate-of-change operator, U is the sea current mean field with components (U, V, W), K is the diffusivity tensor which parameterizes the turbulent effects, and r j (C) are the M transformation rates that modify the tracer concentration by means of physical and chemical transformation processes.
where C 1 is the oil concentration solution solely due to the weathering process, while the final time rate of change of C is given by the advection-diffusion acting on C 1 . Considering the oil volume, density ρ, and unit area A, the surface oil concentration (C S ) and the dispersed oil concentration (C D ), with units of kg m −2 , can be defined as where vs. and V D are the surface and dispersed oil volumes. When the surface oil arrives close to the coasts, defined by a reference segment L C , the oil can be absorbed and the concentration of oil at the coasts, C C , is defined as where V C is the absorbed oil volume. In order to use the [26] transformation algorithms, the surface volume, V S , is subdivided into a thin part, V TN , and a thick part, V TK , as where Here, A TN and A TK are the areas occupied by the thick and thin surface slick volume, and T TN and T TK are the thicknesses of the thick and thin surface slicks.
To solve the advection-diffusion process (Equation (6)) and compute C S (Equation (7)), C D (Equation (7)), and C C (Equation (9)), the surface volume is broken into N constituent particles that are characterized by a particle volume index, v(n k ,t), particle status index, σ(n k ,t), and particle position vector x k (n k ,t) defined as where n k is the particle identification number. Following uncoupled Langevin equations, the time evolution of the particle position vector x k (n k ,t) is given as where the tensor A(x k ,t) represents the deterministic part (i.e., U in Equation (4)), while the tensor B(x k ,t) and ξ(t) characterizing random motion and a random motion factor represent the stochastic part of the flow field. Furthermore, the particle volume index, v(n k ,t), is subdivided into evaporative, v E (n k ,t), and non-evaporative, v NE (n k ,t), particle volume attributes as follows: The particle volume index, v(n k ,t), is updated using empirical formulas that relate to the time rate of change of oil slick volume state variables. The particle status index, σ(n k ,t), identifies the four particle classes (i.e., particles on the surface, subsurface or dispersed particles, sedimented particles, particles on the coasts). Further details of the expansion of Equation (10) and the description of the time rate of oil slick volume state variables are explained in De Dominicis et al. [17]. The above processes are summarized in a flowchart as shown in Figure 4.
The particle volume index, v(nk,t), is updated using empirical formulas that relate to the time rate of change of oil slick volume state variables. The particle status index, σ(nk,t), identifies the four particle classes (i.e., particles on the surface, subsurface or dispersed particles, sedimented particles, particles on the coasts). Further details of the expansion of Equation (10) and the description of the time rate of oil slick volume state variables are explained in De Dominicis et al. [17]. The above processes are summarized in a flowchart as shown in Figure 4.

Tidal Elevation and Current
To validate the model in the target area ( Figure 2e) prior to the construction of the LNG terminal, Table 2 lists the comparisons of the amplitudes and the phases of the tidal constituents between the computed and observed data [21]. Overall, the comparisons are in close accordance with absolute relative errors of 2.0% and 2.5% at T1 and T2, respectively. In Table 2, noticeable discrepancies of tidal amplitudes at T2, near the bottom of the reclaimed land (Figure 2e), may be due to the effects of surf-zone waves and associated intermittent breaking, which is not considered in the model. Breaking-induced splash-ups and droplets often cause spikes in the output signals, for example [27]. Nevertheless, Figure 5 shows that tidal elevations are generally in close match, even at T2, implying that the model is capable of predicting the total tidal elevation, its peaks and its temporal variations fairly accurately.  Table 2 lists the comparisons of the amplitudes and the phases of the tidal constituents between the computed and observed data [21]. Overall, the comparisons are in close accordance with absolute relative errors of 2.0% and 2.5% at T1 and T2, respectively. In Table 2, noticeable discrepancies of tidal amplitudes at T2, near the bottom of the reclaimed land (Figure 2e), may be due to the effects of surf-zone waves and associated intermittent breaking, which is not considered in the model. Breaking-induced splash-ups and droplets often cause spikes in the output signals, for example [27]. Nevertheless, Figure 5 shows that tidal elevations are generally in close match, even at T2, implying that the model is capable of predicting the total tidal elevation, its peaks and its temporal variations fairly accurately.  Tables 3 and 4 and Figure 6 show the comparisons of the computed and observed eastward (U) and northward (V) current velocity components for the major tidal constituents M2, S2, K1, and O1. For all the gauge stations, the amplitude approximation error, Eapp, i.e., − / × 100, where Aobs and Acomp are the observed and computed amplitudes, averaged over the four tidal constituents, was below 10 %, which indicates that, in general, the computed amplitudes and phases of the tidal currents did not deviate intensely from the observed data except for PC2.  Figure 6 show the comparisons of the computed and observed eastward (U) and northward (V) current velocity components for the major tidal constituents M 2 , S 2 , K 1 , and O 1 . For all the gauge stations, the amplitude approximation error, E app , i.e., A obs − A comp /A comp × 100, where A obs and A comp are the observed and computed amplitudes, averaged over the four tidal constituents, was below 10%, which indicates that, in general, the computed amplitudes and phases of the tidal currents did not deviate intensely from the observed data except for PC2.

Residual Current
To better investigate the hydrodynamic model's predictability of time-invariant tidal characteristics, the trajectories of the residual currents and the resultant average velocity magnitudes and directions were compared at the surface, middle and bottom layers (1,3,5) between the compared and observed data [21]. Residual currents are basically defined as the currents averaged over corresponding tidal periods [28]. Figures 7 and 8 depict the comparison of the trajectories and the resultant magnitudes and directions of the average velocities of the residual currents between the computed and observed data at PC2 and PC3, respectively. The comparison clearly shows that the EFDC model employed in the current study is successfully applicable in predicting residual currents not only on the surface but also at subsurface layers.

Water Circulation through Permeable Revetments
Port developments and reclamations are prevalent in South Korea due to the high population density. Therefore, revetment constructions are also in high demand around coastal regions, including the current Songdo LNG terminal area. Since the oil dispersion modeling in the current study needs accurate current prediction as an input in the presence of revetments, as depicted in the three scenarios (Figure 3), variation in the current fields in conjunction with revetment construction stages needs to be investigated. Since most of these revetments are permeable, even if the revetments are completely closed, the sea water circulates in and out of the embankments. Therefore, tidal elevation in areas surrounded by revetments also varies depending on tidal cycles and should be considered. To account for this permeability of the sea dike, Dupuit parabola [29] is employed as follows

Water Circulation through Permeable Revetments
Port developments and reclamations are prevalent in South Korea due to the high population density. Therefore, revetment constructions are also in high demand around coastal regions, including the current Songdo LNG terminal area. Since the oil dispersion modeling in the current study needs accurate current prediction as an input in the presence of revetments, as depicted in the three scenarios (Figure 3), variation in the current fields in conjunction with revetment construction stages needs to be investigated. Since most of these revetments are permeable, even if the revetments are completely closed, the sea water circulates in and out of the embankments. Therefore, tidal elevation in areas surrounded by revetments also varies depending on tidal cycles and should be considered. To account for this permeability of the sea dike, Dupuit parabola [29] is employed as follows where A C denotes a closed area surrounded by sea dikes, h 0 is the tidal elevation outside A, Q is the circulation flow rate per unit length of dike, and L is the total length of dike ( Figure 9). In LHS of Equation ( where AC denotes a closed area surrounded by sea dikes, h0 is the tidal elevation A, Q is the circulation flow rate per unit length of dike, and L is the total length ( Figure 9). In LHS of Equation (15, the time rate of change in h0 is used because the ability coefficient varies depending on the elevations both inside and outside of th ments. Data for the time rate of change in h0 were provided by Korean Ministry of and Fisheries.  Figure 10 shows the selected locations of the inside and outside of the western ment used for comparison. Three passages for an additional seawater exchange we sidered, as depicted in Figure 10, to reduce the possible instability of the revetment by uneven hydrostatic pressure on each side, especially under harsh weather con Table 5 shows the modeling results with and without the installation of the three pa When the additional passages were installed, the difference in the predicted tida tions at T3 and T4 are decreased approximately to half, indicating that the passage lation is effective in enhancing revetment stabilities.  Figure 10 shows the selected locations of the inside and outside of the western revetment used for comparison. Three passages for an additional seawater exchange were considered, as depicted in Figure 10, to reduce the possible instability of the revetment caused by uneven hydrostatic pressure on each side, especially under harsh weather conditions. Table 5 shows the modeling results with and without the installation of the three passages. When the additional passages were installed, the difference in the predicted tidal elevations at T3 and T4 are decreased approximately to half, indicating that the passage installation is effective in enhancing revetment stabilities.

Tidal current velocity fields
Considering the permeability of the sea dike, the current velocity fields were predicted through the EFDC model. The three initial conditions (i.e., case 1: without dike; case 2: with dike only; and case 3: with both dike and passage) were applied in the model. Figure 11a and 11b shows the current velocity fields in the enlarged view, as in Figure 2e, at the pre-construction stage before the sea dikes are built. The flood current fields are dominantly directed northeast, while the ebb current is generally directed in the opposite direction to southwest. Due to the combined effects of geometry and bathymetry, the primary current direction near the construction site is east-west in direction (Figure 11c and  10d). The maximum velocities in the detailed domain near construction site during the initial dispersion phase of the 12-h time frame were 0.12-0.59 m/s on the surface (layer 1). At layers 3 and 5, the maximum velocities decreased to 0.11-0.56 m/s and 0.09-0.49 m/s, respectively (current fields not shown for brevity in Figure 11).

Tidal Current Velocity Fields
Considering the permeability of the sea dike, the current velocity fields were predicted through the EFDC model. The three initial conditions (i.e., case 1: without dike; case 2: with dike only; and case 3: with both dike and passage) were applied in the model. Figure 11a,b shows the current velocity fields in the enlarged view, as in Figure 2e, at the pre-construction stage before the sea dikes are built. The flood current fields are dominantly directed northeast, while the ebb current is generally directed in the opposite direction to southwest. Due to the combined effects of geometry and bathymetry, the primary current direction near the construction site is east-west in direction (Figure 11c,d). The maximum velocities in the detailed domain near construction site during the initial dispersion phase of the 12-h time frame were 0.12-0.59 m/s on the surface (layer 1). At layers 3 and 5, the maximum velocities decreased to 0.11-0.56 m/s and 0.09-0.49 m/s, respectively (current fields not shown for brevity in Figure 11). Figure 12 shows the current fields assuming that only the revetments were constructed completely without allowing passages for additional direct seawater exchanges (case 2). The dominant directions of the tidal currents were more or less similar in direction, indicating that the effects of revetments on the major current direction are not crucial in the study area (Figures 11b and 12a). On the other hand, the maximum velocities near the construction site were significantly reduced to 0.04-0.27 m/s on the surface, as depicted by Figure 12c,d. Note that variations in the maximum values of tidal currents considering the permeable sea dikes are significant, especially at the northern end of the waterway, with and without considering the revetments.  Figure 12 shows the current fields assuming that only the revetments were constructed completely without allowing passages for additional direct seawater exchanges (case 2). The dominant directions of the tidal currents were more or less similar in direction, indicating that the effects of revetments on the major current direction are not crucial in the study area (Figure 11b and 12a). On the other hand, the maximum velocities near the construction site were significantly reduced to 0.04-0.27 m/s on the surface, as depicted by Figure 12c and 12d. Note that variations in the maximum values of tidal currents considering the permeable sea dikes are significant, especially at the northern end of the waterway, with and without considering the revetments.
In Figure 13, the current fields in case 3, additionally considering the effects of the three passages for seawater exchange on case 2, again shows that the resultant impact to the dominant current direction is minimal (Figure 13a,b). For the current velocity fields, the maximum velocities on the surface were slightly increased to 0.05-0.28 m/s compared with those in case 2. These calculated current fields considering seawater exchange via permeable revetments are directly used as inputs in predicting the oil spill based on MED-SLIK-II. In Figure 13, the current fields in case 3, additionally considering the effects of the three passages for seawater exchange on case 2, again shows that the resultant impact to the dominant current direction is minimal (Figure 13a,b). For the current velocity fields, the maximum velocities on the surface were slightly increased to 0.05-0.28 m/s compared with those in case 2. These calculated current fields considering seawater exchange via permeable revetments are directly used as inputs in predicting the oil spill based on MEDSLIK-II.

Oil Spill Dispersion
3.2.1. Scenario 1 Figure 14 shows the spatiotemporal variation of the current velocity fields in Scenari 1, where the total volume of the spilled oil was the largest among the three prescribe scenarios. Scenario 1 features an initial stage of revetment construction, with only limite parts of the revetments installed. Four hours after the oil spill incident was occurred, th oil slick had stretched and advected further offshore to the west away from the Wester revetment (Figure 14b). This is expected from the dominant current direction in the we near the construction site, as previously shown in Figure 11d. The oil slick reached th farthest location to the west after 6 h, at the end of the semidiurnal tidal cycle, and subse quently the advection direction is reversed toward the east (Figure 14c,d). Note that o thickness can be directly estimated from the oil concentration fields, based on a simp conversion. Assuming a homogeneous thickness over the spilled area, the oil thicknes can be calculated by dividing oil concentration by oil density. For example, with this con version relationship, the oil thickness in the spilled area of Figure 15c, where concentra tion is about 500 tonne/km 2 and LSFO density is 910 kg/m 3 , may be estimated as 0.55 mm

Oil Spill Dispersion
3.3.1. Scenario 1 Figure 14 shows the spatiotemporal variation of the current velocity fields in Scenario 1, where the total volume of the spilled oil was the largest among the three prescribed scenarios. Scenario 1 features an initial stage of revetment construction, with only limited parts of the revetments installed. Four hours after the oil spill incident was occurred, the oil slick had stretched and advected further offshore to the west away from the Western revetment (Figure 14b). This is expected from the dominant current direction in the west near the construction site, as previously shown in Figure 11d. The oil slick reached the farthest location to the west after 6 h, at the end of the semidiurnal tidal cycle, and subsequently the advection direction is reversed toward the east (Figure 14c,d). Note that oil thickness can be directly estimated from the oil concentration fields, based on a simple conversion. Assuming a homogeneous thickness over the spilled area, the oil thickness can be calculated by dividing oil concentration by oil density. For example, with this conversion relationship, the oil thickness in the spilled area of Figure 15c, where concentration is about 500 tonne/km 2 and LSFO density is 910 kg/m 3 , may be estimated as 0.55 mm, largely comparable to typical values of oil slick thickness. Approximately 10 h after the initiation of the spill, the oil slick entered the expected reclamation area (Figure 14e). In 12 h, the oil slick is advected to the inner reclamation area and dispersed (Figure 14f). At the same time, oil slick with relatively lower oil concentration tends to be advected north following the alongshore current observed in Figure 11a. Note that the permeable revetment only allows a limited oil slick near the north edge of the western revetment (Figure 14f).

Scenario 2
Unlike Scenario 1, in Scenario 2, the initiation of the oil spill occurs near the northern part of the waterway, as shown in Figure 15. In this case, the oil slick is initially attached along the coasts within the waterway (Figure 15a). As the tidal cycle approaches the maximum ebb, the oil slick is stretched off the waterway to the north and further advected to the southwest (Figure 15c) following the dominant tidal currents, as depicted in Figures 11b and 12b. The oil slick is advected more to the west until the reversal in the direction of the tidal currents occurred (Figure 15c). Then, due to the relatively weak permeability of the revetment, the oil slick is concentrated to form a more circular shape (Figure 15d). Moreover, the oil slick is not advected inside the reclamation area until 12 h after the oil spill initiation. After 8 and 12 h, the oil slick is stretched and attached along the coasts to cover approximately 5.2 and 9.8 km in length, respectively.
Samaras et al. [30] modified the MEDSLIK-II to better reflect the local coast type and beach properties, which is one of the most important many interrelated environmental drivers in predicting the oil behavior at the shoreline [31]. Specifically, Samaras et al. [30] revised the release criteria for beached particles, proposed a more detailed shoreline discretization, and adopted the Oil Holding approach to improve the prediction of C C (Equation (9)) in MEDSLIK-II. Samaras et al. [30] explicitly showed that the oiled shoreline modeled with the modified version better matched satellite observation. With these modifications applied, the current model has the potential to enhance the predictability of oil behavior near the coastline in the future works. The accurate prediction of the expected length of the oiled shoreline is important to evaluate the clean-up costs correctly. For example, Etkin [32] reported that about 6000 USD/tonne are needed to clean-up an oiled shoreline of 2-5 km, while about 11,000 USD/tonne are needed for the recovery of an oiled shoreline of 8-15 km.

Scenario 3
In Scenario 3, the oil spill originated at a location close to that in Scenario 2, near the northern part of the waterway. In 4 h, the oil slick is advected to the southwest with a sharp turnover immediately off the waterway, mainly due to the maximum ebb current (Figure 16b). Here, the shape of the oil slick is similar to that in Scenario 2, except for the oil concentration. The oil concentration during the maximum ebb in Scenario 3 is generally lower because the total leakage volume is less than half of that in Scenario 2. Moreover, comparison of Figures 15f and 16f showed that the oil slick only entered inside the expected reclamation area in scenario 2, where the revetment construction was partially completed. This suggests that an emergent oil-fence installation may be required in scenario 2 to prevent possible hazards caused by the oil penetrating during the further reclamation process. Nevertheless, the overall tendency of advection, diffusion, and transformation of the surface oil is very similar during the entire duration of 12 h ( Figure 16).

Conclusions
In the current study, an oil spill modeling system was proposed by coupling the EFDC for computing tidal elevation and current fields, which were inputted to MEDSLIK-II [17] to predict oil slick advection-diffusion and transformation processes. Wind velocity fields, computed based on extreme value analysis utilizing the probability density func

Conclusions
In the current study, an oil spill modeling system was proposed by coupling the EFDC for computing tidal elevation and current fields, which were inputted to MEDSLIK-II [17] to predict oil slick advection-diffusion and transformation processes. Wind velocity fields, computed based on extreme value analysis utilizing the probability density function of a 5-year easterly wind, were also fed into MEDSLIK-II, accounting for the atmospheric condition. The current oil spill modeling system was applied to the three ship collision scenarios generated based on the proposed construction plan of the Songdo LNG terminal on reclaimed land in South Korea. The numbers, locations, types, and the sizes of the construction vessels in the three scenarios were varied to closely reflect the three distinguished construction phases (phase 1, 2, and 3) depending on the progressive time frames of construction. Prior to the MEDSLIK-II simulation, the current fields and the residual current (i.e., current field averaged over a tidal cycle) were successfully verified based on the observation data. The effects of permeable revetments, required for reclamation, on the currents were also investigated. For each scenario, the simulated results showed distinct patterns in the advection, dispersion, and transformation of the oil slick. It can be deduced from the simulated results that the oil absorption into the coast, which causes immense damage to coastal communities, is highly dependent on the tidal currents, volume of oil spilled, and nearby construction activities.
The current modeling results have a potential to be further refined to better reflect the local beach properties by revising the release criteria for beached particles and by discretizing shorelines in more detail based on the categorized coast types. An indoor laboratory experiment is also required to precisely identify the weathering process of spilled oil and thereby determine the parameters involved in the modelling procedure. In this regard, a depth-varying oceanographic flow condition can be practically reproduced by adopting the state-of-art wavemaker theory and bottom treatment [33,34]. Through accurate prediction of the oiled shoreline, a hazard mitigation plan, such as an emergent oil-fence installation, can be effectively incorporated into the modeling results to minimize the environmental impact of oil spill accidents. Consequently, the proposed modeling system turns out to successfully provide an operational tool for simulating likely schemes of oil spill accidents within the construction time frame.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.