Estimation of the Mechanical Recovery Potential of Spilled Oil at Sea Considering the Spatial Thickness Distribution

: Recovery modeling and countermeasures for oil spilled at sea have been extensively researched, but research remains insu ﬃ cient on recovery potential estimation methods. It is required to access the mechanical recovery potential by considering the relationship between oil behavior, environmental conditions, and the performance of clean-up activities. Two response-planning models were developed in this study. One is a spatially uniform recovery model for estimating recovery potential that reﬂects weathering, oil properties, and equipment e ﬃ ciency. The other is a spatially nonuniform recovery model that considers not only the above characteristics but also local thickness reduction by skimming. A comparison between the two models and an analysis of their e ﬀ ects on response was carried out through the calculation using an accident scenario. It is possible to analyze the e ﬀ ect of the thin slicks, natural dissipation, and the quantiﬁcation of deployable skimming systems with the spatially nonuniform recovery model. Finally, we analyzed interrelationships among residual oil volume on the sea, response time, and the number of skimming systems.


Introduction
Global oil consumption has increased [1,2]. Although the frequency of oil spill accidents at sea has decreased, large and small accidents still occur [3][4][5]. Thus, as long as oil transport and consumption are maintained or increased, we must continue to be prepared for and respond to catastrophic oil spill accidents. Researches pertaining to oil spills at sea can include remote sensing, trajectory modeling, spill modeling, and countermeasures [6][7][8][9]. The present study focuses on recovery at sea through countermeasure planning. Numerous studies have been conducted on spilled oil properties [10][11][12][13][14] and weathering [15][16][17][18][19][20][21]. Mackay and Matsugu [16] researched the evaporation rate of spilled oil, and Delvigne and Sweeney [19] investigated natural dispersion. These studies associated the sea state and the properties of the spilled oil with the weathering process. The properties of spilled oil undergo continuous change with time due to weathering processes such as spreading, evaporation, emulsification, and dispersion [22][23][24]. Mackay et al. [25] developed an oil spill behavior model that encompasses weathering and changes in oil properties, and Berry et al. [26] modeled oil transport and fate processes. However, these studies did not suggest the prediction of the response potential.

Methods for Estimating Recovery Potential of Skimming System
The recovery capacity calculation method considering oil property and equipment efficiency change was implemented in both the SUR model and SNR model. Both models of the present study consisted of weathering calculation and recovery calculation. The recovery calculation made a change in sea surface oil volume, and then the volume change made an effect in weathering calculation, and vice versa. The main difference between the two models is that the SNR model distinguishes the skimmed and the unskimmed areas in the oil slick and calculates the effect of recovery on thickness variation of each zone separately.

Spatially Uniform Recovery (SUR) Model
To calculate the recovery potential, the encounter rate (ER) of the emulsion, which is a mixture of oil and water, for a single skimming system could be obtained via Equation (1) [44,46,48]: where w is the boom swath, v is the tow speed, and h em is the thickness of the emulsion. The amount of recoverable encountered emulsion depended on the oil and sea conditions. The amount of recovered emulsion per unit time, emulsion recovery rate (ERR), is determined by multiplying the throughput efficiency (TE) [9,44]. In addition, the oil recovery rate (ORR) could be determined by subtracting the amount of water from the amount of emulsion [43]: where ERR is the emulsion recovery rate, TE is the volume ratio of the recovered emulsion to the encountered emulsion, ORR is the oil recovery rate, and Φ is the water fraction in the emulsion. The SUR model considered oil type and sea conditions and calculated the slick area and weathering over time. The area of the oil slick is expressed as the product of characteristic lengths l1 and l2: where l1 is the spreading length under calm sea conditions, and l2 is the length increase due to spreading and wind [49]. The geometry of slicks was determined by the physical characteristics of oil films and by environmental parameters. Ermakov et al. [50] developed a model of the spreading of spills accounting for the surface stresses induced by wind waves. In this study, a simple spreading calculation method was used. In a calm sea, spreading undergoes three steps: an initial step, in which gravity and inertia are important (Equation (5)); an intermediate step dominated by gravity and viscosity (Equation (6)); and a final step during which surface tension and viscosity balance each other [51].
where t is time in seconds; a 3600 s time step is used in the calculation. t 0 is the transition time from Equation (5) to Equation (6) and is a function of spilled volume, density, and viscosity [51]. After calculating the spreading on a calm sea, the transport due to wind was considered assuming a constant wind direction during the time step. Movement by the wind was assumed to occur at 3% of the wind speed [52]: where Z(t) is the distance that the oil slick has traveled due to wind at time t.
Oil thickness was calculated every hour using Equation (8). The emulsion thickness could be obtained by inserting (1 − Φ) in the oil thickness equation: where h oil is the average oil thickness, V is the remaining volume of oil on the sea surface, and A mod is the slick area. Area modification, A mod , considered a statistical estimate of the relative thickness distribution along the downwind axis with the droplet spreading algorithm [47]. The remaining oil volume over time is defined as the initial spilled volume minus the naturally removed volume and mechanically recovered volume from the start of the spill until time t: where V 0 is the instantaneous initial oil spill volume, and W is the evaporated and naturally dispersed volume at time t.

Spatially Nonuniform Recovery (SNR) Model
The spatial and temporal modifications were applied in the SNR model to take account of the skimmed zone and local oil thickness reduction by mechanical recovery. Applied methods and expected effects are described below. Concrete results using an accident scenario are detailed in Section 3.

Spatial Improvements
When calculating the amount of recovered oil by the skimming systems, studies of the SUR model assumed that the mechanically recovered volume was a reduction in the total volume of the spilled oil slick, as shown in Figure 1a. Although the relative thickness distribution by wind and droplet was already considered with statistical estimation [47], it was not easy to differentiate the local reduction of thickness by skimming activities. The thickness reduction by mechanical skimming showed at the same rate across the whole oil slick because SUR model studies assumed that oil was recovered at the same rate in all areas. However, as shown in Figure 1b, oil was mechanically recovered from only specific parts of the spill area rather than over the entire oil slick. In other words, the thickness of oil slick was decreased only where the skimming systems were operated. The actual oil slick thickness changed at different rates in skimmed areas and other areas. Therefore, the spatial variation in oil slick thickness should be considered, as the recovery potential of the skimming systems was affected by the oil slick thickness, as shown in Equation (1). where hoil is the average oil thickness, V is the remaining volume of oil on the sea surface, and Amod is the slick area. Area modification, Amod , considered a statistical estimate of the relative thickness distribution along the downwind axis with the droplet spreading algorithm [47]. The remaining oil volume over time is defined as the initial spilled volume minus the naturally removed volume and mechanically recovered volume from the start of the spill until time t: where V0 is the instantaneous initial oil spill volume, and W is the evaporated and naturally dispersed volume at time t.

Spatially Nonuniform Recovery (SNR) Model
The spatial and temporal modifications were applied in the SNR model to take account of the skimmed zone and local oil thickness reduction by mechanical recovery. Applied methods and expected effects are described below. Concrete results using an accident scenario are detailed in Section 3.

Spatial improvements
When calculating the amount of recovered oil by the skimming systems, studies of the SUR model assumed that the mechanically recovered volume was a reduction in the total volume of the spilled oil slick, as shown in Figure 1a. Although the relative thickness distribution by wind and droplet was already considered with statistical estimation [47], it was not easy to differentiate the local reduction of thickness by skimming activities. The thickness reduction by mechanical skimming showed at the same rate across the whole oil slick because SUR model studies assumed that oil was recovered at the same rate in all areas. However, as shown in Figure 1b, oil was mechanically recovered from only specific parts of the spill area rather than over the entire oil slick. In other words, the thickness of oil slick was decreased only where the skimming systems were operated. The actual oil slick thickness changed at different rates in skimmed areas and other areas. Therefore, the spatial variation in oil slick thickness should be considered, as the recovery potential of the skimming systems was affected by the oil slick thickness, as shown in Eq. (1).

Figure 1.
Schematic drawing on the correlation between oil slick behavior and the skimming zone in (a) the spatially uniform recovery (SUR) model and (b) the spatially nonuniform recovery (SNR) model. The yellow area represents a zone where skimming is in progress. Aun: area of the unskimmed zone; As: area of the skimmed zone. The legend shows the amount of oil remaining on the sea surface, which can also be thought of as the oil thickness variation. The yellow area represents a zone where skimming is in progress. Aun: area of the unskimmed zone; As: area of the skimmed zone. The legend shows the amount of oil remaining on the sea surface, which can also be thought of as the oil thickness variation.
Furthermore, the SUR model was limited in that it considered neither the space occupied by an individual skimming system in the spatial scope of oil slick nor the interference between multiple skimming systems. Thus, when calculating the response capacity using a SUR model, the number of response resources applied to the oil slick area was virtually unlimited; the SUR model was allowed to apply the infinite quantity of skimming systems. However, there is a practical limit to the number of skimming systems that can be deployed to the limited size of an oil slick. Two spatial improvements were applied to the SNR model to overcome these two limitations.
First, to distinguish spatial differences in thickness, we differentiated the skimmed and unskimmed zones herein; A un denotes the area of the unskimmed zone and A s denotes the area of the skimmed zone. Figure 1a,b illustrate mechanisms from the SUR model and the SNR model, respectively. As shown in Figure 1a, the skimmed zone was not distinguished in the SUR model. The SNR model, however, distinguished both zones, as shown in Figure 1b. The yellow section denotes a zone in which skimming is in progress at time t. The skimming zone at time t transformed into a skimmed zone at next time, t + ∆t. In other words, an area recovered by one skimming system became a skimmed zone in the next time step. The unrecovered oil among the encounter rate where the skimming system passed was regarded as the remaining oil of the skimmed zone. The amount of remaining oil in a skimmed zone equals ER minus ERR in Equation (2), by the definition of TE.
The oil thickness in the skimmed zone could be represented differently from that in the unskimmed zone because each zone was calculated separately. Furthermore, oil properties and behaviors such as weathering processes, including spreading phenomena, were calculated separately in the unskimmed and skimmed zones; it was assumed that any oil recovery during the next time step occurs only in an unskimmed zone.
Second, to consider the space occupied by an individual skimming system and interference between skimming systems, the area recovered by one skimming system and the quantity of skimming systems that can be deployed in the oil slick were computed and applied in the improved model.
Skimmers and oil booms are usually operated counter to the wind direction to take advantage of the oil movement due to wind during recovery [6]. The SNR model was thus developed to recover oil in the l2 direction. A s , the area skimmed by a single skimming system in ∆t, was calculated similarly to the volume rate in Equation (1), as shown in Equation (11). As the mechanical recovery was assumed to occur in the l2 direction, the l2 of the skimmed zone was identical to the l2 of the entire oil slick. The l1 length (l1 s ) of the skimmed zone was calculated by dividing the skimmed area A s by l2: where l1 s is the l1 length of the skimmed zone. A spatial margin was added to prevent the collision of skimming systems, as shown in Equation (12). After this change, the area occupied by a single skimming system could be defined by Equation (13): where A occu is the area occupied by a single skimming system and A margin is the margin area. The full A occu did not represent the recovered area but was instead used to consider the interactions between skimming systems. The recovered area A s is equal to A occu minus A margin . l1 margin is the marginal length in the l1 direction, which is assumed to be half of l1 s . Finally, the maximum number of skimming systems, x(t), that could be deployed in the unskimmed oil slick area at time t could be obtained via: After oil recovery at time t, the l1 length of the unskimmed zone was updated via: l1 un (t) = l1 un (t) − l1 s (t) (15) where l1' un (t) is the l1 length of the unskimmed zone after skimming, and l1 un is the l1 length of the unskimmed zone before recovery. If there was no skimming, l1 s is zero and l1' un is equal to l1 un . Each space is illustrated in Figure 2.
In spatially uniform recovery studies, differences in thickness due to mechanical recovery by skimming could not be resolved. The first improvement of the SNR model allowed us to distinguish skimmed and unskimmed zones, which in turn facilitates the individual calculation of the area of and oil thickness in each zone. As a result, weathering could be calculated differently for each zone, as well. where l1'un(t) is the l1 length of the unskimmed zone after skimming, and l1un is the l1 length of the unskimmed zone before recovery. If there was no skimming, l1s is zero and l1'un is equal to l1un. Each space is illustrated in Figure 2.
In spatially uniform recovery studies, differences in thickness due to mechanical recovery by skimming could not be resolved. The first improvement of the SNR model allowed us to distinguish skimmed and unskimmed zones, which in turn facilitates the individual calculation of the area of and oil thickness in each zone. As a result, weathering could be calculated differently for each zone, as well. The thickness difference between the SUR model and the SNR model is illustrated in Figure 3. The SNR model supposed that encounter areas where the skimming systems passed were recovered only one time. The thickness in the SUR model tended to be slightly thinner than that of the unskimmed zone in the SNR model. A total of two zones ('unskimmed' and 'skimmed') existed at 2 h in the SNR model because one skimming system began to recover at 1 h. The recovered zone is thinner than the unskimmed zone. This zone could be computed separately from the unskimmed zone due to the improvement mentioned below.
The second improvement allowed the SNR model to describe the space occupancy by skimming systems and prevent interaction between them. The deployable quantity of skimming systems could also be calculated over time. In other words, the number of skimming systems that can be deployed was limited by considering interference between skimming systems. Figure 4 shows the number of skimming systems used in the calculations and the maximum number of skimming systems, x(t), that can be deployable in the unskimmed oil slick of the SNR model. In both calculations, 30 skimming systems were used as an input value. This value was applied continuously in the calculation of the SUR model. The quantity of skimming systems was 30 until 9 h in the SNR model. However, this value decreased from 10 h in the SNR model as the maximum quantity of applicable skimming systems decreased. The SNR model could consider the space occupied by an individual skimming system and the interactions between multiple skimming systems. These improvements allowed the model to reflect actual conditions better. The thickness difference between the SUR model and the SNR model is illustrated in Figure 3. The SNR model supposed that encounter areas where the skimming systems passed were recovered only one time. The thickness in the SUR model tended to be slightly thinner than that of the unskimmed zone in the SNR model. A total of two zones ('unskimmed' and 'skimmed') existed at 2 h in the SNR model because one skimming system began to recover at 1 h. The recovered zone is thinner than the unskimmed zone. This zone could be computed separately from the unskimmed zone due to the improvement mentioned below.
The second improvement allowed the SNR model to describe the space occupancy by skimming systems and prevent interaction between them. The deployable quantity of skimming systems could also be calculated over time. In other words, the number of skimming systems that can be deployed was limited by considering interference between skimming systems. Figure 4 shows the number of skimming systems used in the calculations and the maximum number of skimming systems, x(t), that can be deployable in the unskimmed oil slick of the SNR model. In both calculations, 30 skimming systems were used as an input value. This value was applied continuously in the calculation of the SUR model. The quantity of skimming systems was 30 until 9 h in the SNR model. However, this value decreased from 10 h in the SNR model as the maximum quantity of applicable skimming systems decreased. The SNR model could consider the space occupied by an individual skimming system and the interactions between multiple skimming systems. These improvements allowed the model to reflect actual conditions better.

Temporal Improvement
Oil spreading is a function of the remaining oil volume and time, and calculations of oil properties, weathering, and mechanical recovery were carried out at every time-step. The oil volume changed continuously with time due to weathering and clean-up activities, including skimming recovery. Weathering is a function of slick area (i.e., weathering may increase or decrease with the area). However, the SUR model does not explicitly reflect changes in oil volume in calculations of the spreading area, as shown in Equations (5) and (6).
The SNR model uses a remaining volume that changes continuously, instead of the initial spill volume, to reflect the changes in slick volume over time when calculating the spreading. Equations (5) and (6) in the SUR model thus have been modified to Equations (16) and (17) in the SNR model: where V(t) is changed by weathering processes, such as oil evaporation, dispersion, and mechanical recovery.
Next, the area change over time should be calculated separately in each zone because the spatial improvement distinguished the skimmed zone and unskimmed zone. The l1 of the unskimmed zone is updated to l1' un as in Equation (15). The l1 un (t + ∆t) is defined by adding the increase by spreading and wind, ∆l1, to the length l1' un (t): where ∆l1 un (∆t) is the difference between the two l1 un values calculated at time t and t + ∆t in Equation (16) or Equation (17) and can be interpreted as the length increase due to spreading. Likewise, the l1 of the skimmed zone at t + ∆t could be calculated as the sum of ∆l1 s and the l1 s (t) at the previous time step. The l1 s (t) value is obtained by dividing the skimmed area A s by l2 as in Equation (11), and the spreading at this value is: The change in the length of l1 during ∆t (from t to t + ∆t), ∆l1 tot (∆t), is equal to the increase in length due to the total spreading in the unskimmed and skimmed zones: The l2 lengths of the unskimmed and skimmed zones could be calculated using the same equation. As shown in Equation (7), l2 is defined by considering the length added due to wind to l1, which is the length considering spreading. Therefore, the spreading length ∆l1 in ∆t is added to the length at the previous time step. The increase in length due to wind, ∆Z, is also added: An area comparison between the two models is shown in Figure 5. The area of the unskimmed zone, assumed to be where recovery occurred in the SNR model, is smaller than the SUR model area. Thus, the thickness of the SNR model is thicker than that of the SUR model. Furthermore, it is possible to calculate the areas of the skimmed and unskimmed zones separately described above, using the volume of remaining oil instead of the initial spill amount. As a result, the area changes in the skimmed and unskimmed zones showed different behaviors. To separately calculate the skimmed and unskimmed zone areas, the unskimmed zone is defined as the outside the working range of the skimming systems. The spreading of the unskimmed zone during the next time step is calculated in this area. Therefore, if the skimmed area to be excluded becomes more significant than the area increment due to spreading at time t, the area of the unskimmed zone at t + ∆t decreases.
On the other hand, if the area increased due to spreading is larger, the area of the unskimmed zone will increase.

Comparison of Calculation Results Using a Hypothetical Accident Scenario
In this section, the calculation results of two models (SUR model and SNR model) are compared and analyzed. The hypothetical accident scenario was created. The inputs used for the calculations are listed in Table 1. In this scenario, 500 kl of Bunker C oil was assumed to be spilled in the offshore area of Busan, Republic of Korea. The simulation time was 24 hours after the accident. All 20 skimming systems were deployed from 1 h to 24 h. The three-year average wind velocity and water temperature in the Busan area [53] between 2014 and 2016 were used. Information of the Komara-50 skimmer, which is installed on the response vessel of the Busan Coast Guard, was used in the simulation. A TE was assumed to be 75%. Recovery efficiency (RE), which is defined as the percent of emulsion excluding water in the total recovered volume, was calculated over time. The size of the storage tank was assumed to be large enough to exclude time for emptying the storage tank.
The SUR model considered reduction by both mechanical recovery and natural weathering, which includes evaporation and dispersion. The SNR model not only considered removal by recovery and weathering but also assumed that thin slicks (with a thickness of thinner than 0.6 μm) had been removed, as such slicks dissipate naturally from a response viewpoint [7,54].
As shown in Figure 6, the SNR model showed a total removal of 452 kl of oil during 24 h by three mechanisms, namely mechanical recovery, weathering, and sheen out due to thin thickness. Of this, 54% is due to recovery, 7% is due to weathering, and 39% is due to dissipation. There was no oil dissipation at 12 h. This could be explained by the thin skimmed zone generated from 17 h onward, which has a thickness of less than 0.6 μm. For example, 80 skimmed zones had a thickness of less 0.6 μm at 17 h, and the total volume of these 80 zones was approximately 29.7 kl. The total oil volume removed by dissipation during 24 h was 178 kl with 347 skimmed zones, and 20 skimmed zones were left. The same criteria were applied to the SUR model. However, no reduction in thickness to less than 0.6 μm was observed, as oil decreased over the whole oil slick body rather than in specific areas. In other words, the SUR model considered the thickness variation not by the mechanical recovery Area difference between SUR model and SNR model

Comparison of Calculation Results Using a Hypothetical Accident Scenario
In this section, the calculation results of two models (SUR model and SNR model) are compared and analyzed. The hypothetical accident scenario was created. The inputs used for the calculations are listed in Table 1. In this scenario, 500 kl of Bunker C oil was assumed to be spilled in the offshore area of Busan, Republic of Korea. The simulation time was 24 h after the accident. All 20 skimming systems were deployed from 1 h to 24 h. The three-year average wind velocity and water temperature in the Busan area [53] between 2014 and 2016 were used. Information of the Komara-50 skimmer, which is installed on the response vessel of the Busan Coast Guard, was used in the simulation. A TE was assumed to be 75%. Recovery efficiency (RE), which is defined as the percent of emulsion excluding water in the total recovered volume, was calculated over time. The size of the storage tank was assumed to be large enough to exclude time for emptying the storage tank. The SUR model considered reduction by both mechanical recovery and natural weathering, which includes evaporation and dispersion. The SNR model not only considered removal by recovery and weathering but also assumed that thin slicks (with a thickness of thinner than 0.6 µm) had been removed, as such slicks dissipate naturally from a response viewpoint [7,54].
As shown in Figure 6, the SNR model showed a total removal of 452 kl of oil during 24 h by three mechanisms, namely mechanical recovery, weathering, and sheen out due to thin thickness. Of this, 54% is due to recovery, 7% is due to weathering, and 39% is due to dissipation. There was no oil dissipation at 12 h. This could be explained by the thin skimmed zone generated from 17 h onward, which has a thickness of less than 0.6 µm. For example, 80 skimmed zones had a thickness of less 0.6 µm at 17 h, and the total volume of these 80 zones was approximately 29.7 kl. The total oil volume removed by dissipation during 24 h was 178 kl with 347 skimmed zones, and 20 skimmed zones were left. The same criteria were applied to the SUR model. However, no reduction in thickness to less than 0.6 µm was observed, as oil decreased over the whole oil slick body rather than in specific areas. In other words, the SUR model considered the thickness variation not by the mechanical recovery but by natural spreading. The SUR model showed a total removal of 251 kl of oil during 24 h due to two mechanisms. Of this, 230 kl was due to recovery, 21 kl was due to weathering.

Oil spill scenario Input values
Spilled amount (kl) 500 Oil    Figure 7 shows the number of skimming systems used in the SNR model calculation and the quantity of skimmed zones over time. Twenty skimming systems began response 1 h after the accident. Hence, 20 skimmed zones and one unskimmed zone existed at 2 h. Likewise, 20 skimmed zones appeared again due to recovery activities at 2 h, leaving 40 skimmed zones and one unskimmed zone at 3 h. In this way, the zones recovered at time t were converted to skimmed zones at t + ∆t and calculated separately from the unskimmed zone. Twenty skimmed zones were newly generated at 17 h, because 20 skimming systems were operational at 16 h. At the same time, however, 80 of the existing skimmed zones dissipated due to the reduction of the thickness (i.e., <0.6 µm). Therefore, the number of skimmed zones decreased by 60 at 17 h compared to that at 16 h. The thickness of the unskimmed zone was reduced to less than 0.6 µm at 21 h. As a result, no more skimmed zones were generated. Though the applied number of skimming systems was initially constant at 20, it then dropped below 20 at 17 h, as the number of deployable skimming systems decreased with the unskimmed area ( Figure 8).
The size of the unskimmed zone (where recovery units can be operated) could be grasped in the SNR model. The SNR model took into account the space occupied by skimming systems and distinguished unskimmed and skimmed zones to consider spatial changes in thickness. Therefore, the maximum number of deployable skimming systems could also be calculated at a given time.
In contrast, the SUR model did not explore the maximum number of deployable skimming systems, so 20 skimming systems were applied continuously in all calculations; spatial differences could not be distinguished, either. existing skimmed zones dissipated due to the reduction of the thickness (i.e., < 0.6 μm). Therefore, the number of skimmed zones decreased by 60 at 17 h compared to that at 16 h. The thickness of the unskimmed zone was reduced to less than 0.6 μm at 21 h. As a result, no more skimmed zones were generated. Though the applied number of skimming systems was initially constant at 20, it then dropped below 20 at 17 h, as the number of deployable skimming systems decreased with the unskimmed area ( Figure 8). The size of the unskimmed zone (where recovery units can be operated) could be grasped in the SNR model. The SNR model took into account the space occupied by skimming systems and distinguished unskimmed and skimmed zones to consider spatial changes in thickness. Therefore, the maximum number of deployable skimming systems could also be calculated at a given time. In contrast, the SUR model did not explore the maximum number of deployable skimming systems, so 20 skimming systems were applied continuously in all calculations; spatial differences could not be distinguished, either. Figure 8 illustrates the slick areas of the SUR model and the unskimmed zones in the SNR model. In the SUR model, the slick area was calculated as a function of the initial spill volume, as in Equations (5) and (6), and it continuously increased without reflecting the volume changes by the mechanical recovery. In contrast, it was possible to reflect volume changes since the SNR model distinguished skimmed and unskimmed zones, as described above. As the volume continuously decreased and the unskimmed zone was changed to skimmed zones beginning at 2 h by oil recovery, the unskimmed area was smaller than the SUR model area. Furthermore, area decrease could be observed in the unskimmed area from 6 h to 20 h. It was because the area converted from unskimmed to skimmed was larger than the area increases of the unskimmed zone due to spreading during this time.  The oil thickness of the unskimmed zone in the SNR model can be seen in Figure 9. In the SUR model, the volume decreased due to oil recovery in all oil slick areas. In contrast, the SNR model distinguished skimmed zones that became thin due to oil recovery from the unskimmed zone. Therefore, the thickness of the unskimmed zone was not directly affected by the mechanical removal in the SNR model. If no recovery activity was undertaken, any zone could remain unconverted into the skimmed zone. Then l1s in Eq. (15) became zero, and the difference between l1 values from the SUR model and l1un from the SNR model appeared by only Equations (16) and (17). The volume reduction over time was reflected in these equations. The l1un and area in the SNR model were smaller than the l1 and area in the SUR model. The thickness of the unskimmed zone in the SNR model became thicker than that in the SUR model. However, the oil remaining in the unskimmed zone and  Figure 8 illustrates the slick areas of the SUR model and the unskimmed zones in the SNR model. In the SUR model, the slick area was calculated as a function of the initial spill volume, as in Equations (5) and (6), and it continuously increased without reflecting the volume changes by the mechanical recovery. In contrast, it was possible to reflect volume changes since the SNR model distinguished skimmed and unskimmed zones, as described above. As the volume continuously decreased and the unskimmed zone was changed to skimmed zones beginning at 2 h by oil recovery, the unskimmed area was smaller than the SUR model area. Furthermore, area decrease could be observed in the unskimmed area from 6 h to 20 h. It was because the area converted from unskimmed to skimmed was larger than the area increases of the unskimmed zone due to spreading during this time.
The oil thickness of the unskimmed zone in the SNR model can be seen in Figure 9. In the SUR model, the volume decreased due to oil recovery in all oil slick areas. In contrast, the SNR model distinguished skimmed zones that became thin due to oil recovery from the unskimmed zone. Therefore, the thickness of the unskimmed zone was not directly affected by the mechanical removal in the SNR model. If no recovery activity was undertaken, any zone could remain unconverted into the skimmed zone. Then l1 s in Equation (15) became zero, and the difference between l1 values from the SUR model and l1 un from the SNR model appeared by only Equations (16) and (17). The volume reduction over time was reflected in these equations. The l1 un and area in the SNR model were smaller than the l1 and area in the SUR model. The thickness of the unskimmed zone in the SNR model became thicker than that in the SUR model. However, the oil remaining in the unskimmed zone and the thickness was decreased, as mechanical recovery progressed and the unskimmed zone was converted into skimmed zones. At 21 h, the unskimmed zone dissipated as its thickness decreased below 0.6 µm, and weathering and recovery of that were no longer calculated.  Figure 10 shows oil recovery rate and cumulative recovery volume. The recovery rate featured a trend similar to thickness, and the oil recovery rate was more significant in the SNR model than in the SUR model until 14 h. However, the recovery rate of the SNR model was lower than that of the SUR model at 15 h result from the reduced thickness in the SNR model. The amount of oil recovered during 24 h was 230 kl in the SUR model and 243 kl in the SNR model. Recovery volume of the SNR model was larger by 24 kl than that of the SUR model, despite the smaller quantity of skimming systems applied from 17 h onward.  Figure 10 shows oil recovery rate and cumulative recovery volume. The recovery rate featured a trend similar to thickness, and the oil recovery rate was more significant in the SNR model than in the SUR model until 14 h. However, the recovery rate of the SNR model was lower than that of the SUR model at 15 h result from the reduced thickness in the SNR model. The amount of oil recovered during 24 h was 230 kl in the SUR model and 243 kl in the SNR model. Recovery volume of the SNR model was larger by 24 kl than that of the SUR model, despite the smaller quantity of skimming systems applied from 17 h onward.
The remaining oil volume on the sea surface over 24 h is illustrated in Figure 11. The total remaining volume of the SNR model was equal to the sum of unskimmed and all skimmed zones with a thickness of 0.6 µm or larger. The amount of oil remaining (compared to the initial spilled volume of 500 kl) was 50% in the SUR model and 10% in the SNR model at 24 h. However, the recovered volume of the SNR model was calculated using less than 20 skimming systems from 17 h onward. To compare the recovery time rather than volume, half of the initial spilled volume, 250 kl, is indicated in Figure 11. It took 24 h and 14 h for the volume of the remaining oil in the SUR model and the SNR model respectively to reach half of the initially spilled volume. Thus, less time was required to remove half of the initial spilled volume in the SNR model. Figure 10 shows oil recovery rate and cumulative recovery volume. The recovery rate featured a trend similar to thickness, and the oil recovery rate was more significant in the SNR model than in the SUR model until 14 h. However, the recovery rate of the SNR model was lower than that of the SUR model at 15 h result from the reduced thickness in the SNR model. The amount of oil recovered during 24 h was 230 kl in the SUR model and 243 kl in the SNR model. Recovery volume of the SNR model was larger by 24 kl than that of the SUR model, despite the smaller quantity of skimming systems applied from 17 h onward.  The remaining oil volume on the sea surface over 24 h is illustrated in Figure 11. The total remaining volume of the SNR model was equal to the sum of unskimmed and all skimmed zones with a thickness of 0.6 μm or larger. The amount of oil remaining (compared to the initial spilled volume of 500 kl) was 50% in the SUR model and 10% in the SNR model at 24 h. However, the recovered volume of the SNR model was calculated using less than 20 skimming systems from 17 h onward. To compare the recovery time rather than volume, half of the initial spilled volume, 250 kl, is indicated in Fig. 11. It took 24 h and 14 h for the volume of the remaining oil in the SUR model and the SNR model respectively to reach half of the initially spilled volume. Thus, less time was required to remove half of the initial spilled volume in the SNR model. Figure 11. Comparison of the remaining oil volume of the SUR model and unskimmed zone and total zone in the SNR model.

Comparison of Calculation Results Using Actual Accident Information
The calculation was performed using the information of the actual accident; the inputs used for the calculations are listed in Table 2. In this scenario, 900 kl of Basrah light oil was assumed to be spilled [55]. The simulation time was 72 hours after the accident. Information of Normar 200TI skimmer, which was deployed on response vessel Hwangryong 208, was used in the simulation. All 10 skimming systems were applied from 1 h to 72 h. The size of the storage tank was assumed to be large enough to exclude time for emptying the storage tank.

Comparison of Calculation Results Using Actual Accident Information
The calculation was performed using the information of the actual accident; the inputs used for the calculations are listed in Table 2. In this scenario, 900 kl of Basrah light oil was assumed to be spilled [55]. The simulation time was 72 h after the accident. Information of Normar 200TI skimmer, which was deployed on response vessel Hwangryong 208, was used in the simulation. All 10 skimming systems were applied from 1 h to 72 h. The size of the storage tank was assumed to be large enough to exclude time for emptying the storage tank. The volume of remaining oil on the sea surface over 72 h is illustrated in Figure 12. The volume of remaining oil on the sea surface over 72 h is illustrated in Figure 12. The amount of oil remaining after 72 hours was 390 kl in the SNR model and 494 kl in the SUR model. Recovered oil volumes were 97 kl and 109 kl for the SUR model and the SNR model, respectively. It took 140 h and 28 h, for the volume of the remaining oil in the SUR model and the SNR model respectively to reach half of the initially spilled volume. The comparisons of the number of skimming systems and recovery amounts were carried out. To remove the remaining surface oil, the SUR and SNR models required a different number of skimming systems. The comparison results among two models and response field information are summarized in Table 3. About 34 skimming systems were mobilized in the actual accident [55]. There was a significant difference between the SUR model and field data. On the other hand, the SNR model showed a similar value of field data.  The comparisons of the number of skimming systems and recovery amounts were carried out. To remove the remaining surface oil, the SUR and SNR models required a different number of skimming systems. The comparison results among two models and response field information are summarized in Table 3. About 34 skimming systems were mobilized in the actual accident [55]. There was a significant difference between the SUR model and field data. On the other hand, the SNR model showed a similar value of field data.

Discussion
In terms of the response, various options can be employed at oil spill accidents. A skimming system can be a primary choice among mechanical recovery, chemical dispersion, and in situ burning. Numerous studies have been conducted to forecast spilled oil behaviors and establish response strategies for oil spill accidents at sea. However, current research remains insufficient on recovery potential estimation methods that consider oil weathering, slick behavior, and equipment efficiency. Many previous studies that aim to maximize the recovery focused on the window-of-opportunity [6,7,[31][32][33][34][35][36][37], the performance of skimmers [38][39][40], and clean-up potential [32,41,42]. Nevertheless, it is still required to study the response in terms of recovery capacity. The mechanical recovery potential estimation model of the present study can help decision-makers to set up the mobilization and deployment plan of the skimming systems. It is beneficial to estimate the potential of mechanical recovery systems in terms of recovery capacity in the response-planning phase. The present study is not designed to command the field response activities but to support the expert group by providing the response planning tool.
Two models were developed in this study. One is the SUR model that does not consider skimmed zone and local thickness variations, and the other is a model considering skimmed zone and thickness distribution (the SNR model). Three improvements that allow the model to reflect more specific conditions were implemented in the SNR model. The calculation of two models was conducted through the hypothetical and actual accident scenarios to compare the effects of the spatial and temporal improvements.
The calculation of the hypothetical accident case study revealed that the total recovered oil volume was 230 kl in the SUR model and 243 kl in the SNR model for 24 h. There was little difference between the SUR and SNR models in the oil recovery rate and volume from Figure 10. On the other hand, the remaining oil volume in the sea surface showed a considerable difference between the two models in Figure 11. Furthermore, it took 24 h and 14 h for the volume of remaining oil in the SUR and SNR models, respectively, to reach half of the initially spilled volume. Thus, less time was required to remove half of the initially spilled volume in the SNR model. It seems that the improvement of the SNR model well facilitates the sheen out of the remaining thin oil in the skimmed zones. Moreover, the actual accident case study showed a significant difference in calculation results. The SNR model calculated a similar number of skimming systems and recovery amounts with the field data.
Nevertheless, it is better to clarify the limitation of the SNR model. Skimming systems do not work like vacuum cleaners in actual spills. It means that the field operation of the skimming systems cannot precisely follow the mathematical model of the present study. The SUR and SNR models can mathematically calculate any size of the spill and recovery. However, that does not mean that both models can be applied to any size and shape of the slick. The spreading of oil is a very complicated process. The present study used the spreading models of Fay [51] and Galt and Overstreet [47] in calculating the slick area. Therefore, the uncertainty of the present study is the same as that of the previous models.
To identify the uncertainty of the SNR model, the effect of the parameter on the calculation results was analyzed. Mathematical modeling of the SNR model was developed to distinguish spatial differences in thickness between the skimmed and unskimmed zones. In addition, a spatial margin was added to prevent the collision of skimming systems, as shown in Equation (12). This marginal length might make some uncertainty on the number of skimming systems during the calculation. The uncertainty due to this marginal length was quantitatively analyzed for the 500 kl spill case study of Section 3.1. In the normal condition, the marginal length in the l1 direction, l1 margin , is assumed to be half of l1 s . The uncertainty of the calculation results was analyzed by varying the marginal length from the half to the double of the normal size condition. The number of skimming systems were compared in Figure 13. As the marginal length increased, the recovery time when a total of 20 skimming systems were operated decreased. Figure 14 shows the number of skimmed zones. The smaller the marginal length, the more the skimmed zones. The remaining oil volume of uncertainty study is summarized in Table 4.       Unfortunately, the skimming systems cannot entirely cover a wide geographic area in the case of very large-scale catastrophic spill accidents. It is essential to be guided with aerial observation and remote sensing. The deployment of the skimming systems could be directed to the thickest oil with aerial support to optimize recovery. The approach of the present study can be used to start with one slick and then tracks portions of it as skimming varies spatially and over time. That said, it is practically difficult to handle all the slicks in the same manner with the present model. Most numerical oil spill models are based on a Lagrangian approach. The Lagrangian spill model, aerial observation, and remote sensing could give the present recovery model the information on the location and thickness of slick. Then, the present recovery model can be used to estimate the mechanical recovery potential based on the collected information, but there are many technical challenges to handle this issue. It will be future work.
Finally, even though the SNR model cannot perfectly reflect all details of field conditions, it can allow a better understanding of the effect of nonuniformity in recovery activities on the response effectiveness and the amount of oil remaining floating. It will be helpful to make a response plan. Given enough accident information, this model can be utilized to estimate and reserve the necessary response resources in the event of oil spill accidents.

Conclusions
The recovery capacity of the mechanical skimming systems should be carefully considered in the response-planning phase. To make a suitable response strategy, it is necessary to assess mechanical recovery capacity. The SNR and SUR models of the present study consisted of weathering calculation and recovery calculation. The recovery calculation made a change in sea surface oil volume, and then the volume change made an effect in weathering calculation, and vice versa. To evaluate the spatial and temporal improvement of the SNR model, the recovery potential and its effects on the responses were compared and analyzed using case studies. The main findings are described below.
First, spatially uniform studies are limited in that the whole remaining oil in the sea surface is considered a single oil slick, with no spatial variations in thickness by mechanical clean-up activities. The first improvement in the SNR model considers these spatial changes in thickness by distinguishing areas where oil recovery was performed from those where it was not. Through this change, a more realistic mechanical recovery volume can be calculated using the thickness of each zone; the area of and weathering in each zone can also be calculated separately. This zone distinction in the SNR model allowed the model to represent thin slicks, which could not be easily distinguished in the SUR model, and reflect natural oil dissipation effects by clean-up activities. In the calculation result through the SNR model, a total of 367 skimmed zones were distinguished. In all, 347 zones of total skimmed zones and the unskimmed zone with thin thickness dissipated during 24 h, leaving 20 skimmed zones at 24 h in the case study.
Second, the SUR model takes into account neither the space occupancy by a skimming system nor spatial inference with other skimming systems. Therefore, any quantity of skimming systems could be deployed by the user during the calculations with no spatial constraints. That leads to an unaffordable response strategy in decision-making. The second improvement of the SNR model involved the calculation of both the space occupied by a skimming system and the interactions between skimming systems. This allowed the improved model to estimate the maximum applicable number of skimming systems over time and use this number in the recovery potential calculation. The unskimmed area in the SNR model decreased over time. As a result, the quantity of deployable skimming systems also decreased, and this decrease was reflected in the recovery volume calculation. The upper limit to the number of skimming systems depends on the slick area, and thus the number of skimming systems applied in the calculation changes with time. This enables the application of more realistic quantities of available skimming systems in the model. It will be very helpful in recovery strategy decision-making.
Third, spatially uniform studies calculate oil spreading using the initial spill volume, and changes in the residual oil volume due to weathering processes (such as evaporation and natural dissipation) and mechanical recovery were not considered. That can be a critical limitation in accessing the remaining oil mass balance, as the oil volume and area are essential factors in recovery potential calculations. The third improvement applied in the SNR model made it possible to consider the residual oil volume changes at the sea surface with time.