Modelling Tidal Flow Hydrodynamics of Sunda Strait, Indonesia

In the past years, Indonesian people put more attention to Sunda Strait located between Java and Sumatra Islands, one of the busiest straits occupied with residential, recreational, fisheries, transportation, industrial and mining activities. Previous works on numerical modelling of tidal flow hydrodynamics of the Sunda Strait have resulted in good agreement against field data; however, the calibration of the models used was not described in detail. This paper presents the process of setting up the model, extensive calibration, validation and prediction of tidal currents for the Sunda Strait. A two-dimensional tidal-driven model is constructed using Delft3D, an opensource developed by Deltares. Four different bathymetry datasets, four different boundary condition configurations, and various bed roughness values are used, and their suitability in predicting tidal water level and current are investigated. It is found that changing the bathymetry and boundary conditions improve the model validation significantly. GEBCO_2019 bathymetry dataset outperforms the Batnas, even though it has a coarser resolution. For boundary conditions, the combination of water level and current velocity results in a better validation compares to using water level or current velocity only. However, the bed roughness shows an insignificant influence in predicting tidal conditions. The averaged current velocity is lower at the Southern than the Northern side of the strait due to a larger cross-section, consequence of deeper water. High tidal currents of magnitude around 2 m.s-1 are seen at the bottleneck of the strait.


Introduction
Numerical modelling has become a common tool to understand tidal hydrodynamics in a region. Compared to field measurements, numerical modelling is relatively cheaper and capable in presenting the spatial and temporal variations of parameters. Once proven to be reliable, the model can be easily modified to recreate the real-time, forecast, or preceding data. One challenge that need addressing is the model's calibration which aims to tune and prove model's reliability.
This paper aims to perform a sensitivity analysis to enhance model reliability. Sunda Strait, located between Sumatra and Java Island of Indonesia, is chosen as the study location. This strait is included in the national strategic region since it is highly occupied by human activities such as for residential, recreational, fisheries, transportation, industrial and mining (Authority of Lampung Province, 2010; Authority of Banten Province, 2017).
Previously, some numerical modelling studies have been undertaken in the Sunda Strait using MIKE21, Princeton Ocean Model (POM), and Delft3D. Both MIKE21 models by Welly et al. (2012) and Novico et al. (2015) investigated the potential of tidal energy and presented an acceptable validation against field measurement. Orhan et al. (2019) also studied the energy potential using Delft3D incorporating the density induced flow. The model validation against water level showed a good result, and the comparison between the baroclinic and barotropic surface flow presented significant differences. Rahmawati (2017) used POM and modelled the ocean circulation and land-based pollution to determine the fisheries loss around Sunda Strait. The water level validation presented showed a good agreement with around 10 cm discrepancy.
The above models presented a decent validation against field data; however, the model calibration was not highlighted. Therefore, in this study, the sensitivity of model inputs, i.e. bathymetry, boundary conditions, and bed roughness is investigated using Delft3D. The tidaldriven model is run depth-averaged and with several bathymetry datasets, boundary conditions configurations, and bed roughness alternatives. The pattern of the tidal flow is explored. Further, the findings are expected to aid other researchers in enhancing the model performance in general and understanding the hydrodynamics of Sunda Strait.

Materials and Methods
In this paper, numerical hydrodynamic modelling was carried out using Delft3D-Flow of Deltares (Deltares, 2020). The method is described into three sections. Model description presents the overview of the study location and model preparation. Model schematization explains the variation applied to investigate the sensitivity of bathymetry, boundary conditions, and bed roughness inputs. Model validation describes the validation procedure and the available field data.

Model description
The area of study was at the Sunda Strait, Indonesia which is separating Banten Regency in Java Island and Lampung Regency in Sumatra Island with the strait stretch for around 28 km. See Figures 1a and 1b. The strait links the deep Indian Ocean and the shallow Java Sea which makes the region interesting.
The numerical modelling is performed using Delft3D-Flow which has been widely used for various study cases in Indonesian waters (Kurniawan et al., 2017;Orhan et al., 2019;Takagi et al., 2019). Delft3D-Flow is an open-source developed by Deltares for hydrodynamic modelling, which solves the Navier-Stokes equations for an incompressible fluid, within shallow water and the Boussinesq assumptions (Deltares, 2020).
The model uses the geographic or spherical coordinate for the horizontal coordinate system. In the vertical space, this model only simulates the depth-averaged velocity. The modelling domains are shown in Figure 1b. Three domains were used which have different resolutions, denoted as coarse (900 x 900 m 2 ), medium (300 x 300 m 2 ), and fine domains (100 x 100 m 2 ). Each was represented by the area inside the green, blue, and red dashed line in Figure  1b. The fine domain is intended to represent the small islands at the bottleneck of the strait, see Delft Dashboard (DDB), an integrative tool for quick setup a coastal hydrodynamic model is used to construct the initial model. With domain coverage in Figure 1b, the initial model applies GEBCO_08 for the bathymetry, TPXO 7.2 for water level boundary conditions, and default manning 0.020 for bed roughness (Nederhoff et al., 2016). Also, a measured bathymetry dataset around Panjurit and Sangiang Island are incorporated in the model. The model applies cold starts and other constant values, such as 9.81 m.s -2 for gravity, 1025 kg.m -3 for water density, and default 1 m 2 .s -1 for eddy viscosity. The Courant Friedrich Lewy number is designed lower than 10 to achieve stability (Bijvelds, 2001). Thus 15 seconds is used for the computational time step.

Model schematization
The modelling aims to investigate the sensitivity of three parameters, which are the bathymetry, boundary conditions, and bed roughness. The testing of these three parameters was performed sequentially, left to right, as shown in Figure 2. As mentioned, the initial model applies the default model setting in DDB, which is GEBCO_08 for bathymetry, water level boundary conditions, and 0.020 manning coefficients.
There were four bathymetry datasets tested. Besides GEBCO_08, the GEBCO_2019, SRTM15+, and Batimetri Nasional (Batnas) are utilized. The GEBCO data were provided by the International Hydrographic Organization (IHO) and the Intergovernmental Oceanographic Commission (IOC) (GEBCO Compilation Group, 2019). SRTM15+ was prepared by Tozer et al. (2019). Batnas was made available by the Indonesian Geospatial Information Agency or BIG (Badan Informasi Geospasial, 2018). Each dataset was applied to the model and run individually. Bathymetry with the best validation will be used for the next stage.
As shown in Figure 1b, there were four sides of boundaries, the North (N) and East (E) are in the Java Sea and the South (S) and West (W) close to the Indian Ocean. At the initial model, all boundaries are prescribed in water level. Further, in the sensitivity analysis, three more configurations were introduced, as shown in Figure 2, which shows the combinations of water level (WL) and current velocity (CV). The harmonic constituents of water level and current velocity were generated using TPXO 7.2. Again, the best configuration was used for the next stage.  In Delft3D-FLOW, user can specify whether manning, chezy, or white-colebrook formulation to be used (Deltares, 2020), and in this study, manning formulation was employed. The manning 'n' coefficient was used, which is in m -1/3 .s unit. There were five values tested from 0.016 to 0.032, each with 0.004 differences. Later, the best validation among all will be denoted as the final model.

Model validation
Validation was conducted to quantify the reliability of the models, using root mean square error (RMSE) and correlation coefficient (r). RMSE value close to 0 indicates a small discrepancy between model and field data. For the r, the desired value is close to 1, meaning the model and field data have a positive linear correlation. Meanwhile, 0 and -1 indicate a no linear correlation and a negative linear correlation, respectively.
There were four water level observation data used to be compared with the model output, provided by BIG. The data are available hourly from June 1 to 15, 2019. The stations were Ciwandan, Marina Jambu, Sebesi, and Panjang shown in Figures 1b and 1c as black circle denoted as WL1 to WL4, respectively.
Two current velocity data were also used to validate the model. The primary velocity data denoted as red circle CV1 in Figure 1c was obtained from acoustic doppler current profiler (ADCP) measurement at Anyer Water, which also used in Welly et al. (2012). The depth-averaged data was hourly and available from October 9 to 14, 2010. Since the data is not ideal with a short period (Firdaus et al., 2019), an additional velocity data from Indonesian Navy Tide Tables, called Dishidros or DISH is also incorporated. The tide table data is a result of prediction and available hourly throughout the year 2019 (Dishidros, 2019). The station was located at Ketapang Water shown as red circle CV2 in Figure 1c. To accommodate the field data, the models were run in two-time frames, October 7 to 15, 2010 and May 30 to June 15, 2019. The model was run two days earlier to allow the model spin-up.

Result and Discussion
The validation of the initial model is presented in Figure 3. As seen in Figures 3a to 3d, the water level data shows a good agreement. While in Figures 3c and 3d, the current velocity is not showing a satisfying correlation. Model tuning by applying different bathymetry, boundary conditions, and bed roughness is expected to improve the validation result.

Influence of bathymetry
The calibration summary for each model is presented in Table 1. Previously researches experienced the inaccuracy of GEBCO in the Java Sea (Koropitan and Ikeda, 2008;Firdaus et al., 2019), however, in this model, it is seen that the GEBCO_2019 improves the current velocity agreement significantly. Despite a finer resolution, the Batnas fails to show a better result at CV1 station and it is expected due to the seabed elevation inaccuracy around CV1. Meanwhile, for the water level errors, the numbers do not notably change compared to the initial model.

Influence of boundary conditions
As presented in Table 2, Configuration C (Conf_C) produces the best result where water level boundary conditions are applied at the South and North side while current velocity are prescribed for the West and East side. Despite a higher RMSE at CV2 station, the Conf_C is more preferred since it improves and gives the best RMSE and r values at CV1 station, which is primary current velocity data.

Influence of bed roughness
With Conf_C as boundary conditions, five manning roughness values are tested. As seen in Table 3, the optimum bed roughness value obtained is 0.024 for this case. In comparison, Novico et al. (2015) obtained that the optimum roughness is 0.031 used in MIKE21 model. While this study is still limited in the uniform roughness, other options are available by using varying bed roughness based on water depth (Dias and Lopes, 2006) and bedforms (Pramono, 2005) which is potential to be used in the next research. Table 4 summarizes the comparison between the validation of the initial and final models. In water level' RMSE and r, comparable numbers are seen. However, the current velocity numbers are significantly improving from the initial to the final model. Figure 4 presents the current velocity validation of the final model. Compared to Figures  3e and 3f, it is seen that the calibration improves the model quality significantly. However, the result at CV2 in Figure 4b indicates that the model underestimates the amplitudes in field data, though the trend is nicely represented.
In addition, Figures 5a to 5c are dedicated to displaying the velocity validation at CV1 in a scatter plot, u velocity and v velocity time-series graphs respectively. The results are satisfying and demonstrates that the model is reliable by nicely replicating the amplitude, phase, and direction of the measured current velocity. Nevertheless, this tidal-driven depth-averaged model is still open for more development such as running threedimensional and adding more constituent in the model, i.e. density variations (Orhan et al., 2019).

Tidal flow
Time series of water level and current velocity is presented in Figure 6, and, it is seen that there is a different trend in velocity between high and low tide. Figure 7 shows the spatial distribution of averaged peak velocity at high and low tide from June 1 -30, 2019. At high tide, water flows to the North with velocity up to 1.5 m.s -1 . While at low tide, water flows to the South with velocity reaching 1.2 m.s -1 . Both conditions indicate that the velocity in the northern region is higher than the southern region, and it is expected due to the depth characteristic. The northern region is shallower, less than 100 meters since it is a part of Sunda Shelf.
Meanwhile, water depth in the southern region is closing to 1,000 meters and going deeper in the direction to Java Trench. The velocity is amplified at the bottleneck, between Sangiang and Panjurit Islands (refer to Figure 1d.).         Figure 8 shows the distribution of maximum depth-averaged velocity each for flow to North and South. Flow to the North is seen to be more energetic, with more region having velocity over 1.5 m.s -1 . Waters around Panjurit and Sangiang Islands are especially containing current velocity over 2 m.s -1 . The previous study shared a similar finding, indicating the regions are having velocity up to 2.6 m.s -1 (Welly et al., 2012).

Conclusion
The default model provides satisfying results on water level validation; however, the velocity is inaccurate. After going through the sensitivity analysis the RMSE and r values on ADCP data improved significantly, from 0.209 to 0.093 and 0.908 to 0.920. The bathymetry, boundary conditions, and bed roughness were found to be enhancing the RMSE value of measured ADCP data for around 30%, 31%, and 6% respectively. Further, the model showed that the tide flowed to the North at high tide and vice versa. The model improvement and the findings are important to minimize the uncertainty of operational oceanography in Sunda Strait. In the future, tidal energy harvesting and construction of Sunda Strait bridge also may take place where this study will find its use. Flow around Panjurit and Sangiang Islands shows a potential with speed over 2 m.s -1 .