Tides and Their Dynamics over the Sunda Shelf of the Southern South China Sea

A three-dimensional Regional Ocean Modelling System is used to study the tidal characteristics and their dynamics in the Sunda Shelf of the southern South China Sea. In this model, the outer domain is set with a 25 km resolution and the inner one, with a 9 km resolution. Calculations are performed on the inner domain. The model is forced at the sea surface by climatological monthly mean wind stress, freshwater (evaporation minus precipitation), and heat fluxes. Momentum and tracers (such as temperature and salinity) are prescribed in addition to the tidal heights and currents extracted from the Oregon State University TOPEX/Poseidon Global Inverse Solution (TPXO7.2) at the open boundaries. The results are validated against observed tidal amplitudes and phases at 19 locations. Results show that the mean average power energy spectrum (in unit m2/s/cph) for diurnal tides at the southern end of the East Coast of Peninsular Malaysia is approximately 43% greater than that in the East Malaysia region located in northern Borneo. In contrast, for the region of northern Borneo the semidiurnal power energy spectrum is approximately 25% greater than that in the East Coast of Peninsular Malaysia. This implies that diurnal tides are dominant along the East Coast of Peninsular Malaysia while both diurnal and semidiurnal tides dominate almost equally in coastal East Malaysia. Furthermore, the diurnal tidal energy flux is found to be 60% greater than that of the semidiurnal tides in the southern South China Sea. Based on these model analyses, the significant tidal mixing frontal areas are located primarily off Sarawak coast as indicated by high chlorophyll-a concentrations in the area.


Introduction
The South China Sea (SCS) is a semi-enclosed tropical sea, located between several land-masses that include Peninsular Malaysia, Borneo, the Philippines and East Asia. The SCS has a complex bathymetry with a depth ranging from over 1000 m in the middle and northern parts to less than 100 m in the continental shelf (Fig 1). The southern South China Sea (SSCS) is bounded by Peninsular Malaysia's eastern continental shelf, the Gulf of Thailand and the sea off Borneo and the southern coast of Vietnam. It is connected to the Java Sea through the Sunda Shelf in the south. Monsoonal winds have great influence on the sea circulation in the SSCS [1][2][3][4][5][6][7][8][9][10][11]. Furthermore, sea surface and seabed also have different significant impacts on  Table 2 for station identifications), where the tidal harmonic constants are available for model verification. The dashed line marked with T1 at (b) indicates the pathway between the mouth of SSCS and the Java Sea, transect of A1A2 is used to evaluate the bi-monthly (January-February) variations of density and nutrients (such as phosphate) in the water column. (c) Bathymetry of the Strait of Malacca is marked by the dashed line T2 at the head of the Strait to estimate the resonance frequency.
wind-induced circulations as well as on the distribution, propagation and dissipation of tidal energy flux [12] with complicated tidal dynamics, especially at the bottom [13][14][15]. In this region, several successful 2D numerical studies have been performed by Ye and Robinson [16] on M 2 and K 1 tidal constituents with approximately 34 km resolution, Fang et al. [17] on M 2 , S 2 , K 1 , and O 1 with 28 km resolution. Zu et al. [18] used simulation model of 10 km resolution and ETOPO5 with a model integration of 240 days to investigate characteristics and dynamics of M 2 , S 2 , K 1 , O 1 , N 2 , K 2 , P 1, and Q 1 . Most recently, Green et al. [19] used the Oregon State University Tidal Inversion Software (OTIS) without data assimilation but with a realistic tidal conversion scheme to demonstrate that the modelled dissipation levels are overestimated over the entire SCS, noting that the discrepancies are far larger in regions with steep bathymetry.
All these studies illustrate that the existing tidal currents within the continental shelf are strong and complex. In contrast to the northern region of the South China Sea, the dynamics of the SSCS have yet to be thoroughly investigated. The lack of a detailed modelling study and sufficient observations in the SSCS region justifies the need for such a study. Hence, the aims of this study are to understand and simulate tidal characteristics and their dynamics in the region using the Regional Ocean Modelling System (ROMS), as well as to validate the model against the observed values from tide gauges (TGs). After describing the model set up, we discuss the model validation and the basic tidal features, inclusive of its dynamics and mixing fronts. The final section summarizes our study.

Model configuration
A two-domain, one-way nested model is configured for the SSCS (Fig 1). The model is based on ROMS (refer to: https://www.myroms.org/wiki/index.php/Documentation_Portal), which is commonly used for coastal applications [20], and the developed version (ROMS AGRIF) by the Institut De Recherche Pour Le Développement (refer to: http://www.romsagrif.org) is utilized for simulation here. The bathymetry of the outer and inner domains is based on ETOPO2 (refer to: http://www.ngdc.noaa.gov), which is derived from the depth soundings and satellite gravity observations [21]. However, the bathymetry is smoothed to reduce the pressure gradient error to an acceptable level by using the relative bathymetric gradient (r = rh/h) to be 0.2 [6][7]11]. The outer domain has a horizontal resolution of 0.25°×0.25°which is approximately 25 km and is vertically separated into 30 S-levels following the bathymetry but with a minimum depth (hmin) setting of 5 m at the shore. The outer domain covers 20°S to 30°N and 90°E to 140°E, thus encompassing the eastern Indian Ocean and the western Pacific Ocean. The inner domain covers 5°S to 14°N and 99.5°E to 117°E, thus encompassing the SSCS (Fig 1B). It shares the same vertical levels as the outer domain but utilizes a finer horizontal grid which is 0.083°×0.083°, each representing approximately 9 km. Overall, the dimensions are 199×207×30 and 239×207×30 for the outer and inner domains respectively.
The K-profile parameterization scheme, which includes important physics of the upper ocean mixing, is used for the vertical mixing processes [22]. The lateral boundary conditions for the inner domain are provided by integrating the outer domain according to the time step [23], thus making it into a one-way nested model. The leapfrog integration of outer and inner domains are based on the Adaptive Grid Refinement in FORTRAN (AGRIF) [24], as provided by the Institute De Recherche Pour Le Développement (refer to: http://www.romsagrif.org). It uses boundary conditions from the coarse outer domain for the fine resolution in the inner domain. The four open boundaries are specified at the north and south, east and west of both domains respectively (see Fig 1A). The simulation included both barotropic and baroclinic modes. The barotropic mode calculates the 2D momentum and elevation fields, whereas the baroclinic mode computes the 3D velocity and tracer (such as temperature and salinity) fields. Flather's [25] and Chapman's [26] boundary conditions are used for the 2D momentum and elevation fields respectively. The Orlanski [27] radiative boundary condition is used for the 3D fields.
The tidal signals are added to the primitive equation model through the outer boundaries which ignore the earth tides and astronomical tides potential flow [28]. The tidal force is computed using the TOPEX/Poseidon Global Inverse Solution (TPXO7.2) [29] with a resolution of 0.25°from the global barotropic tidal model of Oregon State University (refer to: http://volkov. oce.orst.edu/tides/global.html). Eight major tidal constituents (M 2 , S 2 , K 1 , O 1 , N 2 , K 2 , P 1, and Q 1 ) are used. Their harmonic phases (in degrees) are presented in Universal Time (UT). The phase relation between Universal Time (UT) and Local Time (LT) is UT = LT-8ω, where the desired value of ω (in units of degrees/hour) is 28.984, 30.000, 15.041, and 13.943 for the M 2 , S 2 , K 1 , and O 1 tides, respectively. Apart from the above, atmospheric and oceanic forces (such as wind stress, net heat and freshwater fluxes) from the Comprehensive Ocean-Atmosphere Data Set (COADS) [30], (refer to: http://iridl.ldeo.columbia.edu/SOURCES/.DASILVA/. SMD94/.climatology/) with a horizontal resolution of 0.5°×0.5°are also injected into the ROMS model. During the hourly integration of baroclinic tides model for the selected 90 days period (15 th December-15 th March), the model is initialized to the climatological monthly mean salinity [31] and temperature [32]

Model Validation
From the last two months run for the inner domain of ROMS, u-and v-components as well as sea surface height are used to estimate the current and elevation for each tidal constituent based on the T-Tide harmonic analysis [33]. The estimates of simulated tidal phase and amplitude are compared with the tide gauges (TGs) data. Data from 19 TGs provided by the Department of Survey and Mapping Malaysia are used to validate the modelled sea surface elevation. Seventeen of these TGs are in the Malaysian waters, mostly distributed along the east coast and the southern end of Peninsular Malaysia, as well as Borneo. The remaining two gauges (Laut Island and Natuna Islands) are located within SSCS ( Fig 1B). Tidal harmonic constants for these tidal gauges are computed and used for model validation.
The root mean square error (RMSE) differences in terms of amplitude and phase from global inverse tide model (TPXO7.2) and ROMS are computed with respect to those from tidal gauges ( Table 1).
The RMSE for amplitudes between modelled and observed is computed using the following equation: A mi and A oi are the respective modelled and observed amplitude at the station i. For the phases, the RMSE is computed by finding the pair-wise differences of the phases around the circle in degrees [34]. The inner domain of ROMS has the smallest RMSE in the estimated tidal amplitude and phase in comparison to those of the outer domain and TPXO7.2 (as evident in Table 1). SSCS is known to have complex bathymetry and coastlines. Complex bathymetry may lead to non-linear interaction between neighbouring tidal frequencies and baroclinic instabilities, which inject energy to the upper layers and modify the surface tidal signals [35].
As there are many small scale features in the region, higher resolution and smaller RMSE differences in amplitude and phase of the inner domain thus provide better simulation. More importantly, there is an additional need to accurately reproduce the stratification and baroclinicity of SSCS due to its error sensitivity towards bathymetry and model resolution. To attain this simulation, it can be demonstrated by analyzing the annual variations of average potential density and buoyancy frequency [6][7]11] in the different selected latitudes of the SSCS with respect to the following equation (Eq 2) as shown in Fig 2.
where N is the Brunt-Väisälä frequency/buoyancy frequency (Hz), g is the acceleration of gravity (m/s 2 ), and ρ the potential density (kg/m 3 ). The stratification of layers and hence its indirect assessment of vertical displacement of water parcels are very clearly shown in Fig 2. Also, the existence of the maximum frequency of oscillatory fluid in the water column could not only lead to the formation and propagation of internal waves but also amplify the sea surface waves. The simulated and TGs elevations of the East Coast of the Peninsular Malaysia (ECPM) at Chendering, Tioman Island and the coastal regions of Sabah at Kota Kinabalu, and the Labuan station located off north-western Sabah are shown in Fig 3. The comparison between the TGs surface elevations and the simulated elevations shows that the spring-neap, diurnal and semidiurnal tides of the simulated model are reasonably reproduced by the model at all stations. This implies that the boundary conditions are accurately specified. However, small differences between the TGs and the simulated results can be expected due to the complex bathymetry [19,[35][36]. At the Labuan station, the elevation is simulated reasonably well in terms of phase and particularly, the tidal range. Results indicate that the largest tidal range occurs at Tioman Island (Station 15) and Labuan (Station 19), whereas the tidal ranges at Chedering (Station 13) and Kota Kinabalu (Station 17) are lower.
The skill of the model is assessed with a detailed comparison of the semidiurnal (M 2 and S 2 ) and diurnal (K 1 and O 1 ) tidal constituents for the simulated model and TGs as listed in Tables  2 and 3. Absolute error for the amplitudes and phases is calculated based on the differences between those of the model and those of the observed tides. Some stations indicate large absolute errors for both the amplitude and phase (see Tables 2 and 3). These large absolute errors may be due to the sensitivity of the tidal constituents to the bathymetry, especially in regions of steep bathymetry where the largest turbulent dissipation occurs [19,[35][36].
The d value which is the difference between the two sets of harmonic tides (observed and modelled) estimated as the distance in the complex plane according to the following formula where A o , A m , p o , and p m represent the observed and modelled amplitudes and phases, respectively. The obvious differences for M 2 occurred at Johor Bahru (Station 11) with 47 cm and at Kukup (Station 14) with 50 cm, both at the southern mouth of the Strait of Malacca. Such difference can be attributed to the insufficient model resolution, leading to the inaccurate representation on the dissipation of energy caused by bottom friction and baroclinicity [28]. Also, the tidal waves response from the model may be very sensitive to the local bathymetry and coastline. The estimated amplitudes and phases from the simulated model for the semidiurnal (M 2 and S 2 ) and diurnal (K 1 and O 1 ) tidal constituents are compared against the corresponding values observed from the TGs (Figs 4 and 5). The estimation is done using the following linear regression equation with 95% confidence intervals, where n represents the number of stations, x i and y i are the model and observed variables at the station i. The regression coefficients of β 0 and β 1 are known as the y-intercept and the slope respectively, and e i ¼ y i Àŷ i is the error term. It can be seen that the estimated semidiurnal and diurnal tidal phases and amplitudes match reasonably with the observed values as these values lie between the confidence bounds (Figs 4 and 5). There are only two points which are significantly away from the confidence bound for semidiurnal (M 2 and S 2 ) tides. These points  Tables 2  and 3).

Basic Tidal Features M 2 tide
The spatial pattern and magnitude of amplitude and phase lags from the simulated M 2 tide are generally similar to those found by Fang et al. [17], Zu et al. [18] and Green and David [19] ( Fig 6A). The M 2 amplitude is generally higher (> 0.7 m) along the southern Vietnam coast and the north-western coast of East Malaysia due to the refraction of the waves from the coastlines. For the ECPM it is approximately 0.5 m higher in the southern coast and lower in other areas as a consequence of strong shoaling and narrowing effects from the deep basin [18].
Based on the works of Fang et al. [17], a simulated elongated nodal band is found to originate from the southern tip of the Gulf of Thailand, running parallel to the ECPM. This node ends somewhere between the southern Peninsular Malaysia and northwest coast of East Malaysia. Along this nodal band, there exits two amphidromic points located at 8.5°N, 104.5°E and 2°N, 108°E, respectively, with clockwise rotation. The simulated amphidromic point at 8.5°N, 104.5°E is noted to be almost consistent with those done by Yanagi and Takao [37], Fang et al.  [17], Jan et al. [38] and Zu et al. [18]. However, the amphidormic point at 2°N and 108°E is not featured in either Fang et al. [17] and Zu et al. [18]. The existence of amphidormic points in shallow waters is due to the fact that the bottom stress is computed using the velocity of the layer nearest to the seabed [39][40][41] in the three dimensional model instead of using uniform mean velocity as in the two dimensional model. The simulation shows that the tidal waves from the deep basin reach the Sunda shelf in the form of standing waves. Propagation speed, C (m/s), is estimated by the following equation: where L denotes the distance between two neighbouring co-phase lines, and D g= 360 the time needed for the waves to travel across this distance. Δg is the difference between two neighbouring co-phase lines in degrees and T M 2 the period of the M 2 tide. Using the values between 30°and 330°of the co-phase lines (Fig 6A) in the Sunda Shelf, Δg = 300°, T M 2 ¼ 12:4 hours, and L ffi 445 km, C in the continental shelf area is found to be~12 m/s. The wave speed for the M 2 in the SCS deep basin is approximately 164 m/s [18]. This implies that the wave propagation in the deep basin is faster than that in the continental shelf. Apart from the wave speed, the tidal currents are strong in the continental shelf but weak in the deep basin (as evident in Fig 6B). Furthermore, the large and strong currents coincide with the areas of high amplitudes. Tides over Sunda Shelf

K 1 Tide
The spatial pattern of the K 1 amplitudes is much different from that of the M 2 even though it shows similar tendency of having higher values over the continental shelf than in the deep open sea region (Fig 7A).
In comparison with the M 2 tide, two semi-circular cores of high tidal amplitudes are no longer visible off south-eastern Vietnam and western East Malaysia. High K 1 amplitudes are instead found closer to southern Sumatra in the Karimata Strait as shown in Fig 7A. As the tidal currents propagate from the deep basin into the continental shelf, they diverge partly into the Gulf of Thailand and turn counter-clockwise in contrast with that of M 2 . The contrast in rotation here may due to its specific tidal phase [37]. The tidal currents also flow southwards into the Java Sea through the Karimata Strait in similar counter-clockwise direction (Fig 7B). In general, the M 2 co-phase lines are densely clustered with smaller lines as compared with those of the K 1 tidal waves.

S 2 and O 1 Tides
The characteristics of two tidal constituents, the principal solar semidiurnal (S 2 ) and lunar diurnal (O 1 ) are shown in Fig 8. Similar co-amplitude and co-phase patterns in M 2 and K 1 are also noted in S 2 and O 1 . In terms of amplitude, the S 2 is smaller than the M 2 while the O 1 is smaller than the K 1 . These results are consistent to those obtained by previous studies Tides over Sunda Shelf [17][18]42]. Only four amphidromic points can be seen in the S 2 at 9.5°N, 105°E; 5°N, 105°E ; 3°N, 106.5°E and 1.7°S, 109°E. These points are consistent with those obtained in Mao et al. [42].

Tidal Dynamics Tidal resonance
Tidal resonance occurs when the tidal force excites one of the resonant modes of the coastal sea, especially when the continental shelf width is of a quarter wavelength. Consequently, tidal energy strengthens due to reflections between the coast and the vicinity of the shelf edge, leading to the production of a very high tidal range at the coast. Assuming that the tidal wave length from the mouth of the SSCS (see Fig 1B, dashed line (T1) with the length of x between L0-L1), i.e., the head of the channel which is the pathway between the mouth of SSCS and the Java Sea is much larger than the mean depth of the channel, shallow water equations are therefore applicable to simulate the tidal waves as described in the following equation:  wave shown in the following equation, ( where u max is the maximum tidal current, A ¼ su max gk is the wave amplitude, k ¼ 2p l is the wave number, s ¼ 2p Tl is the wave frequency, T and λ are the wave period and wave length, respectively. In terms of the corresponding wave dispersion relation (σ 2 = ghk 2 ), the wave length (λ) is ffiffiffiffi ffi gh p . If the tidal amplitude at the head of the channel (x = L 0 L) is A x, then according to a study by Bowden [43], the elevation and amplitude along the T1 can be written as η = A x cos (σt) and A = A x /cos(kx), respectively.
Resonance will occur if cos(kx) ffi 0, with the arguments of cosine (kx ¼ p 2 ; 3p 2 ; 5p 2 ; . . . ; ð2nÀ1Þp 2 , n = 1, 2,. . .). The first resonance mode (n ¼ 1 ! kx ¼ p 2 ) with the length x ¼ l itw 4 leads to the Helmholtz resonance in which λ itw is defined as the length of the incoming tidal wave. Nevertheless, a computational estimation based on the following Eq (8) which is associated with the resonant angular frequency (ω 0 ) and period (T 0 ), by assuming the length x = L 0 L ffi 500 km and average depth of the channel h ffi 50 m at the entrance of the channel, the resonance period The above estimated value is approximately equal to the periods of K 1 (23.93 hours) and O 1 (25.82 hours) tidal constituents. This implies that the amplitude of the diurnal tides increases similar to the Helmholtz resonance after the tidal wave propagates from the SSCS into the Java Sea through the Karimata Strait.

Tidal energy flux and its dissipation
The strength of the tidal constituents' variations (energy) in the continental shelf of the SSCS is assessed using power spectral density function (m 2 /s/cph). It shows that the diurnal frequencies are most dominant at the selected stations (as evident in Fig 9). The mean power energy for diurnal tides at Labuan (Station 19) is approximately 43% less than that at Chendering (Station 13), while the mean power energy for the semidiurnal tides at Labuan (Station 19) is greater by approximately 25% than that at Chendering (Station 13). This implies that the diurnal tide is more energetic along the ECPM, while both diurnal and semidiurnal tides are almost equally energetic at coastal East Malaysia.
To assess the tides and its energy dissipation in the continental shelf of the SSCS, the tidal energy flux (J = W/m) is computed based on the following equation: : where E is the total tidal energy per unit area, ρ 0 (ffi1025 kg/m 3 ) is the mean density of the seawater, h, the depth of the water column (m), u and v (m/s) are the vertically averaged zonal (semi-major) and meridional (semi-minor) components of the tidal currents, respectively. The group velocity C g (m/s) is given by; is the wave celerity in m/s, k ¼ 2p l , the wave number, λ = TC denotes the wave length (m), and T (in sec) is the period for each tidal constituent.
The general patterns of the tidal energy flux for the diurnal tides (the K 1 and O 1 constituents) are dominant as compared with the semidiurnal tides over the continental shelf of the SSCS (Fig 10). The high semidiurnal tidal energy flux of approximately 0.4×10 5 W/m is seen near the coastal regions of Sabah and Sarawak, especially near Tanjung Datu (Station 6), Kuching (Station 9) and Kuala Paloh (Station 10) which are located in the coastal areas of Sarawak (Fig 10A). This indicates that these areas are primarily influenced by the semidiurnal tides. Similarly, the strong tidal energy flux of both the diurnal and semidiurnal tides are found stretching along the coastal areas of Sabah (station Kota Kinabalu), Sarawak (station Miri), and Labuan. However, it was noted that the dissipation of the tidal waves with strong tidal energy flux of approximately 0.4×10 5 W/m tends to be at Miri (Station 8) located at the coast of Sarawak.
The diurnal tidal energy flux dissipates in the Java Sea after traversing the SSCS through the Karimata Strait (Fig 10C and 10D). Here, the mean diurnal tidal energy flux is 6% greater than the semidiurnal tides in the area between 105°E-109°E and 5°S-1°N. In general, the total mean diurnal tidal energy flux over the whole domain between 103°E-116°E and 5°S-8°N is 60% greater than that of the semidiurnal tides, implying the dominance of the diurnal tides in the region. These results are in good agreement with those previous studies which have noted that the dissipation of tides from the Pacific Ocean into the SCS is primarily diurnal [44][45].

Tidal Mixing Fronts
The mixing parameter (given by log 10 ð h = U 3 Þ) from Simpson and Hunter [46] is used to assess the contribution of tidal waves in shallow water on key biological processes that lead to the transport of nutrients and phytoplankton blooms across a tidal mixing fronts, where h (m), is the water depth and U, the tidal current amplitude (m/s).Yanagi et al. [47] shows the existence of the tidal mixing fronts in the SSCS for the northern Gulf of Thailand and the offshore area of the ECPM. Previous findings have pointed out that the regions with smaller parameter values of~<3 are locations of tidal mixing fronts [47][48][49]. Hence, coastal areas of Sarawak, namely Tanjung Datu (Station 6), Kuching (Station 9), and Kuala Paloh (Station 10), are found to be the locations of tidal mixing fronts due to the existence of strong M 2 tide because of its Tides over Sunda Shelf maximum tidal current and strong tidal energy flux. For the diurnal tides, the maximum tidal current is clearly correlated with the strong energy flux in the region between southern tip of the ECPM and western East Malaysia (Fig 10C and 10D). To distinguish the effect of tidal mixing from the other atmospheric and oceanic forces, the bi-monthly (January-February) averages of the concentration of chlorophyll-a (in mg/m 3 ) from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua for the periods 2005-2013 (http://podaac.jpl.nasa.gov/dataset/MODIS_Aqua_L3_CHLA_Daily_4km_R), which is an appropriate proxy of phytoplankton biomass along the coastline is used to assess its response to the various corresponding forces. It is found that areas with mixing parameter values~< 3 are well correlated with the areas of strong tidal current and high tidal energy flux which are supported by the high concentrations of chlorophyll (as evident in Fig 11).
Tidal mixing processes for this region can be assessed from the distributions not only of the bi-monthly wind stress curl (derived from COADS [30]) and simulated density difference between the sea surface and its subsurface, but also of the potential temperature (y ¼ T þ R P 0 P GðS; T; PÞdP which is a function of Pressure (P), where Γ is adiabatic lapse rate, S and T are salinty and temperature, respectively) as computed from the model. Accordingly, both low values of density difference and potential temperature upon the area of zero curl (implying no Ekman transport) (Fig 12A and 12B) in the vicinity of coastal Sarawak indicate that this area is well mixed. Hence, the buoyant tidal fronts become significant and enhance the instability (see Fig 12C, given by E ¼ N 2 g $ À1 r @r @z ). This causes the cold water at the sea surface to induce vertical convection and subsequent mixing to a relatively deep mixed layer depth (as evident in Fig 12D) based on calcualtion by Lorbacher et al. [50]. The density profile along the transect A1A2 (see Fig 13A) reveals that denser water is uplifted to the sea surface in Tides over Sunda Shelf conjunction with the well-mixed and uniform concentration of nutrients (such as phosphate) from the sea surface to an approximate depth of 25 m (see Fig 13B). It is also to be noted that this area has low incoming shortwave radiation (Fig 14). Therefore, the high chlorophyll concentration at the sea surface in the vicinity of coastal Sarawak is the outcome of the M 2 tidal mixing front generated by the high tidal energy flux there. In contrast, along the western coast of Borneo, the high incoming shortwave radiation (see Fig 14) warms the upper water layers, leading to the reduction of mixed layer depth (as evident in Fig 12D) and decreasing temperature below it. The subsurface cooling below the warmed mixed layer reduces vertical mixing in the upper thermocline and redistributes nutrients apart from the heat and salt [51].

Tides over Sunda Shelf
It is interesting to see that the chlorophyll concentration in the Strait of Malacca is very prominent too. From Eq (8), the resonance frequency at the head of the Strait of Malacca (see Fig 1C, dashed line: T2) can be estimated to be around 11.3 hours. This value is roughly close to the periods of M 2 (12.4 hours) and S 2 (12 hours) due to the narrow channel (approximately 350 km width and 120 m deep) feature here. As such, very strong mixing is effected by semidiurnal tides, mostly by M 2 and to a lesser extent by S 2 tides. Further work on this aspect will be done in future.

Summary and Conclusions
Results of the simulated tides in the Continental Shelf Area of the SSCS from a 3D, one-way nested regional ocean modelling system indicate that the modelled tides compare well with the observations at 19 tidal stations. Existence of high tidal elevation at the southern tip of the ECPM and East Malaysia reflect the significant role of tides in these regions. Moreover, the M 2 tidal current tends to flow towards ECPM due to the convergence of the counter clockwise tidal current from the coast of Vietnam and the clockwise tidal current from East Malaysia. Depth of the bathymetry can alter significantly the phase change. In general, the region between the southern end of the ECPM and coastal Sarawak is dominated by the diurnal tide with a counter clockwise phase rotation. However, along the coastal areas of East Malaysia, especially the coastal areas of Sarawak in places like Tanjung Datu, Kuching and Kuala Paloh are affected by the M 2 tide with the high energy flux and strong tidal current rotating in a clockwise direction. The patterns of the co-tidal lines (lines of constant tidal phase) for the two main tidal constituents, namely M 2 and K 1 are the result of their differing amplitudes, leading to different contributions of tidal energies for each tidal constituent in the SSCS. The largest tidal energy flux in the study area is related mostly to the diurnal tides. There is a good correlation between the areas of high tidal energy flux and strong tidal current, causing the areas to be well mixed by the tidal mixing fronts. The areas affected significantly by the tidal mixing fronts are the coastal areas of Sarawak as exemplified by high chlorophyll concentrations. Nonetheless, more work needs to be done by considering the use of more refined model and bathymetry resolutions as well as more robust harmonic analysis with favourable hydrographic conditions.