Integrated 3D numerical modeling during La Niña and El Niño events using Regional Ocean Modeling System (ROMS) in Makassar Strait: preliminary study

This research showed the difference of thermocline layer depth during La Niña and El Niño events in the Makassar Strait, the main passage connecting the Pacific and Indian Oceans. La Niña event represented by 2011 and El Niño event represented by 2015. Three-dimensional numerical modeling using Regional Ocean Modeling System (ROMS) is used to determine ocean dynamic characteristics. The tides and atmospheric factor used as generating forces for coupled Ocean-Atmosphere-Wave-Sediment simulation model. Tidal elevation verification results from January to February 2017 showed a high confidence with RMSE = 0.148 m and MAPE = 13.51 %. Ocean circulation showed that water flows from north to south of the Makassar Strait in both conditions, El Niño and La Niña, indicating a transport from the Pacific to the Indian Oceans during those time. Two-dimensional seasonal sea surface temperature (SST) was used to determine the variability of thermocline layer depth. In general, based on simulation results, the depth of thermocline during La Niña was deeper than that during El Niño. The upper threshold of thermocline layer during El Niño phenomenon was shallower (average 27.29 m) compare to the threshold during La Niña (average 39.25 m). So do the lower threshold, during El Niño was shallower (average 160.47 m) compare to the threshold during La Niña (average 165.70 m). The thermocline thickness during El Niño (133.19 m) was found to be thicker compared to La Niña (126.43 m).


Introduction
Makassar Strait, geographically, is water that connects the Pacific Ocean and the Indian Ocean. Many unique phenomena, such as Indonesian Troughflow (ITF) [1,2], tidal mixing caused by steep slope bathymetry [3], internal Kelvin wave propagated from Lombok Strait [4], and effects of ENSO events to water characteristics [5,6]. Makassar Strait is also affected with seasonal trade wind that helped mixing in surface layer between sea water and fresh water from river [7].
To understand its complexity, a wind-wave-tide integrated numerical model needed to describe physical phenomenon in Makassar Strait. One of high level model that many others used in previous study is the Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) [8,9,10]. Figure 1 showed a schematic high level model COAWST which is comprised of the Model Coupling Toolkit to exchange data Figure 1. Basic schematic of the COAWST modeling system, include their 'link' between each other, ROMS as ocean model, SWAN as wave model, and WRF as atmospheric model [10].
Objective of this study is to understand the variability of ENSO phenomena, for both El Niño and La Niña conditions, to understand the characteristic water masses changes in Makassar Strait. To know its difference between two conditions, two variables were used as determinant factors; they are thermocline depth and sea surface temperature (SST). We compared these variables among three conditions, the La Niña in year of 2015, normal condition in year of 2013, and El Niño in year of 2011. The interest of this research area showed using figure 3. We compared the model elevation result with tide observation as verification

Atmospheric Forcing
Atmospheric forcing such as surface wind, air pressure, air temperature, air humidity, net fresh water flux, rainfall rate, net long-wave radiation flux, and solar shortwave radiation flux extracted from European Centre for Medium-Range Weather Forecasts (ECMWF) (http://apps.ecmwf.int/datasets/) and applied to model for each scenario. Surface wind is velocity of wind 10 meters above sea surface elevation and used as generating force for surface current circulation. Surface wind was imposed every 3 hours for a-year long simulation. Atmospheric forcing were computed by ROMS using bulk-fluxes formulation internally and turbulent fluxes for wind, heat, and moisture are computed using Monin-Obukhov similarity theory [13]. These atmospheric forcing were imposed to the model every 3 hours time step. Figure 3 showed bathymetry that was used for the model area. It was extracted from global topography data fusion of NASA Shuttle Radar Topography Mission (SRTM) [14] land topography with measured and estimated seafloor topography (SRTM15_PLUS) (ftp://topex.ucsd.edu/pub/srtm15_plus/). This data is corrected by sounding [15] and gravity data [16] and modified from SRTM30 product distributed by USGS EROS data center. The grid resolution is 30 second which is roughly one kilometer. Land data are based on the 1-km averages of topography derived from the USGS SRTM30 gridded DEM data product created with data from the NASA Shuttle Radar Topography Mission. GTOPO30 data are used for high latitudes where SRTM data are not available. Version 1 of SRTM15 is based on V11 of the SRTM30_PLUS grid. Basically the land data is completely new based on a combination of SRTM, ASTER, and CryoSat-2 ice sheet topography. The ocean data are the same as SRTM30_PLUS but required more extensive editing to remove the bad points mostly along edges of the swath data. Ocean data are based on the global seafloor topography from satellite altimetry and ship depth soundings [17]  . It is produced on a 2° × 2° grid with spatial completeness enhanced using statistical methods. This monthly analysis begins in January 1854 continuing to the present and includes anomalies computed with respect to a 1971-2000 monthly climatology. The newest version of ERSST, version 5, is based on optimally tuned parameters using the latest datasets and improved analysis methods [19,20,21]. We compared model results from year 2017 model using field observation data and used the setting to be applied to the three scenarios of model. We compared 3-month average temperature profiles (SST and vertical temperature profiles) between each scenario. We collected 7 random points in the model from north part to south part to see the vertical profiles and compared them between La Niña, El Niño, and normal scenario. Thus thermocline layer calculate using gradient criterion 0.05 o C/m to determine the depth of upper and lower threshold of thermocline layer [22].

Bathymetry
To see the correlation between model and observation data, we used Root Mean Square Error (RMSE) and Mean Percentage Absolute Error (MAPE) formula that shown by equation (1) and (2) [23] below 2 , , with t is time step model in hour, n is total data, A is actual value of observation data, and B is model result.

Hydrodynamic Model ROMS
The simulation of hydrodynamic process in domain area were conducted by using the numerical model Regional Ocean Modeling System (ROMS). This model was quite popular among modeler, scientist and researcher who wanted to study the coastal application and developed from Princeton Ocean Model (POM).
Many literature described this model capability especially in the regional ocean domain. ROMS is a threedimensional, free surface, terrain-following numerical model that solves finite difference approximation of the Reynolds-averaged Navier-Stokes (RANS) equation using the hydrostatic and Boussinesq assumption with a split-explicit time stepping algorithm [24,25,26]. It uses a horizontal curvilinear Arakawa "C" grid and vertical stretched terrain-following coordinates [26]. This model also can be configured depending of users application which has several choices for advection schemes, pressure gradient algorithms, turbulent closure, and many types of boundary conditions.
with the continuity equation: and scalar transport: These equations are closed by parameterizing the Reynolds stresses and turbulent tracer fluxes as where KM is the eddy viscosity for momentum and KH is the eddy diffusivity. Eddy viscosities and eddy diffusivities are calculated using one of five options for turbulence-closure models in ROMS: (i) Brunt-Väisälä frequency mixing in which mixing is based on the stability frequency; (ii) a user-provided analytical expression such as a constant or parabolic shape; (iii) the K-profile parameterization [27], expanded to include both surface and bottom-boundary layers [28]; (iv) Mellor-Yamada level 2.5 (MY2.5) method [29]; and (v) the generic length-scale (GLS) method [30] as implemented in [31] that also includes the option for surface fluxes of turbulent kinetic energy due to wave breaking. and for this study, we applied option (v) to

Model Design
This model have 160 east-west grid cells and 208 north-south grid cells with 2.6 km grid spacing uniformly. This model have 413.4 km wide and 538.2 km long with width area is about 222.49 10 3 km 2 . We used masking to determine which one is land and which one is ocean using 1 and 0, respectively. This model used sigma vertical coordinate with 30 layers and the shallowest depth is -5 m and deepest is -2734 m. Vertical depth references is mean sea level and no wet and dry scenarios used. We used 3 momentum open boundary (north, south, and east, see figure 4) with tidal elevation as major forcing input and coupled with atmospherics data from ECMWF as mention in Section 2. No nesting model used and no river discharge in domain model used. Time step was calculated using Courant-Friedrichs-Lewy condition (CFL) [32] and obtained 123.32 second. We plotted ocean current pattern 3-months averages based on ONI for 3 scenarios (El Niño, La Niña, and Normal) and so do for SST. Also we plotted temperatures vertical profiles and compared them to see upper threshold and lower threshold for thermocline layer.   Figure 4 showed a tidal verification between observation data (red line) and ROMS model results (blue line). It showed that model still under-estimated with observation data have higher value than output model. This verification is still good with pattern of tidal is match and RMSE value is about 0.148 m and MAPE is around 12.22 %.

Surface Current Circulation in Makassar Strait
Ocean pattern circulation in Makassar Strait was influenced mostly by wind, tidal forcing, and gradient pressure difference between West Pacific Ocean and East Indian Ocean. Figure 5   Model results showed that current always move along from north part to the south throughout the year. In July, August, and September with ONI is -0.7 at 2011 (weak La Niña condition) (figure 5 (c)), average current velocity is about 0.72 ms -1 . In April, May, June, where ONI value is -0.3 at normal condition, figure  6 (b), average current velocity is about 0.88 ms -1 . And when weak El Niño (January, February, March) with ONI value is 0.6, model showed ocean current velocity is about 0.64 ms -1 . It showed that velocity when La Niña is slower than El Niño condition but still less powerful than normal condition. When strong La Niña happens with ONI -1.1, depth-average value of velocity is about 1.36 ms -1 . And when normal condition at October, November, December with ONI -0.2, average velocity is about 0.89 ms -1 . And maximum velocity showed at strong El Niño condition with ONI -2.5 at ONJ with value 1.56 ms -1 (darker green in figure 7 (d)). Current velocity showed also that higher velocity is placed at deeper bathymetry thus far from north Water masses passing trough from west Pacific Ocean to Makassar Strait using narrow gap at Talok Cape and move to Sangkulirang Bay, East Borneo ( figure 9). When west monsoon occurs, there are counter current that move from Java Sea to west part of Makassar Strait. This fresher water masses mixed with saltier water masses from Pacific Ocean therefore it became fresher and disappeared at southern part of Makassar Strait.   This result is not quite correct because SST supposes to be warmer during La Niña than normal condition with warm pool should move to west Pacific and effecting Makassar Strait water mass. This is caused by model input is not good enough to reproduce a warmer SST during La Niña condition. The atmospheric forcing from peak of La Niña period which is OND on previous year is neglected. Another explanation possible is during La Niña SST supposes to be warmer but peak location of La Niña is do not match in the model area in Makassar Strait possibly will get a delay in SST signal. So, SST during La Niña is colder than normal period which is model result do not explain the realistic condition. Reversely, during El Niño condition, SST is colder than La Niña and normal condition with 25.7 o C -31.23 o C. There are some area that always colder than any other area where placed at south-east part of the model and believed that was caused by input of atmospheric parameter such as heat flux. A warmer area was found at west part of the model where is water mass from Java Sea passing trough along coastal area to Makassar Strait. Higher SST found when strong La Niña 2011 occurred with ONI = -1, January February March (JFM) period. Colder SST found when weak El Niño 2015 occurred with ONI = 2.2, October, November, December (OND) period.

SST vs ENSO
To see distribution and variability of thermocline depth, we took 7 sample points and plot them in vertical profiles of temperature ( figure 13). Seasonal vertical profiles of temperature showed by figure 13 (a) for La Niña, figure 13 (b) for normal condition, and figure 13 (c) for El Niño. We calculated upper threshold (UT) and lower threshold (LT) of thermocline layer using gradient criterion ≥ 0.05 o C/m [21]. This criterion is good to used at South China Sea and around seas with 200 m depth or more [21]. We also calculated thermocline thickness during 3 different conditions with simply minus upper threshold to lower threshold. Thermocline layer at Makassar Strait found between 1.  Table 3 showed variability of thermocline thickness, upper and lower threshold during 3 different conditions. Model result showed that thermocline layer during La Niña condition in Makassar Strait is found deeper than during El Niño condition. During La Niña, the thermocline is stable and very stratified and sea level height is lower than normal over the eastern Pacific, resulting in an increase slope of the ocean surface across the basin. Water masses from west Pacific Ocean pushed down the water masses at Makassar Strait hence the thermocline layer affected by it and became deeper as in line with strong La Niña occurred.