http://journal.uir.ac.id/index.php/JGEET E-ISSN : 2541-5794 P-ISSN :2503-216X Journal of Geoscience, Engineering, Environment, and Technology Vol 5 No 1 2020 Ondara, K., et al./ JGEET Vol 5 No 1/2020 25 RESEARCH ARTICLE Hydrodynamics Features and Coastal Vulnerability of Sayung Sub- District, Demak, Central Java, Indonesia Koko Ondara1,*, Ruzana Dhiauddin1, Ulung J. Wisha1, Guntur A. Rahmawan1 Research Institute for Coastal Resources & Vulnerability, Ministry of Marine Affairs and Fisheries, Indonesia. * Corresponding author : kornkoko@yahoo.com Tel.:+628976929648 Received: Sept 19, 2019; Accepted: March 20, 2020. DOI: 10.25299/jgeet.2020.5.1.3996 Abstract The Sayung sub-district is an abrasion area in Demak Regency that is mostly affected by sea level rise. The purpose of this research is to determine the features of hydrodynamics and coastal dynamics occurrence in the sub-district of Sayung.Collecting field data/information and modeling approach (tides, waves, currents, weather and coastline changes) have been done in Sayung, Demak. The wave height in the eastern coast is the highest. The significant wave height in 2004 was greater than March 2016 showing that in 2004 the wind energy transfers were larger than 2016. The refraction coefficient in 2016 for all directions was the greatest from the west at the depth of 8 m and the smallest one was identified in the south. The refraction coefficient in 2004 for any direction yielded the largest value in the southwest at the depth of 2 m and the smallest onewas identified the south as well. During a cycle of tidal fluctuation, it occurs twice flood and ebb events. The maximum depth is 6.5 m located about 3.8 km from the coastline. The sediment thickness reached 564,886.39 m3. Coastline changes occurred in 2003 andstarted to gain sedimentation in 2015. Data and information produced can be useful as a basisfor further developments to mitigate abrasion andto create policy-briefin managing coastline affected abrasion even thoughsomeimprovement efforts have been made. Keywords: hydrodynamics, abrasion, sedimentation, coastal vulnerability 1. Introduction Coastal development is taking place around the North Coast of Java (Pantura)resulted in changed water conditions and coastlines around the north coast which is greater compared to the south coast. Threats including tsunamis, tectonic earthquakes, the dynamics of large waters, and others will potentially take place. In addition, the north coast of Java region is often used as the center of government and community activities such as Jakarta, Cirebon, Tegal, Pekalongan, Semarang, Surabaya, and others. Sayung is one of the sub-districts in Demak Regency which experienced the worse tidal flooding compared to other sub- districts.The coastal areas in Sayung sub-district affected by tidal flood reaching 0.25 meters are Sriwulan, Surodadi, Bedono and Timbulsloko Village (Desmawan and Sukamdi, 2012a). A total of 83 slums affected by abrasion and almost 2,036 slums affected by inundation (Marfai, 2012). the highest tidal flood occurs in the west season (Widada, Sugeng, Baskoro Rochaddi, 2012). During 1999-2009 in Demak Regency, Sea level rose around 0.72 mm/year and 7.9 mm/year in relative sea level rise. The sea level rise inundated an area of 26.83 kmΒ² in 8 villages. Significant coastline changes in Wedung District resulted in the maximum accretion of 233,994 meters and maximum abrasion of 141,037 meters during 2012 - 2015 (Dewi et al., 2016). This shoreline changes causes the settlements to draw closer to the sea, resulting in increased coastal vulnerability due to tidal flood events(Wahyudi et al., 2012). Abrasion that occurs in Sayung subdistrict is influenced by waves that occur continuously (Ondara and Wisha, 2016) with the type of soil composed of silt and sandy silt and silt (Subardjo, 2004) causing the population of people living in the area to decrease (Damaywanti, 2013). Due to frequent tidal floods and abrasion, the Sayung Subdistrict community supplies fresh water for daily activity from other areas (Desmawan and Sukamdi, 2012). Fig. 1. Administration Map of Demak, Central Java (Source: Bappeda Demak,2010) This study aims to review the process of hydro- oceanography and coastal dynamics in Sayung, Demak to determine the coastal vulnerability. Data/information and model of this research can be useful as a basis in creating policy brief related to the further developments in the future. 2. Material and Method Java sea semarang Semarang regency Grobogan Regency Pati Regency Kudus Regency Jepara Regency http://journal.uir.ac.id/index.php/JGEET 26 Ondara, K., et al./ JGEET Vol 5 No 1/2020 The research was conducted during March - August 2016 in the Sayung coast which is administratively located in four villages that are Sriwulan, Bedono, Timbulsloko, and Surodadi. Those four villages were chosen because of its susceptibility to the changes in settlement areas caused by abrasion and tidal flood inundation. This research begun with a comprehensive literature study from various previous-related research sources such as scientific journals, technical reports, proceedings and research reports carried out to identify coastal issues in the north coast of Java especially in the Sayung coast and surrounding. The secondary data and information were also employed consist of several aspects of GIS, oceanography and coastal vulnerability. Furthermore, data collection was conducted directly in the field to retrieve primary data as the basis for analysis and to validate the potential coastal resources and its threatening factors. Characteristics of wind including magnitude and direction was analyzed. Wind features were measured using an anemometer, a tool consisting of a spinning mug or blade. The average wind speed over a given period of time is calculated by counting the number of rounds over a given time period. Conversion of wind speed can be made using Shore Protection Manual 1984 or 1993 API formulation as follow: π‘ˆπ‘‘ π‘ˆ3600 = 1.2 + 0.2 β‹… π‘‘π‘Žπ‘›β„Ž {0.9π‘™π‘œπ‘”10 ( 45 𝑑 )} ; 0< t < 3600 sec. (1) π‘ˆπ‘‘ π‘ˆ3600 = 1 + {3.0 + 𝑙𝑛 ( 3 𝑑 ) 0.6 }{0.15( 1 2 ) βˆ’0.125 } ; t ο‚£ 60 sec. (2) Where: Ut = Wind speed with a duration of more than 3 seconds U3600 = Wind speed with a duration of 1 hour Wind data were used to hind cast significant wave heights and periods that will be corrected against elevation, stability, location effects and drag coefficient to obtain wind stress factor (UA). Wind data used were the maximum monthly wind data that will induce the maximum wave height. Wind stress factor (UA) was calculated using following equation: π‘ˆπ‘ˆ = 0,71π‘ˆπ‘ˆπ‘ˆ1,23 (3) π‘ˆπ‘ˆ Is the wind speed in the sea obtained from the equation: π‘ˆπ‘ˆ = π‘ˆπ‘ˆπ‘ˆπ‘ˆπ‘ˆ (4) whereπ‘ˆπ‘ˆ is the measured wind velocity and π‘ˆπ‘ˆ is the relationship between winds above sea and wind in the nearest land. To determine the maximum wave height and period that occurs basedon the depth change from 0.1 m to 10 meters, a relationship between maximum wave height and maximum wave period (10 years of input data) has been calculated using hindcasting method to obtain the equation from the correlationgraph which then used to calculate the period of wave in the deep ocean. To calculate the wavelength of the deep sea used a formula as follow: π‘ˆ0 = 1,56 π‘ˆ2 (5) where L0is length between wave crests or troughs and T is time it takes a wave crest to travel one wavelength (units of time). The value of Refraction Coefficient (Kr) is obtained using the equation as follow: Kr = √ cos Ξ±0 cos Ξ± (6) The coefficient of shoaling is calculated by using the equation of: Ks = √ n0L0 nL (7) where Kr is the refraction coefficient, and Ks is the coefficient of shoaling, 𝛼 is angle of incidence relative to the shoreline or some prescribed depth contour and n values varies from 0,5 for a deep-water wave up to 1.0 for a shallow-water wave. 3. Result and Discussion To facilitate the reading of wind data for 10 years, the wind rose diagram was created shown in Figure 2. Fig 2. wind rose in 2005 - 2015 The direction information of incoming wave is obtained by classifying the 2005-2015 wind data and grouped by the Beaufort scale. It can be seen in Figure 2 that the dominant and maximum waves are sourced from the Southeast. This is because the eastern region is open water and topographically there is no obstacle factor hampering the wind blows and energy transfers. Fig 3. The maximum daily wave chart on March 2016 Fig 4. Correlation graph of wave height and wave period for east direction After calculating the coefficient of refraction and shoaling, we obtained the new wave height value by using the following equation: π‘ˆ = π‘ˆ0.π‘ˆπ‘ˆ.π‘ˆπ‘ˆ (8) where H is representing the wave height at some coastal location, and Hois the wave height at some offshore locations. Ondara, K., et al./ JGEET Vol 5 No 1/2020 27 Table 1. Shoaling coefficient calculation z (m) 6 4 2 0.1 n 0.4871 0.4194 0.2028 0.9258 Ks 1.0129 1.089 1.5231 1.1464 H 0.042 0.0451 0.0631 0.0475 Table 2. Breaking wave calculation z (m) 6 4 2 0.1 H'o 0.0414 0.0414 0.0414 0.0414 H'o/gt2 0.0023 0.0023 0.0023 0.0023 According to the calculation of shoaling and refraction coefficients, it is obtained that the wave breaks at 0.04 meters with a depth of 6 meters. 3.1. Bathymetry Bathymetry data was measured directly using Echo Sounder. The depth of the horizontal waters is processed into the depth contours of the seafloor topography to show the depth variation. The maximum depth is 6.5 m which at about 3.8 km from the coastline. The result of bathymetry must be corrected using the existing the sea level adjusted with measurement time and the distance between surface and transducer sensor in order to obtain the actual depth (Wisha and Heriati, 2017). Fig. 5.Water Depth Contour Map 3.3 Tide Fig 6. Tidal Data Verification As one of the prerequisites of modelling, tidal data validation was done by comparing field data and model data (Gao et al., 2009). as in Figure 6, by using the root mean square error formula it was obtained the RMSE of 0.3.Figure 6 also shows that between the two surface elevation data, the same phase during the spring and neap tidal conditions is identified which is slightly different at the neap period might be caused by the influence of wave which affects the surface elevation formed (Garrett and Kunze, 2007). 3.2. Sea current The result of current field data processing is depicted in the form of scatter plot, stick diagram, current rose, current velocity vertical profile, and hydrodynamic map which illustrate the features of currents in various depths. The current speed ranges from 0.06-0.5 m/s which predominantly moves towardsoutheast and northwest.According to (Wisha et al., 2015), surface currents are influenced by monsoon winds. During the 1st transitional season, where the transition of the rainy season to the dry season took place, the wind blows from Australia to Asia. This theory corresponds to the current rose that the direction of the surface current shows the same pattern with the wind dominant direction. Fig 7. The direction of the dominant current profile Fig 8. Vertical east velocity current profile to depth The current moves southeastward and northwestward with a speed of 0.056-0.44 m/s. Furthermore, the currents patterns vertically vary at each depth. Figure 8, 9 and 10 show In the current vertical profile results explaining that the deeper the water column, the weaker velocity observed, this is due to the influence of bottom friction and the role of density (Wisha et al., 2016). During the high spring tidal phase, the current speed ranged from 0 - 0.05 m/s moving landward. While during the low tidal condition, the current flowed seaward with the speed of 0 - 0.08 28 Ondara, K., et al./ JGEET Vol 5 No 1/2020 m/s. This condition shows that the current speed id higher during ebb tides. This event is closely related to the transport mechanism in the coastal area where the water mass flowing toward the coast will be forwarded to the bottom of the waters (downwelling) and vice versa for the ebb event, where the water mass below the water level will rise to fill the void area in the surface. Fig 9. Vertical up velocity current profile to depth Fig 10. Vertical north velocity current profile to depth The deeper the water depth, the greater the density and the smaller the current speed. The movement of the insignificant currents is caused by the basic friction and density effects so that the energy generated by the ocean currents becomes weak which the inhibition of current movement makes it calmer. Thus, the pattern and velocity of the current are more randomly and faster in the surface where the obstacle is reduced and the wind-energy transfer is maximum. According to Wisha et al. (2016) the vertical profile of current velocity components is influenced by bottom topography, the dynamics of sea level pressure, bottom friction, density and winds. In the coastal areas of Sayung, there are several permeable dams in the form of piles of wood and stone. The beach building was built with the aim to reduce the impact of abrasion by supporting mangrove growth around the coastal area (Pribadi et al., 2010). The occurrence of vector accumulation around the permeable dam during high and low tidal conditions is shown in Fig. 12. Due to the narrow channel of permeable dam, the current flow passing the structures becomes more rapidly. However, the current velocity becomes weaker among the permeable dam. It can be seen from the graph of the current velocity in Fig. 11. The current profile will be calmer around the structures and more rapidly in the estuary and tend to be stable offshore (Cummins et al., 2012). Fig 11. Graph of ocean current velocity: front of permeable (black); Between permeability and river (blue); Between permeable (green); River mouth (light blue); Open sea (red) Based on the observation, during the high tidal condition, the oscillation of current toward the land and sea will bring the sediment into the permeable structure. However, during the low tidal condition, the water mass movement tends toward the sea, in this condition the sediment material will be trapped in the permeable structure. Thus, the greater the mass of water entering the permeable structure, the larger the sediment material will betrapped. 3.4 Sedimentation 3D map of field bathymetry measurements consists of channel 1 (high frequency) and channel 2 (low frequency), which was analyzed to obtain sedimentation volume in the study area. The result of sedimentation volume calculation (Fig. 14) shows the suspension of suspended sediment is high within the bay. Tidal current regime has a role in evoking the suspended sediment concentration. At the high tidal condition, the sediment concentration in the estuary will intensely accumulated. The cohesive sediment in Sayung that is highly susceptible to salinity contributes to the higher suspended sediment inducing sedimentation. If salinity becomes high, the clumping process will take place, so sediment tends to be deposited(Marfai, 2012). The result of overlapping Channel 1 and Channel 2 results in the figure of sediment thickness during bathymetry measurement. Volume calculation was performed using the surfer software which the total volume of sediment is shown in Table 3. Table 3. Comparison of Sediment Thickness Volume Equation Volume (m3) Trapezoidal Rule 573.922.20 Simpson's Rule 564.886.39 Simpson's 3/8 Rule 587.137.34 Ondara, K., et al./ JGEET Vol 5 No 1/2020 29 Fig. 12. Simulation of ocean currents at high tide Fig. 13. Simulation of ocean currents at low tide Fig 14. 3D bathymetry measurement map In Sayung Sub-district there are several rivers that carry suspended solids from the mainland. At the high tidal condition, the water masses moved from the sea toward the interior of the river, the mass of water will bring suspended solids from the sea into the river. During the spring tidal phase, the concentration of suspended sediment will be higher due to the higher astronomical energy transfer evoking the higher current speed and transport mechanism. The difference in surface elevation will cause massive movement of water masses. This is in line with the results of the study by Wisha et al. (2015) which stated that the amount of sediment will increase during the low tidal condition because the sediments carried during flood tides will be settled down in the surface bottom. 3.5 Geographic information system analysis The Fig. 15 below shows the results of satellite image map analysis of 1999 and 2003 based on Landsat-7 image obtained from USGS online data of USA. In the picture above shows that the coastline of 2003 experienced a change that leads to coastline retreat as seen in 30 Ondara, K., et al./ JGEET Vol 5 No 1/2020 the boxes 1 and 2 in Fig. 2 above. The rate of change in these two regions can be seen in Table 4 below. Fig 15. Coastline overlay 1999 (red line) and 2003 (blue line) Table 4. Coastline change between 1999 and 2003 No. The pace of change (m/year) information 1 6.83 Abrasion 2 10.85 Abrasion Total 17.68 Abrasion 3.6 Shoreline changes in 1999 – 2015 The Fig. 16 below shows the results of satellite image map analysis in 1999 based on Landsat-7 image and Landsat-8 2015 image obtained from USGS online data of USA. Overlay of the coastline of 1999 and 2015 over Landsat 7 satellite imagery (1999) shows several locations that have undergone significant changes over the past 16 years. White box shows a change of accretion while the blue box shows a shoreline retreat (abrasion). Table 5.Coastline change between 1999 and 2015 No. The pace of change (m/year) information 1 42.86 Abrasion 2 57.20 Abrasion 3 104.02 Abrasion Total 89.68 Abrasion In general, from the two tables above, the rate of shoreline changes leads to the abrasion phenomenon, where coastlines tend to retreat from year to year. Coastlines of 1999, 2003 and 2015 extracted from the imagery are then overlaid so that a map can be generated as shown in Fig. 17. Fig. 16. Shoreline Overlay 1999 (red) and 2015 (orange) Fig. 17. Map of coastline change in 1999, 2003 and 2015 From the picture above, it can be seen that most of the coastline continues to decline from year to year. However, there are also several accretion sites seen in the white circle on the map above, where coastline retreats occurred in 2003 but again progressed in 2015. Changes that occur in this coastline can be caused by natural conditions such as sediment type of water base; aquatic dynamics; climate and weather and human Ondara, K., et al./ JGEET Vol 5 No 1/2020 31 activities such as land clearing, and over-exploitation of water resources. 4. Conclusion The dominant and maximum waves that occurred in Sayung District came from the southeast, northwest and north. The dominant current direction moved toward Southeast and Northwest. The wave is greater in the east part of study area. The significant wave height in 2004 was higher than on March 2016, indicating that in 2004 the wind blew larger than 2016. The largest refraction coefficient in 2016 is identified in the west at depth of 8 m and the smallest one is in the south. While, the largest refraction coefficient in 2004 in the southwest at the depth of 2 m and the smallest one is in the south. The maximum depth is 6.5 m around 3.8 km from the coastline. The sediment thickness reached 564,886.39 m3. Coastline retreat occurred in 2003 but again progressed in 2015. Changes that occur on this coastline can be caused by natural conditions such as sediment types baseline water; Aquatic dynamics; Climate and weather and human activities such as land clearing, and over- exploitation of water resources. Coastal protection buildings that can be used to protect the coastal areas of Sayung sub-district arebreakwater, revetment, groynes and combinations. Data/information and model of this research can be useful as a basis in creating policy brief related to the further developmentsin the future. Acknowledgement Gratitude is given to Research Institute for Coastal Resources and Vulnerability (RICRV) for research funding in 2016 and for those who assisted in the completion of this work. References Cummins, S.J., Silvester, T.B., Cleary, P.W., 2012. Three- dimensional wave impact on a rigid structure using smoothed particle hydrodynamics. Int. J. Numer. Methods Fluids. https://doi.org/10.1002/fld.2539 Damaywanti, K., 2013. Dampak Abrasi Pantai terhadap Lingkungan Sosial (Studi Kasus di Desa Bedono , Sayung Demak). Pros. Semin. Nas. Pengelolaan Sumberd. Alam dan Lingkung. 363–367. Desmawan, B.T., Sukamdi, S., 2012a. Adaptasi Masyarakat Kawasan Pesisir Terhadap Banjir Rob Di Kecamatan Sayung, Kabupaten Demak, Jawa Tengah. J. Bumi Indones. Desmawan, B.T., Sukamdi, S., 2012b. Adaptasi Masyarakat Kawasan Pesisir Terhadap Banjir Rob Di Kecamatan Sayung, Kabupaten Demak, Jawa Tengah. J. Bumi Indones. 1. Dewi, R.S., Bijker, W., Stein, A., Marfai, M.A., 2016. Fuzzy classification for shoreline change monitoring in a part of the Northern coastal area of Java, Indonesia. Remote lens. https://doi.org/10.3390/rs8030190 Gao, J., Lan, X., Fan, Y., Chang, J., Wang, G., Lu, C., Xu, C., 2009. CFD modeling and validation of the turbulent fluidized bed of FCC particles. AIChE J. https://doi.org/10.1002/aic.11824 Garrett, C., Kunze, E., 2007. Internal Tide Generation in the Deep Ocean. Annu. Rev. Fluid Mech. https://doi.org/10.1146/annurev.fluid.39.050905.11022 7 Marfai, M.A., 2012. Preliminary assessment of coastal erosion and local community adaptation in sayung coastal area, central java -indonesia. Quaest. Geogr. https://doi.org/10.2478/v10117-012-0028-2 Ondara, K., Wisha, U.J., 2016. Simulasi Numerik Gelombang (Spectral Waves) Dan Bencana Rob Menggunakan Flexible Mesh Dan Data Elevation Model Di Perairan Kecamatan Sayung, Demak. J. Kelaut. Indones. J. Mar. Sci. Technol. https://doi.org/10.21107/jk.v9i2.1694 Pribadi, K., Diposaptono, S., Agung, F., 2010. Climate change and coastal zone management in Indonesia: example of adaptation at Demak Coast, Java Island, Indonesia. Res. Publ. Serv. 1, 185–200. Subardjo, P., 2004. Studi Morfologi Guna Pemetaan Rob di Pesisir Sayung, Kabupaten Demak, Jawa Tengah. Indones. J. Mar. Sci. 9, 153–159. https://doi.org/10.14710/ik.ijms.9.3.153-159 Wahyudi, S.I., Ni’am, M.F., Le Bras Gilbert, 2012. Problems, Causes and Handling Analysis of Tidal Flood, Erosion and Sedimentation in Northern Coast of Central Java: Review and Recommendation. Int. J. Civ. Environ. Eng. Widada, Sugeng, Baskoro Rochaddi, and H.E., 2012. Pengaruh Arus Terhadap genangan Rob Di Kecamatan Sayung, Demak. Bul. Oseanografi Mar. 1, 31–39. https://doi.org/10.14710/buloma.v1i2.11218 Wisha, U.J., Heriati, A., 2017. Bathymetry and Hydrodynamics in Pare Bay Waters During Transitional Seasons (SeptemberOctober). Omni-Akuatika. https://doi.org/10.20884/1.oa.2016.12.2.98 Wisha, U.J., Husrin, S., Prasetyo, G.S., 2016. Hydrodynamics of Bontang Seawaters: Its Effects on the Distribution of Water Quality Parameters. ILMU Kelaut. Indones. J. Mar. Sci. https://doi.org/10.14710/ik.ijms.21.3.123-134 Wisha, U.J., Husrin, S., Prihantono, J., 2015. Hydrodynamics Banten Bay During Transitional Seasons (August- September) (Hidrodinamika Perairan Teluk Banten Pada Musim Peralihan (Agustus–September)). ILMU Kelaut. Indones. J. Mar. Sci. https://doi.org/10.14710/ik.ijms.20.2.101-112 Β© 2020 Journal of Geoscience, Engineering, Environment and Technology. All rights reserved. This is an open access article distributed under the terms of the CC BY-SA License (http://creativecommons.org/licenses/by- sa/4.0/). http://creativecommons.org/licenses/by-sa/4.0/ http://creativecommons.org/licenses/by-sa/4.0/ http://creativecommons.org/licenses/by-sa/4.0/