The application of coupled 3d hydrodynamic and oil transport model to oil spill incident in karawang offshore, indonesia

This study presents a coupled hydrodynamic and oil transport numerical model to study the spread of Karawang oil spills at sea due to well-kick failures. This model uses the 3D configuration of ROMS-CROCO in the Java Sea. The model has a resolution of 1 km, 25 vertical layers, and runs from January 2019 to September 2019. Temperature, salinity, sea surface height, ocean currents, and harmonic tides are derived from global models and applied to open boundaries. Hourly atmospheric datasets during model simulation are forced as flux input in the surface. The 3D profile of the flow generated by the model is converted to the GNOME oil transport model format as mover type input to disperse the oil. The hydrodynamic model shows that the result has a good agreement with in-situ data and observation with mean of correlation exceeding r>0.8 for sea surface height and sea surface temperature. Compared with radar satellites, oil spill dispersion shows the same scattered trend as satellite data. Backward modelling shows oil particles returning to the initial spill location. The oil spill was moving westward, and some are stranded on the coast between Karawang and Bekasi.


Introduction
After the marine environmental pollution incident due to leaking the underwater pipeline from the Lawe-Lawe terminal to the PT Pertamina refinery facility in Balikpapan Bay, East Kalimantan in 2018, the oil spill incident happened again in Karawang, West Java, on July 12, 2019. The oil and gas bubbles leak from the offshore drilling location of Rig YYA1, owned by PT Pertamina Hulu Energi Offshore North West Java (PHE ONWJ), Karawang. As a result, the oil spreads and extends westward to Bekasi, West Java, even to Jakarta's Seribu Islands. According to Vice President Relations of PT Pertamina Hulu Energi (PHE), as the mother company of PHE ONWJ, during a press conference in Jakarta, Wednesday (17/7/2019), the gas bubbles appeared due to a well kick occurring in the YYA-1 reactivation well [1,2]. In general, a well kick is an event where formation fluid enters the well during drilling activities due to the formation pressure being higher than the hydrostatic pressure of the mud.
The incident potentially harmed marine ecosystems and coastal communities living on the coast of Karawang. Crude oil pollutes several beaches, and marine tourism sites are closed because seawater is polluted with oil, and tourists cannot enjoy swimming on the beach. The community might lose their livelihood as fishermen because the fish catchment location was polluted by oil and gas deposits. Meanwhile, salt cannot be produced. Polluted seawater kills fish and shrimp at the ponds. The fish traders also suffer losses because no one wants to buy stingy fish. The Local Agency of Environment and Sanitary in Karawang Regency noted that 232,000 mangrove trees were exposed to oil. In addition to the loss of catch, oil and gas deposits containing hazardous and toxic materials (B3) also trigger adverse effects on coastal residents' healths. MMAF stated that as many as 7782 fishermen spread across 12 villages in Karawang Regency and 2200 fishers and fish farmers in Bekasi Regency, especially in Pantai Bakti and Pantai Bahagia villages, were affected by this oil pollution [9][10][11].
Various studies on oil spill pollution in Karawang have been carried out, i.e., detection using C-band SAR Sentinel-1 [4], detection using C-band SAR Sentinel-1 with the adaptive threshold [5], detection using multispectral image Sentinel-2 [6], detection using multispectral image Aqua-MODIS, Landsat-8 and Sentinel 2B [7], oil spill trajectory by combining SAR Images and ocean currents and wind [8], oil spill trajectory using numerical model 2D [9], the impact of oil spill on fishes [10], social impacts towards the people living in coastal areas [11], absolute responsibility for environmental loss of salt farmer, fishermen's economic loss, compensation on marine tourism business [14]. However, we do not see any studies on oil spill trajectory using a numerical model of coupled hydrodynamic and oil transport by using realistic oceanic and atmospheric forcing. In this model, oil particles are transported by ocean currents. In practice, however, most oil spill models only consider barotropic ocean currents and ignore the effect of baroclinicity (3D). This study aims to use a realistic coupled three dimensional hydrodynamic ROMS-CROCO model and GNOME oil spill model to assess the transport of oil spill sources in Karawang waters based on the 2019 event. The model results will be compared to oil spill studies conducted using SAR satellites.

CROCO Model
This study uses the CROCO version (www.croco-ocean.org) of the ROMS-UCLA model [15] which is based on ROMS AGRIF [16] and developed by the French Institute of Research for the Development (IRD). CROCO is a free surface hydrostatic C-grid model with a terrain following-coordinate system and an efficient split-explicit scheme for separating barotropic and baroclinic terms. CROCO also provides MATLAB-based pre-processing tools [17] that can be used to create a model grid, interpolate atmospheric and oceanic data as boundary and forcing input, and set tides and rivers input in the grid model. CROCO supports MPI and OPENMP computations.
We use the default CROCO regional configuration by enabling tides without river forcing. For vertical advection, the momentum term employs the 3rd-order upstream biased advection scheme and the 4th-order centered advection scheme with harmonic averaging. The K-profile (KPP) parameterization is used for vertical mixing, and a split and rotated 3rd-order upstream biased lateral tracer advection scheme is used for lateral tracer advection. Flather open boundary conditions for barotropic velocities and radiative open boundary conditions for both baroclinic velocities and tracers are imposed on the lateral boundary. Vertical coordinate use NEW_S_COORD [18], which corresponds to Vstretching=4 in ROMS. Inter-annual climatology scripts are also available, which can help speed up the operational process. The CROCO model has been applied to the Indonesian sea to study the circulation, eddies, and internal tides [19][20][21][22][23].

GNOME Model
GNOME is an open-source oil spill-trajectory model developed by the Hazardous Materials Response Division (HAZMAT), Office of Response and Restoration1, National Oceanic and Atmospheric Administration (NOAA) [24,25]. In this study, the GNOME desktop version 1.3.11 integrated with the GIS module was used to pre-process and display the results. In order to manipulate the risk of oil spills, the Python version (PyGNOME) can be used with a more comprehensive petroleum database, weathering module and can be used for automation. The basic GNOME algorithm is based on the lagrangian trajectory. If the oil spill is detected from satellite data or field evidence, it can be modeled through a reverse model to find the original source. After the original location is known, the forward trajectory predicts the oil spill [26,27].
The Lagrangian algorithm provides both best estimates assuming that all model inputs are true and minimal regret solutions for uncertain predictability. The best guess solution is generated with the assumption that all the inputs to the model are correct. The uncertainty in GNOME is a statistical set of trajectories that sample any prediction errors related to the trajectory prediction by using a random walk number. The required input is ocean current and wind data. The ocean current input is a "mover" variable that can be extracted from the output of hydrodynamic models such as ROMS, HYCOM, and FVCOM. Wind mover is entered as a constant or time variable in date, speed, and direction.

Model input and Scenario 2.3.1. Model domain and bathymetry
Domain spans along the Tangerang coast to the Karawang coast in West Java, with geographic coordinates of 106.30E -106.80E and 5.30S -6.40S. The horizontal grid resolution is 1 km and divided into 25 vertical layers. The water depth in the grid model is obtained from the spline interpolation of the GEBCO 2020 data (www.gebco.net) with a resolution of 15 arc seconds. The wetting and drying algorithm is not activated in the model. Thus, to avoid "dry" grids during ebb tide, the coordinate formation's minimum depth (hmin) is set to 3 meters. Maximum depth is found 75 meters on the north coast of Tangerang ( Figure 1). The vertical grid is set to be a uniform layer from the surface to 2 meters (hc = 2). Grid generation is created and smoothed using CROCO-TOOLS with a denser vertical mesh parameter on the surface (theta_s = 6, theta_b = 1).

Model initial condition and forcing input
The open and lateral limits of baroclinic ocean currents (U), (V), deep integrated ocean currents (ubar, vbar), salinity, sea level (zeta), and temperature are extracted from the NEMO model of daily data with resolution 0.083 degree, which is part of the Copernicus Marine Service (CMEMS). The same data was used to create the initial conditions for December 1, 2018. CMEMS has a partial stepwise vertical coordinate system that requires interpolation and smoothing to S coordinate vertical grid in CROCO. For surface forcing, hourly ERA5 data is used for heat flux, surface stress, precipitation, evaporation, height, and direction of significant waves, and these data are interpolated into the first layer of the model. Time filtering is necessary to smooth out the time variation of ERA5 to avoid blow-up when running the model. To the best of our knowledge, our study is the first to publish the model with hourly ERA5 and CMEMS baroclinic forcing in the CROCO model.
Tidal forcing is obtained from TPXO9-atlas [28], including 10 main components such as M2, S2, N2, K2, K1, O1, P1, Q1, Mf, and Mm. CROCO requires the amplitude and phase of each tidal component of the sea surface elevation and the axis of the primary and secondary tidal flow of the barotropic tidal flow, respectively. River forcing is not used, and we are still developing to include daily freshwater flux from satellite data.

Model scenario and coupling
In order to obtain the real flow conditions at the time of the oil spill, we conducted a 13-month CROCO simulation (from December 1st 2018, to December 31st 2019). We take the first month as a spin-up period to stabilize the water mass and tides in the model. The effective OPENMP computation in CROCO allows performing calculations on a 1km grid without performing Agrif multi-grid calculations.
The simulation was run using an Intel i7 PC with 16GB of memory. Computation time is 30 minutes/month simulation with baroclinic and barotropic timestep 600 seconds and 10 seconds, respectively. The oil spill model scenario is based on the oil spill incident at (6.09417° S;107.6257 °E) on the north shore of Karawang regency. The initial oil spill was detected around the oil platform (YYA1) managed by Pertamina Hulu Energy (PHE).
On the GNOME desktop, there are multiple oil-type options. Several types of hydrocarbons can be used, such as gasoline, diesel, unweathered, and medium crude oil. Our study is limited to examining the distribution of oil particles due to the advection of ocean currents. We use the non-weathered type of oil to neglect the physico-chemical reaction on the fate of the oil particles. However, we did simulation tests for all types of hydrocarbons (results not shown). Furthermore, GNOME uses a regular grid, where the U and V points are located at the same point. At the same time, ROMS uses a staggered C-grid by placing U and V velocities at different points. For GNOME to read the CROCO results, the CROCO velocity points must be interpolated into the GNOME coordinate system. For model coupling, we use python tools developed by The National Oceanic and Atmospheric Administration (NOAA), Office of Response and Recovery (ORR), and Emergency Response Department to convert CROCO output to GNOME grid.
We take advantage on satellite processing image on 18 th July 2019 from [10] as a baseline oil spill pattern to reconstruct our model scenarios. The exact time of the first oil spill incident in YYA1 is not shown in [10], however we use the exact time of the first incident time on July 12 th 01.30 am from [5]. There are forward time and backward time modelling scenarios in this study. The forward modelling runs for 10 days on July 12th 2019 at 01 am to July 22nd 2019 at 23 pm. The backward modelling starting on July 18th 2019 at 01 am with time step back until the first oil spill incident time on July 12th 2019 at 01.30 am. Both scenarios are carried out by inserting 0.5 barrels of oil as initial input in GNOME model. Initial oil spill position for GNOME forward modelling is in YYA1-oil platform and backward modelling use initial oil spill in the area furthest from the initial spill site on 18th July 2010.

Model validation
Model validation is performed by comparing daily average data from model output and satellite observation from January 2019-December 2019. Degree of correspondence and association between the model results (m) and the observational or satellite data (s) are calculated using RMSE and correlation equation (r) as follows: n and overbar denotes number of data and total sum respectively. Taylor diagrams is used to graphically summarize the RMSE, standard deviation and correlation coefficient. The selected location is the average area in Karawang offshore. For the characteristics of sea surface temperature, we compared the model results with 1 km ultra-high-resolution multi-scale data (MURSST) [29] The latest CMEMS data from the global hourly model are used to compare current speeds, which comprise the total data for zonal velocity (U) and meridional velocity (V) by tides, baroclinic pressure, and euler stokes drift. Figure 2 shows that CROCO can well represent and show the same trend of sea surface temperature. Good correlation are achieved between model and satellite data with mean correlation r=0.9, r=0.95, r=0.91,r=0.978 for MUR-SST, GHRSST, OSTIA and CMEMS respectively. RMSE is less than 0.3 for all comparison (Figure 6b). Timeseries plot of CROCO sea level data variable (zeta) and absolute dynamic topography from altimetry satellites (Figure 3) are in good agreement with correlation r=0.805 and RMSE 0.03 ( Figure  6a). Interestingly, CROCO's zonal (U) (Figure 4) and meridional (V) ( Figure 5) velocities are close to the velocity from the global CMEMS model with RMSE=0.2 for (U) velocity (Figure 6c) and 0.08 for (V) velocity (Figure 6d). However, we found a low correlation coefficient r=0.5. Unfortunately, we do not have in-situ measurement data in the model area.    Figure 7 shows that tidal height from the model has a good agreement with BIG tide data. The calculation of the square RMSE value=0.1 and obtained a 99% confidence level.

Current condition during oil spill incident
We use surface current from CROCO model result as mover variable in GNOME. Figure 8 shows the snapshot of sea surface temperature and the daily average current vector. The current pattern when oil spills occurred on July 12, 2019 was moved mainly to the west between 0.1-0.5 m/s. Along the coast, the current is weak due to the influence of bottom friction. The current pattern is an important factor that is responsible for transporting oil spills.  Figure 9 shows the backward trajectory of the oil spill after  8 1 barrel of oil was released on July 18, 2019. The black dots represent the best estimate of the oil particles distribution and the red dots represent the 10% uncertainty of the best estimate of particle distribution.
The trajectory of the oil particles was plotted daily, and it was found that the oil spill position is at the point where the initial position of the oil spill is reported. The initial time at this point is July 12, 2019. GNOME provides sufficient evidence to identify the original location of oil spills. GNOME provided enough evidence to determine the initial location of the oil spill. As seen in satellite images, the oil particle distribution cannot converge to one point. However, more field data and information from the public or fishermen can be used to support the model results. Figure 9. Comparison of satellite analysis and GNOME backward trajectory result. Radar satellite data from [8] (a), oil spill initial condition for backward model at 18 July 2019 (b), Oil spill particle movement every 1-day from 17 July(c) ,16July(d), 15 July (e), 14 July (f), 13 July (g). Position of oil spill trajectory at 12 july 2019 (h). Initial location of oil spill incident marked by yellow star.

Forward trajectory modelling
We first examine the trajectory of oil particles by comparing the model result and satellite image ( Figure  10). Visual comparisons indicate that GNOME has a scattering trend close to the satellite result. However, there is a difference in the details of the orbital pattern. It can arise from the influence of the non-linear process, as well as the impact of oil spill prevention activities carried out by the owner of the oil platform. The results of the simulation indicate that the wind can spread the wider particles. The red dots start to disperse from black dots after 5 days from the beginning of the simulation and eventually stranded on the beach. We assume that the red dots can be imagined as a dispersion process from the value of highly concentrated oil particles to low concentration oil particles The beach on the coast of Karawang (Buntu River) has direct impact soon after the incident. After 5 days, the oil particles continued moving towards the west and on the coast of Bekasi Regency. After 10 days of the case, the oil particles entered the Muara Gembong of Bekasi Regency. This is 60 km from the initial position of the oil spill. Unfortunately, we do not have any scientific data on oil particles stranded along the coast from Karawang to bekasi. Figure 10. Oil trajectory from GNOME result(a,b,d) and radar satellite data from [8] (c).

Conclusion
The results of the CROCO simulation show good validation through observation. The GNOME oil spill and the current CROCO model input are used to very well simulate the oil spill off the coast of Karawang. The results of the hydrodynamic model provide an accurate comparison with satellite data. Additionally, the GNOME trajectory results showed significant results from satellite data. The results of this study can be used as inputs for emergency response units to address risks in oil and gas activities.