**Sensitivity of Storm-Induced Hazards in a Highly Curvilinear Coastline to Changing Storm Directions. The Tordera Delta Case (NW Mediterranean)**

#### **Marc Sanuy \* and Jose A. Jiménez**

Laboratori d'Enginyeria Marítima, Universitat Politècnica de Catalunya·BarcelonaTech, c/Jordi Girona 1-3, Campus Nord ed D1, 08034 Barcelona, Spain; jose.jimenez@upc.edu

**\*** Correspondence: marc.sanuy@upc.edu; Tel.: +34-93-401-7392

Received: 26 March 2019; Accepted: 8 April 2019; Published: 10 April 2019

**Abstract:** Extreme coastal storms, especially when incident in areas with densely urbanized coastlines, are one of the most damaging forms of natural disasters. The main hazards originating from coastal storms are inundation and erosion, and their magnitude and extent needs to be accurately assessed for effective management of coastal risk. The use of state-of-art morphodynamic process-based models is becoming standard, with most being applied to straight coastlines with gentle slopes. In this study, the XBeach model is used to assess the coastal response of a curvilinear sensitive deltaic coast with coarse sediment and steep slopes (intermediate-reflective conditions). The tested hypothesis is that changes in wave direction may cause large variations in the magnitude of storm-induced hazards. The model is tested against field data available for the Sant Esteve Storm (December 2008), obtaining an overall BSS (Brier Skill Score) score on the emerged morphological response of 0.68. Later, the 2008 event is used as baseline scenario to create synthetic events covering the range from NE to S. The obtained results show that storm-induced hazards along a highly curvilinear coast are very sensitive to changes in wave direction. Therefore, even under climate scenarios of relatively steady storminess, a potential shift in wave direction may significantly change hazard conditions and thus, need to be accounted for in robust damage risk assessments.

**Keywords:** XBeach; inundation; erosion; BSS; Sant Esteve 2008; coarse sediment

#### **1. Introduction**

The impact of extreme storms on the coast is one of the costliest forms of natural disasters (Kron [1]; Bertin et al. [2]). In heavily urbanized coastal areas, such as the Mediterranean (in general) and the Catalan coast (in particular), where properties, infrastructures and businesses are located close to the shoreline, this kind of event usually results in the damage or destruction of exposed assets (Jiménez et al. [3]). These effects are the integrated consequences of two main storm-induced coastal hazards: inundation and erosion. In this context, an accurate assessment of the magnitude, location and extension of these hazards is becoming an essential part of the risk management process (e.g., Ciavola et al. [4,5]; Van Dongeren et al. [6], Jimenez et al. [7]; Plomaritis et al. [8], Harley et al. [9]) and, in this sense, the use of process-oriented models to forecast storm-induced morphodynamic changes under given scenarios is now standard (e.g., Roelvink et al., [10]; McCall, [11]; Van Dongeren et al. [12] and references therein, Dissanayake et al. [13]). Most of the studies on testing state-of-art morphodynamic process-based models have addressed cases characterized by straight coastlines and gentle slopes (i.e., conditions close to the comfort zone of the models) (e.g., McCall, [11]; Harter and Figlus [14]). However, applications to estimate costal hazards in hihgly curvilinear environments (e.g., deltaic cuspate coasts) have seldom been tested (e.g., Roelvink et al., [15]; Valchev et al., [16]; and Dissanayake et al. [13]). Furthermore, the effect of testing models based on surf-beat (i.e., the infragravity wave

band) on steep slopes and coarse sediment has been recently undertaken mainly in 1D applications (e.g., Vousdoukas et al. [17]; Elsayed and Oumeraci [18]) but rarely so in fully 2DH (2-dimensional, depth-averaged) simulations.

Within this context, the magnitude of storm-induced hazards on a highly curvilinear coast by using XBeach is assessed in the present study. The relevance and main aim of this work is twofold: first, from a general standpoint, to test the use of Xbeach on a highly curvilinear coast characterized by coarse sediment reflective beaches, and second, from the local standpoint, to analyze the sensitivity of an already identified hotspot, the Tordera Delta (NW Mediterranean) (Jiménez et al. [7]), to assess storm impacts for different storm direction scenarios. Thus, the largest recorded storm in the area is used as base case scenario. It occurred in December 2008 and had the typical incoming direction of current climate conditions where eastern (E) incoming storms dominate (e.g., Mendoza et al. [19]). Existing storminess projections under climate change scenarios for the Western Mediterranean do not predict any increase in wave height (e.g., Lionello et al. [20]; Conte and Lionello [21]), but some projections identify potential changes in wave direction (Cases-Prat and Sierra [22,23]). Due to this and to the great sensitivity of cuspate coastlines to wave direction resulting from their curvature (e.g., Slott et al. [24]; Johnson et al. [25]), the study aims to assess the potential effects of changing wave direction on extreme storm-induced hazards for the Tordera Delta. The hypothesis to be tested is that changes in wave direction may cause large variations in the magnitude of storm-induced hazards. Other studies have included the sensitivity to incoming storm direction in their assessments, such as those by Mortlock et al. [26] in Australia, or de Winter and Ruessink [27] in Holland.

The article is arranged as follows: the second section introduces the study site and the data used, describes the Sant Esteve 2008 event, which is used as the base case storm-scenario, and presents the methodological part, i.e., the used morphodynamic model and the comparative assessment framework descriptions; the third section presents obtained results; and finally, the discussion and concluding remarks are presented in the fourth section.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Tordera Delta is located in the NW Mediterranean Sea approximately 60 km northwards of Barcelona (Figure 1). It is a coarse sand delta (i.e., sediment in the range between 0.8 and 1.6 mm) with an aerial surface of approximately 4.2 km<sup>2</sup> at the end of a small river basin of approximately 879 km2 (Vila and Serra [28]). The Tordera river is dry during most part of the year, with long dry summers and episodic discharges after heavy rainfall (Martin-Vide and Llasat [29]). Coastal storms and heavy rainfall events are usually uncorrelated at the area and, in fact, during the simulated storm, significant river outflow was present after the storm peak, with a phase delay of 1 day (Sanchez-Vidal et al. [30]). It has a cuspate configuration with two well differentiated parts at each side of the river mouth (Figure 1). The northern part is fronted by the S'Abanell beach, which is oriented towards the E. It has a steep nearshore bathymetry without any relevant geomorphological features in shallow waters The northern hinterland presents a higher topography, except for the zone closest to the river mouth outlet. The southern part is fronted by the Malgrat de Mar beach and is oriented towards the S, which serves as natural protection of a lower hinterland. Shallow water bathymetry is characterized by the presence of a longshore bar running parallel to the coast from the river mouth to the SW, which encloses a plateau at 4 m water depth.

The delta coastline has been eroding during the last several decades as the net result of the littoral drift and the decrease of the Tordera river sediment output, with maximum measured retreats of approximately 120 m Jiménez et al. [31]. The combination of a progressively narrowing beach protecting a low-lying hinterland and this being mainly occupied by campsites makes this area a hotspot for storm-induced hazards (Jiménez et al. [7]) with different consequences depending on storm characteristics and beach morphology at the time of impact (see e.g., Sanuy et al. [32]). As Jiménez et al. [31] noted, under former accretive-stable deltaic conditions only extreme storms were able to exceed the capacity of protection provided by wide beaches, but under present medium/long-term erosion conditions, smaller storms have become able to exceed the dissipation capacity of the narrower beaches, increasing the frequency of storm-induced problems (e.g., Jiménez et al. [3]).

**Figure 1.** Tordera Delta study area, location and model set-up. (**a**) SWAN (Simulating Waves Nearshore) grids (red) with boundary wave-spectrum nodes (yellow), wave buoys (blue triangles), and significant wave height wave rose for the period 1948–2009 (Global Ocean Waves, GOW, Reguero et al. [33]) at the Tordera wave buoy location. (**b**) XBeach set-up, curvilinear grid domain (red) and interest areas (orange) and pre-storm topobathymetry.

Wave climate at the NW Mediterranean is characterized in the wave rose in Figure 1 (Global Ocean Waves, GOW dataset hindcast, Reguero et al. [33]). Storm events are defined as events with significant wave height (Hs) exceeding 2 m during a minimum of 6 h (Mendoza et al. [19]). The main incoming wave direction at the site is NE and E, with some events arriving from the S, especially during spring. Within this two main groups, some residual events, can also be found. Thus, in propagated conditions (nearshore area) storm conditions can be characterized with waves in the range NE-S being the two extremes of the range the most frequent situations. Nonetheless, some studies rise concerns at the Catalan coast about future climate-induced changes in storm direction, particularly a frequency transfer from the current main cluster (NE-E) towards the secondary one (S) (Casas-Prat et al. [22,23]). The area is micro-tidal, with storm surges having a limited role in storm-induced processes due to its relative weight when compared to the wave component. Notably, surges are uncorrelated with waves, are most frequently under 25 cm, with some extreme events showing maximum recorded surges around 50 cm (Mendoza and Jiménez [34]).

#### *2.2. Data*

The wave data used in this study to characterize the Sant Esteve storm and validate the models were measured by a directional wave buoy located off the Tordera Delta at approximately 70 m water depth (Figure 1) belonging to the XIOM (Xarxa d'Instruments Oceanogràfics i Meteorològics) network, which is no longer operative (Bolaños et al. [35]). The Tordera Buoy was a Datawell waverider, sampling during 20' every hour, and thus providing sea states and statistics with an hourly time-step. It covers the period from 1984 to 2013.

Wave spectra from deepwaters and wind fields used to force the model chain were provided by Puertos del Estado, from the WAM (WAve Model, version 4, European Centre for Medium-Range Weather Forecasts, Reading, UK) and HIRLAM-AEMET (HIgh Resolution Local Area Modelling version 7, Agencia Estatal de METeorología, Madrid, Spain) models respectively. Both datasets have a temporal time-step of 1 h. Nearshore water levels were obtained from the same source, provided by the HAMSOM model (HAMburg Shelf Ocean Model, barotropic version, Center for Marine and Atmospheric Sciences. Institute of Oceanography, Hamburg, Germany) with the same temporal resolution.

Storm-induced topographic changes have been quantified by using LIDAR-derived topographies obtained before (16 October 2008) and after the storm impact (17 January 2009) by the Institut Cartogràfic i Geológic de Catalunya. The data were provided as high-resolution digital elevation models (DEMs) with a 1-m grid step, a maximum vertical precision of 2–3 cm and overall RSME of 6 cm (Ruiz et al. [36]).

The bathymetry of the study area has been obtained by combining an offshore grid with a spatial resolution of 0.28' derived from the GEBCO (General Bathymetric Chart of Oceans) database [37], and a finer inner topography while nearshore bathymetry was built by combining the LIDAR-derived topography and multi-beam bathymetric data provided by the Ministry of Agriculture, Fish, Food, and Environment, covering the whole area from the +2 to the −50 m with a 5 × 5 m resolution. Multiple bathymetries were available from different years, including 2006 and 2010, and these information was merged to properly fit the 2008 LIDAR shoreline and link the emerged topography with the submerged bed.

#### *2.3. The Sant Esteve 2008 Storm*

The storm of reference used in this study was a V-class (extreme) event, according to the classification of storms in the NW Mediterranean of Mendoza et al. [19]. This storm, known as the Sant Esteve storm, occurred on 26 December 2008 in the northern part of the Catalan Sea. It was created by the presence of a low-pressure center located over the Balearic Sea, with a minimum pressure of 1012 hPa, along with a high pressure center over northern Europe (1047 hPa). This is one of the typical mechanisms of cyclogenesis in the Mediterranean (e.g., Trigo et al. [38]), and it is the most common situation originating extreme storms along the Catalan coast (Mendoza et al. [19]). Under these conditions, the action of very strong NE winds in the Gulf of Lyon (wind velocities up to 20 m s−<sup>1</sup> were recorded at the coast) generated a wave field with a clear spatial pattern, with Hs values and power content decreasing from north to south along the Catalan coast (Jiménez [39]; Sánchez-Vidal et al. [30]). Thus, according to the data recorded by the Palamós buoy (see location in Figure 1), the storm lasted 73 h (Hs > 2 m) reaching a Hs of 7.5 m at the peak of the storm and a maximum wave height (Hmax) of 14.4 m. The associated return period of this storm was approximately 125 years according to the data on extreme climate obtained by Puertos del Estado [40] for this buoy and in light of data previous to the storm. The storm progressively lost strength as it moved south and, off the study site, the values recorded by the Tordera buoy showed a storm with a duration of 66 h, reaching a Hs at the peak of the storm of 4.65 m and Hmax of 8.0 m. Mean wave direction during the storm was E, which corresponds to the dominant direction during extreme storms in the area (Mendoza et al. [19]). Figure 2 shows the recorded values of wave parameters during the storm by the Tordera buoy (see location in Figure 1).

**Figure 2.** Wave records off the study site (Tordera buoy, see Figure 1) during the San Esteve storm. Blue dots indicate SWAN output at the same location.

The impact of the storm produced a significant coastal morphodynamic response in the form of erosion and overwash for most beaches along the northern part of the Catalan coast (Plana-Casado, [41]; Jiménez et al. [42]; Durán et al. [43]). Moreover, the magnitude of the storm was so considerable that many benthic ecosystems were also significantly affected (Sánchez-Vidal et al. [30]; Teixidó et al. [44]; Pagès et al. [45]).

Storm-induced topographic changes in the surroundings of the Tordera river mouth are shown in Figure 3. The observed response was different on both sides of the river, with the largest erosion taking place in the northern part, the s'Abanell beach. This beach is oriented to the East, nearly perpendicular to storm waves, and thus, the beach presented a generalized erosive behavior along its total extension (2.4 km). The volume of sediment eroded from the subaerial part of the beach was approximately 66,000 m3, with a beach-averaged erosion rate of approximately 30 m3/m and a maximum value of approximately 80 m3/m at its northernmost part (section SB-1 in Figure 3). Storm-induced wave overtopping occurred along the entire beach, and in its southernmost part, close to the river mouth, overwash deposits up to 6 m3/m were detected (section SB-3 in Figure 3). These volume changes resulted in a beach-averaged shoreline retreat of 11 m, with a maximum recession of approximately 30 m (Plana-Casado, [41]; Jiménez et al. [42]).

From the river mouth to the south, the coast experienced a different morphodynamic response. This section, the Malgrat de Mar beach, is oriented to the S, resulting in a large obliquity to E incoming waves during the storm. Just south of the river mouth a small post-storm accretion spot of approximately 7000 m<sup>3</sup> was detected. This seems to be related to the alongshore deposition of material eroded from the northern part. South of this area, the coastline shows a nearly generalized moderate erosion together with significant overwash deposits in the subaerial part of the beach (Figure 3). The spatially averaged erosion of the emerged beach was approximately 10 m3/m (one third of that observed for the northern beach) whereas the averaged overwash accumulation was approximately 7.5 m3/m (Plana-Casado, [41]; Jiménez et al. [42]).

**Figure 3.** Morphologic changes after the St Esteve 2008 storm. Dashed line pre-storm and continuous line post-storm profiles. Grey shaded area represents de position of the promenade (north) or road (south).

#### *2.4. Models*

Simulations have been done by using a model setup composed of the SWAN model (Simulating Waves Nearshore, version Cycle III v41.01, Delft University of Technology, Deltares, The Netherlands) (Booij et al. [46]; Ris et al. [47], TU Delft [48]), which propagates waves to the coast, and of the XBeach model (eXtreme Beach behavior, Kingsday version, Deltares, Delft, The Netherlands) (Roelvink et al. [10]), which is used to assess two storm-induced coastal hazards: inundation and erosion.

SWAN (Simulating Waves Nearshore) is a third generation wave model which is based on the wave action balance equation. It simulates short-crested wind generated waves by incorporating wave–water interactions, wind growth and dissipation processes such as whitecapping, bottom friction, and depth-induced breaking. For more detailed insight into these mechanisms controlling energy and wave propagation processes the reader is referred to the SWAN manual and to SWAN scientific technical documents at http://swanmodel.sourceforge.net/. The implemented model version uses the Komen et al. [49] formulation to calculate whitecapping. This version permits counteracting the previously reported under-predictions of significant wave height and wave period in areas characterized by fetch-limited conditions under the influence of transient winds (see e.g., Bolaños [50]; Pallares et al. [51]), which are typical conditions for the NW Mediterranean coast.

The model has been implemented using a nested grid configuration (Figure 1). A coarse grid with a total extension of approximately 80 × 70 km is used to transfer offshore wave conditions to the study area. The bathymetry grid has a spatial resolution of 0.28 whereas the wind field grid has a spatial resolution of 5 . This coarse setup is fed with wave spectra at 15 positions distributed along the offshore boundaries of the grid (Figure 1). The inner fine grid covers a domain of approximately 20 km × 26 km with a spatial resolution of 0.06'. This grid has been generated to properly reproduce the existing sharp changes in the bathymetry between intermediate-shallow waters due to the large steepness of the lower shoreface. This can permit a better simulation of wave propagation in the region close to the XBeach coastal grid (external boundary at 20 m water depth). A limitation of the SWAN model is its inability to simulate storm surges. To characterize storm surges in the study area we use water level predictions obtained with the HAMSOM model implemented by Puertos del Estado (Ratsimandresy, et al. [52]) at three locations close to the study site. The storm surge contribution to the total water level in the study area is of low magnitude and significantly smaller than wave-induced runup during storm conditions (e.g., Mendoza and Jiménez [34]). The validated SWAN model has been used to propagate waves from deep to shallow waters during the entire storm duration. Since a simultaneous time series of measured wave data was available within the modelled domain, it was possible to obtain transfer coefficients for wave conditions (wave height and direction) for any point in the grid with respect to buoy location. With this information, the recorded wave conditions at the buoy have been transferred to selected grid points to be used as input data for Xbeach modelling.

Once the forcing conditions during the storm have been propagated from deep water to nearshore, the remaining task is to propagate these conditions to the coast and to assess the magnitude of storm-induced hazards, i.e., erosion, overwash and overtopping, which is done with the XBeach model. XBeach is an open-source 2D depth averaged model which solves wave propagation, flow, sediment transport and bed level changes for varying wave and flow boundary conditions (Roelvink et al. [10]). It solves the time-dependent short wave action balance on the scale of wave groups, which allows for the reproduction of directionally spread infragravity motions (so called surf-beat) along with time-varying currents. The frequency domain is represented by a single representative peak frequency, assuming a narrow banded incident spectrum. Shallow water momentum and mass balance equations are solved to compute surface elevation and flow. Additionally, to solve the contribution of short waves to mass fluxes and return flows, XBeach uses the Generalized Lagrangian Mean formulation (Roelvink et al. [10]).

Sediment transport rates are calculated from the spatial variations in depth-averaged concentration, which are calculated from advection-diffusion equations with a source-sink term based on an equilibrium sediment concentration. The equilibrium concentration takes into account both the contribution of the suspended and bed loads by means of the Soulsby-Van Rijn formulation (Soulsby [53]) with a limitation of the maximum stirring velocity based on the Shields number at the start of the sheet flow (McCall et al. [11]). For deeper insight into the XBeach description, setup and equations, see Roelvink et al. [10].

The model has been implemented by using a curvilinear grid with variable cell size in both alongshore and cross-shore directions (see Figure 4). The extension of the mesh is approximately 1.5 km in the cross-shore direction, with cell size ranging from 5~6 m at the offshore boundary (20 m depth) to 0.7–0.8 m at the swash zone. In the alongshore direction the model has an extension of 4.5 km and the cell size ranges from 25 m at the lateral boundaries down to 2–3 m around the river mouth. The grid was obtained by means of the Delft3D-RGFGRID module, imposing consecutive maximum cell-size changes ~5–10% to ensure proper smoothing, while maintaining valid orthogonalisation. The final result has 669 alongshore by 568 cross-shore nodes. The model wave boundary conditions consist of wave characteristics specified at four different locations along the offshore boundary (see Figure 1), with a time-step of 1 h, which is the same resolution of the original data used to force SWAN. The use of four different locations aims to capture the difference in wave conditions on both sides of the river due to the different coastline orientation. Water level variations during storms to be simulated are directly introduced in the model from HAMSOM simulations, also with a time-step of 1 h. XBeach model computations (i.e., hydrodynamics and morphodynamcis) are performed with a temporal resolution of 1 s (dtbc = 1).

**Figure 4.** Distribution of cross-shore and alongshore grid resolution over the XBeach domain. Orange boxes show the location of post-processing subdomains, being SN and SS sectors (**left**) for the morpholodynamic analysis and NORTH, SOUTH sectors (**right**) for the inundation assessment.

The morphology of the study area, characterized by coarse sand and steep reflective beaches, makes XBeach modelling a demanding exercise in terms of predicting beach morphodynamic response during storms (see e.g., Vousdoukas et al. [17,54]). Notably, coarse sediment environments are characterized by a lower frequency of the avalanching processes, a greater importance of the bed load over the suspended transport, and higher importance of mechanisms such wave asymmetry, or water infiltrations and groundwater effects. All these are by default configured in XBeach to work for fine sand environments, and must be revised and modified for its application at the study site, modelled with a D50 of 1.3 mm and a D90 of 1.9 mm.

To properly reproduce morphodynamic changes, both the surf-beat and the non-hydrostatic modes of the XBeach model were initially tested. The surf-beat model was observed to under-predict overwash when using typical XBeach-grid resolution for straight and mild-slope coasts. This under-prediction is also caused by the lower contribution of infragravity waves to the total run-up in steep profiles (Wright and Short [55]). However, the surf-beat model accurately reproduced the alongshore patterns of erosion in the entire domain. In contrast, the non-hydrostatic model was observed to have better performance in reproducing wave-by-wave run-up but lower accuracy in reproducing the alongshore morphodynamic patterns of erosion and deposition, with a quite higher computational cost (since it requires denser grids). An additional difficulty was using XBeach with a curvilinear grid to properly reproduce alongshore processes driven by the change of coast orientation at both sites of the river mouth. The final adopted approach was to use the surf-beat model with a higher resolution than the typically used in straight and gentle-slope coasts, which considerably improved the overwash

prediction. The average cross-shore resolution of 5.2 m (typically 20~25 m) at the offshore boundary, going down to 0.7 at the swash zone, ensured a better reproduction of wave propagation. This setting aimed to properly capture the abrupt changes in the bathymetry from the 20 m to the 5 m depth by using a larger number of cells. In addition, the average alongshore resolution around the river mouth, where the largest alongshore gradients are present, was increased up to 2.3 m (typically 5~10 m) (Figure 4).

#### *2.5. Scenario Testing*

Since the recorded Sant Esteve 2008 event had the characteristics of a V-class extreme storm in the area, with the largest recorded Hs and with a typical E direction, it was used as base case scenario (C0). From this, six additional scenarios were defined to test the effect of different incoming wave directions. The procedure consisted of simulating the exact same wave time series as the Sant Esteve 2008 storm (keeping the same Hs, Tp, and directional spreading) and only changing the mean wave direction. This approach permits maintaining the same storm wave intensity as the reference case at the offshore boundary but under a different direction. The tested conditions were two scenarios where wave direction was shifted 20◦ and 40◦ counter-clockwise to the North (C20− and C40−) and four scenarios shifting 20◦, 40◦, 60◦, and 80◦ clockwise to the South (C20+, C40+, C60+, and C80+). Thus, seven different scenarios, including the baseline condition, were simulated and compared to assess the differences in storm-induced hazards under incoming directions ranging from ~60◦ N (C40−) to ~180◦ N (C80+). The scenarios have been chosen to cover the tipical range of incoming conditions at the −20 m depth with directional spans of 20◦.

The magnitude of storm–induced hazards was quantified in different control sectors along the study area to capture main factors potentially affecting the beach response to different incoming directions. Thus, the morphodynamic response was analyzed in five sectors: two 250 m long sectors northward of the river mouth (SN1, SN1), and three sectors southward of the river mouth (Figure 4). These S sectors are as follows: (i) SS1, a 200-m-long stretch between the river mouth and an existing rigid structure at the shoreline; (ii) SS2, a 200-m-long sector, southward of the mentioned existing structure; and (iii) SS3, a 500-m-long stretch at the southernmost end. At each sector, three variables are used to characterize morphodynamic changes: erosion volume (m3/m), overwash volume (m3/m) and profile retreat (m). All of the variables are calculated within sectors from XBeach gridded output, from which sector-averaged values and standard deviations are derived. Erosion volumes are computed in the inner part of the beach, from the subaerial part down to the −2 m level, which roughly determines the water depth where main inner profile changes occur. Overwash volumes are computed as deposited sediment volumes in those parts of the subaerial beach where vertical growth is detected. The profile retreat is measured at three elevations at the beachface (1, 1.5, and 1.75 m above mean water level).

To characterize inundation, just two sectors were selected, one for the region at the north of the river mouth and another to the south (Figure 4). The selected variable to characterize inundation was the inundation surface (Ha) over different thresholds of inundation depth: 0.05, 0.25, 0.5, and 1 m.

#### **3. Results**

#### *3.1. Base Case Scenario (C0). The Sant Esteve Storm*

This base case scenario corresponds to the model validation using the recorded Sant Esteve 2008 Storm event. The SWAN model was validated with wave conditions recorded by the Tordera wave buoy and by comparing the measured and modelled wave conditions (Hs, θ), shown in Figure 2. As can be seen, the model reproduces well wave parameters during the storm, especially during the peak of the event, when simulated variables almost coincide with recorded ones. The obtained root mean squared error (RMSE) of the model during the entire duration of the storm is 0.53 m in Hs and 12.5 degrees in θ, with the largest contributions to the error taking place during the relaxation phase of the storm, whereas error values at the storm peak are significantly lower (Figure 2). Since the largest

storm-induced morphodynamic changes occurred during the storm peak, we can accept that the use of the SWAN model to simulate wave propagation in the study area is acceptable for the purposes of this research.

To calibrate XBeach, the model's results were compared with LIDAR measurements of the emerged beach. The calibration of the model was performed by adopting a double approach: (i) by optimizing the Brier Skill Score (BSS), which quantifies model performance by comparing model output to the real post-storm LIDAR measurements of the emerged profile; (ii) by performing a qualitative assessment of the modelled features, such as alongshore and cross-shore patterns of bed level changes, magnitude, and location of the overwash deposits and validation of the inundation reach according to the available qualitative information on the event. The used BSS to characterize model predictive skill takes into account the measurement error (ΔZe) as in Harley and Ciavola [56], and thus, the BSS score is given by:

$$\text{RSS} = 1 - \left(\sum \left( |\mathbf{z\_{mf}} - \mathbf{z\_{mod}}| - \Delta \mathbf{Z\_{e}} \right)^{2} / \left(\sum \left( \mathbf{z\_{mf}} - \mathbf{z\_{mi}} \right)^{2} \tag{1} \right)$$

where zmf is the final LIDAR measured bed level, zmod the final modelled bed level, and zmi the initial bed level. Here, ΔZe is considered as the LIDAR measurement error (i.e., RSME of the overall LIDAR product), which is 0.06 m. According to Sutherland et al. [57], the classification of models' performance based on the BSS score can be considered to be very good for values over 0.4 and excellent for values over 0.5–0.6. Due to the morphology of the study area, which induces a differentiated morphodynamic response along the coast (Figures 5 and 6), three BSS were calculated: (i) a global BSS, which is calculated for the entire area; (ii) a local BSS at the north of the river mouth; and (iii) a local BSS southward of the river mouth. Since the BSS assessment can only be performed for the emerged profile (where pre- and post-storm topographic data exist), a qualitative assessment of the modelled submerged profile was also performed. To this end, we analyzed the final shape of the modelled submerged profile taking into account the expected typical morphodynamic response under storm conditions at the site.

**Figure 5.** XBeach validation with LIDAR measurements of the emerged morphological changes. BSS: Brier Skill Score.

**Figure 6.** Comparison between XBeach simulated and measured beach profiles after the impact of the storm (see profiles location in Figure 5).

Although different combinations of model parameters resulted in BSS scores over 0.4 for the emerged profile, the qualitative assessment highlighted that some of them produced excessive deposition volumes in the submerged part. The final setup parameters, which resulted in an overall BSS of 0.68 (measured at the northern beach and first 600 m of the southern beach) and a meaningful qualitative simulation of the predicted submerged profile and alongshore bed level change patterns, are shown in Table 1 along with parameter description and tested ranges.


**Table 1.** Calibration parameters. Range of tested values and final parameter setup.

<sup>1</sup> based on best BSS score and qualitative assessment.

The implications of the selected values are as follows. First, the increase in wave-attack on the coast improves the behavior of the model in terms of reproducing observed overwash deposits (gamma and delta). Second, the contribution of avalanching to bed level changes was limited, taking into account the steepness of the site (wetslp). Third, the non-linearity effect on sediment transport for steep profiles (facAs and facSk) may be taken into account as reported in Elsayed and Oumeraci [18]. Fourth, sediment particle mobilization may be limited (sedcal). With respect to this, other authors

have already reported an excess of modelled sediment suspension because the shear stress values required to initiate particle motion are higher than those predicted by using the Shields curve, as noted in Elsayed and Oumeraci [18] or McCall [11]. Fifth and finally, the groundwater module was turned on (gwflow) because the role of infiltration is more significant in coarse sediment environments, such as the Tordera Delta, where grains are close to gravel-size. The value of the permeability factor has been estimated using the Kozeny–Carman formula as described in Carrier [58]. For the grain size in the study area, the value of the permeability factor was estimated to be 0.0058 m/s. The results of the storm-induced morphological changes simulated with the final adjusted set of model coefficients are shown in Figures 5 and 6. As it can be seen, the measured changes in beach elevation at both sites of the river mouth are well reproduced by the model, with a BSS of 0.68 when the entire area is considered. The model also properly reproduces the differentiated response at both sites of the river mouth, mimicking the effects of change in coastal orientation with respect to the storm wave direction and the differences in coastal morphology.

Figure 7 shows the modelled depth-averaged wave-induced circulation during the peak of the storm, where a different circulation pattern is observed at both sides of the river mouth. Mean XBeach output (average conditions at 1 h output time-step) is used to average the current velocity and sediment transport during 4 h around the maximum peak of the event. At the northern part, where the beach is fully exposed to storm waves (i.e., a relatively small obliquity during the storm) and without any submerged morphological features, wave-induced circulation shows a typical quasi-uniform longshore current structure along the beach until the river mouth at the southern end. At the southern part, the coast is partially sheltered from storm waves (i.e., a large wave obliquity during the storm) and there is an alongshore bar running parallel to the shoreline with varying crest levels. This bar delimits a shallow shelf of approximately 4 m water depth to the shoreline. In this area, the induced circulation pattern during the peak of the storm shows a longshore current directed towards the south close to the shoreline and a local inversion of the current towards the N over the shallow shelf (Figure 7).

**Figure 7.** Simulated depth-averaged currents and sediment transport. Mean XBeach output at each 1 h time-step is averaged during the 4 h of maximum storm intensity.

With respect to morphodynamic changes (Figures 5 and 6), the model predicts, for the northern part, a nearly continuous erosion along the beach without significant alongshore gradients in sediment transport, and with significant sediment mobility down to about −5 m depth. The model properly reproduced the observed increase in erosion magnitude from B2 (P4) to B1 (P3), as well as most of the overwash deposits, which increase towards the south as the height of the berm decreases. In summary, the model reproduced well the observed variability in beach erosion and overwash along this northern site with a local BSS of 0.75. Southwards of the river mouth (Malgrat Beach), induced sediment transport and beach erosion present significant alongshore variations, with the largest erosion taking place just southwards of an existing structure located at the northern part of the sector, A2 (P2), which should act as a local boundary condition. This local effect can be seen in Figure 7 where a gradient in the longshore sediment transport, starting at the rigid structure north of A2, is detected. Longshore transport rates are larger than in the north, and mainly concentrated down to −3 m water depth. The model reproduced observed beach topographic changes with large erosion and overwash deposits due to a relatively low beach berm (A2, P2). At the southernmost part, erosion is only taking place at the upper level of the profile, whereas part of the material is deposited around zero level (A1, P1). The obtained local BSS for this sector was 0.60.

#### *3.2. The E*ff*ects of Wave Direction on Storm-Induced Hazards*

As was previously mentioned, once the morphodynamic model was calibrated, it was used to analyze the sensitivity of the area to changes in wave direction during storm impacts. Simulated morphodynamic changes and inundation for cases C20− (~80◦ N) to C40+ (~140◦ N) at northern and southern parts of the study area are shown in Figure 8, whereas the variation of integrated control variables for each sector can be seen in Figures 9 and 10 for inundation and morphodynamic changes respectively.

**Figure 8.** XBeach-simulated inundation depth (m) (**bottom**) and bed level changes (m) northwards (**top**) and southwards (**middle**) of the Tordera River.

With respect to inundation, the beach northward of the river mouth (Figure 9a) experiences an increase of the predicted inundation surface from C40− to C20+, when the surface reaches its maximum extension (5.47 Ha with 1 Ha over 0.5 m depth). Wave direction in this scenario corresponds to a nearly normal wave attack on the local coastline orientation. As incoming wave direction continues shifting towards the south, the predicted inundation surface progressively decreases, reaching values for the C80+ case similar to those observed for the C40− case. At the southern coast, the obtained pattern is significantly different (Figure 9b). Thus, for wave direction scenarios dominated by NE components (C40−, C20−) no significant inundation is observed, and this is consistent with a large sheltering from highly oblique wave incidence. For wave directions shifting from the base scenario to the south (C0 to C80+), the predicted inundated surface increases, reaching a maximum value for scenarios C60+ and C80+ of approximately 74 Ha (for an inundation depth above 0.05 m) and 33 Ha (for an inundation depth of 0.5 m), the C80+ case. Taking into account the total inundation of the study area, we can state that the magnitude of the inundation hazard significantly increases as wave direction shifts to the south. The southern coast is the most affected in terms of the magnitude of expected changes due to the local low-lying topography.

**Figure 9.** Variation of simulated inundated surface for different inundation depths as a function of the simulated storm direction for the northern (**a**) and southern (**b**) control areas (see Figure 4 or Figure 8 for location).

Regarding morphological changes, the northern beach (S'Abanell) shows a relatively low sensitivity to wave direction for scenarios C40− and C20+ when erosion volumes and profile retreats are considered (Figure 10). For this range of wave directions, the local wave-induced circulation is characterized by a southward directed longshore current along the beach (see Figures 7 and 11) which turns north when the wave direction is C20+ (Figure 11). As the incoming wave direction turns southwards, a stronger north-directed alongshore current is induced with velocity increasing along the beach (Figure 11) which results in an increasing erosion and profile retreat. Figure 11 also shows hatching areas in the velocity field which originate in those areas where wave incidence is orthogonal to the nearshore bathymetry, characterized by a heterogeneous bar in front of the coast. This beach erosion increase is particularly observed at the southernmost end of the section (SN1), just northwards of the river mouth, where an existing revetment acts as a boundary condition (barrier) for northwards directed transport. This overall modelled behavior is consistent with local field observations, where the northern part of the beach (out of the domain) experiences a significant sediment deposition under the impact of southern storms, bringing sediment from the south.

**Figure 10.** Variation of simulated morphodynamic parameters (see text for description) at selected control areas the N (**left**) and S (**right**) coasts for tested storm directions (see Figure 4 or Figure 8 for location). Continuous lines denote variable mean, shaded areas represent standard deviation and dashed lines indicate maximum profile retreat at each sector.

**Figure 11.** XBeach simulated depth-averaged currents. Mean XBeach output at each 1 h time-step is averaged during the 4 h of maximum storm intensity.

Overwash deposits along the northern section present a variation pattern with wave direction consistent with modelled inundation. The largest overwash verifies in the southern end of this sector (SN1) due to its lower beach berm. Thus, maximum overwash deposition verifies at SN1 under C20+ and C40+, which were the scenarios producing the largest inundations. As wave obliquity increases (scenarios C60+ and C80+), overwash significantly decreases, which is also in agreement with the observed inundation decrease under these conditions.

Southwards of the river mouth (Malgrat de Mar), the analysed sectors show a differentiated response due to the local effect of the existing building in the coastline which acts as a boundary condition for sediment transport. The northernmost end, SS1, just southward of the river mouth, experiences accretion (i.e., positive profile displacement) under C40− and C20− storms, which would correspond to the deposition of sediments brought by the southward directed longshore transport from the northern beach. As wave direction shifts to the south (from C20+ to C80+) this area becomes slightly erosive. This morphodynamic behavior is consistent with the observed switch in alongshore current direction for scenarios C20+ and C40+ (Figure 11). The central sector SS2, located just southwards of the mentioned obstacle, shows a different morphological response with storm direction to that observed in SS1. Thus, the largest erosion rate verifies under C0 due to combined effect of the orientation of the coast, a gap in the submerged bar and the presence of the hard structure at the northern end of the sector. This can be seen in Figure 7, where a significantly large gradient in sediment transport rates under C0 is observed. As wave direction shifts towards south, erosion rates significantly decrease, reaching a nearly null integrated value under C40+ to C80+ scenarios. The largest erosion under scenarios C20− to C0 in this sector, in comparison with the rest of the southern coast, seems to reflect the effect of the revetment on enhancing local downcoast erosion, similar to the well-known flanking effects in seawalls (e.g., Kraus and McDougal [59]). The southernmost sector, SS3, is far from the local effect of the revetment, so much so that the aforementioned downcoast erosion enhancement is not detected and, under highly oblique storms, erosion rates are significantly lower. Total volumetric changes in this sector present a relative low sensitivity to changes in wave direction for southern wave scenarios (C20+ to C80+).

Regarding overwash, the southernmost section SS3 presents a wave-direction influence similar to the influence observed for inundation, with overwash rates growing as wave direction turns to the South. The sectors closer to the river mouth SS2 and SS1 present a similar response with a maximum for C20+ and C40+ scenarios, when the wave directions are close to local orthogonality (Figure 8), and they are significantly larger for SS2 (~12 m3/m) than SS1 (~2 m3/m).

#### **4. Discussion and Conclusions**

In this study, the potential effects of changing wave direction for the storm-induced hazards on a highly curvilinear coarse sandy coastline have been assessed. This sensitivity test has been selected because although storminess projections under climate change scenarios for the Western Mediterranean do not predict any increase in wave height (e.g., Lionello et al. [20]; Conte and Lionello [21]), some existing projections identify potential changes in wave direction (Casas-Prat and Sierra [22,23]). These changes in wave direction may have significant implications for coastal sediment transport and coastal stability, as has been confirmed for the interannual changes influenced by El Niño (e.g., Barnard et al. [60]). Moreover, regarding cuspate coastlines such as the study area, their greater sensitivity, due to their curvature, results in even more significant implications (e.g., Slott et al. [24]; Johnson et al. [25]).

The tested hypothesis is that changes in wave direction may cause large variations in the magnitude of storm-induced hazards. This effect has also been addressed in other studies such as those by Mortlock et al. [26] and de Winter and Ruessink [27], which specifically analyzed the effects of changes in wave direction on the storm-induced hazards in the SE Australian and Holland coasts respectively. To this end and to isolate the influence of wave direction, we used a recorded long-return period storm as a base case scenario and we built test scenarios just by changing wave direction while maintaining the other wave parameters as recorded during the base storm (wave height and period).

In any case, tested conditions have not been designed to be used as climate change induced projections, as that may require the proper forecasting of regional wave conditions under given climate scenarios (e.g., Casas Prat and Sierra [23]). These have to be considered from the perspective of coastal risk management, in which a set of possible conditions are analyzed to characterize coastal vulnerability and resilience to inform risk management under uncertainty (see e.g., Hinkel et al. [61] for an application of this perspective to analyze sea level rise). In the study area, the current storm wave conditions depend on direction, with largest wave height and power being associated with NE-E waves, whereas S storms are less frequent and present a smaller associated power (e.g., Sánchez-Arcilla et al. [62]; Mendoza et al. [19]). To assess the potential variability on storm-induced hazards, tested scenarios were built by just changing wave direction while the remaining recorded parameters (representative of a worst case scenario, according to recorded conditions) were maintained.

This analysis has been performed by using the SWAN and XBeach models to simulate storm-induced hazards. Both models were calibrated by using data recorded during the impact of an extreme storm recorded in December 2008, which is used as the base case scenario. Although it is desirable to use more than one event to properly calibrate/validate the models (e.g., Ranashinge [63]), data availability during storm conditions was restricted to this event. However, on the positive side, it has to be considered that this storm was the largest event recorded in the area and representative of extreme storms with a very long return period (Mendoza et al. [19]) under current climate conditions. The SWAN model was very successful in simulating wave conditions during the development phase of the storm up to the pass of the peak of the storm, with the larger differences between measured and simulated waves being detected during the relaxation phase of the storm, when most of the induced changes had already occurred. The default parametrization of the XBeach model had to be adapted for application at the side to represent the effects the coarse-sand environment. Sediment transport was limited by using the sedcal parameter, avalanching was limited by increasing the critical slope, wave asymmetry was increased as suggested in literature for steep slopes (Elsayed and Oumeraci [18]) and groundwater effects were included. Gamma and delta wave breaking parameters were also tuned (Table 1). Calibrated parameters setup for XBeach in the study area led to a BSS score of 0.68 in spite of the out-of-comfort tested conditions (i.e., highly curvilinear coast, steep beach, coarse sediment). Although the predictive skill was very good for the northern and southern beaches, the model performance was better in the northern domain (BSS = 0.75) than in the southern one (BSS = 0.60), since this last area presented a significantly larger obliquity to wave direction during the storm, and a more complex bathymetry.

The obtained results show a very high sensitivity of storm-induced processes, i.e., inundation and erosion, to changes in storm wave direction. With respect to inundation, expected changes in hazard magnitude are very significant, especially in the southern part of the study area, since its morphology is characterized by a lower berm, and its low-lying unprotected hinterland makes this area sensitive to storm flooding (Jiménez et al. [7]). Thus, as storm waves turn from the base case (C0) to the south, the inundated surface along this southern beach dramatically increases due to its direct exposure to the south. On the contrary, a potential shift of wave direction to the N will have a positive impact on inundation in this area, since it will be more sheltered from wave action. At the northern beach, the largest increase in inundation hazard verifies under C20+ and C40+ scenarios when waves face nearly orthogonally to the coastline, although due to local morphology, the affected surface is much lower than in the southern beach. The hinterland of the study area is mostly occupied by agriculture land and, in the outer fringe just behind the shoreline, by campsites. In this sense, to transfer the potential change in hazard magnitude to changes in damage risk, it should be important to consider not only the change in direction but also its seasonality. Thus, risk may vary dramatically between the summer season (when the campsite facilities are used by visitors) and the rest of the year when only installations will be affected (e.g., Merz et al. [64]). An analysis of the risk associated with storm-induced inundation for different storm conditions can be seen in Sanuy et al. [32].

Similarly, storm-induced morphodynamic changes are more sensitive to directional changes on the southern beach, where the magnitude of the changes is larger. The beaches at the south of the river mouth present a larger spatial variability than those in the north due to the presence of a local boundary condition in form of a revetment at the shoreline. This revetment, which modifies local longshore transport, significantly enhances downcoast erosion under storm conditions. This induces a southwards directed longshore sediment transport while simultaneously promoting the accumulation of upcoast sediment. This contrasting behavior is particularly observed in the base case scenario which seems to represent the optimum conditions for longshore sediment transport in the area, thus inducing the largest changes in the surroundings of the structure.

In a particular case, under C40− and C20− scenarios, when wave direction turns north, the beach sector just south of the river experiences an important sediment accumulation due to the apparently efficient transfer of sediment from the northern beach across the river mouth and the partial barrier effect of the existing revetment.

The magnitude of the erosional response along the two control sectors in the northern beach is similar, although a higher variability is detected in the area closest to the river mouth. In general, there is a slight increase in erosion rates as wave direction turns south. This variation should be indicative of the role of longshore sediment fluxes during storm conditions. Thus, as the controlled northern area is just besides the river mouth, where there is another structure acting as a boundary condition, the increase in longshore sediment transport as waves turn S (scenarios from C20+ to C80+) will increase sediment losses, which will be transported further to the north. This behavior is currently observed in the northernmost part of this beach (out of the control zone in Figure 10) which experiences sediment accumulation under the impact of southern storms.

As expected, changes in the magnitude of overwash deposits follow observed changes in inundation, i.e., they increase as wave direction turns to the south, with maximum values around C20+ and C40+. The exception to this is the predicted changes in the southernmost sector, which present the largest overwash for C80+ conditions. The spatial variability in the northern beach is significantly lower than in the south, with small variations in magnitude across the tested range. Moreover, and reflecting the observed differences in inundation, the magnitude of overwash deposits is much higher in the southern sector.

Finally, and as a concluding remark, this analysis has shown that storm-induced hazards along a highly curvilinear coast are extremely sensitive to changes in wave direction. This means that even under a climate scenario of relatively steady storminess (wave power and frequency), a potential shift in wave direction may significantly change hazard conditions and, in consequence, need to be accounted for in robust damage risk assessments. To this end, an analysis such as the one presented here also permits an assessment of how coastal geomorphology modulates induced changes. In the study area, the low-lying nature of the southern beach and its orientation with respect to the current dominant storm direction make this area much more sensitive to directional changes. This is especially relevant from the coastal management standpoint because this area has been already identified as a hotspot for storm impacts under current conditions. The use of detailed process-based models has permitted the identification and quantification of the drastic increase in sensitivity when anthropogenic perturbations are present along the coast. These perturbations act as boundary conditions modifying local hydrodynamics and associated transport. For the case study analyzed here, the obtained results clearly identify the hazardous potential of the existing revetment in the southern beach, which has also been identified under current conditions, suggesting that its removal will soften the estimated morphodynamic response.

**Author Contributions:** Conceptualization, J.A.J. and M.S.; methodology, J.A.J. and M.S.; software, M.S.; validation, M.S.; formal analysis, M.S.; investigation, M.S. and J.A.J.; data curation, M.S.; writing—original draft preparation, M.S. and J.A.J.; writing—review and editing, M.S. and J.A.J.; supervision, J.A.J.; project administration, J.A.J.; funding acquisition, J.A.J.

**Funding:** This study was conducted in the framework of the RISC-KIT (Grant No 603458) and the M-CostAdapt (CTM2017-83655-C2-1-R) research projects, funded by the EU and the Spanish Ministry of Economy and Competitiveness (MINECO/AEI/FEDER, UE) respectively. The first author was supported by a PhD grant from the Spanish Ministry of Education, Culture and Sport.

**Acknowledgments:** Lidar data were obtained thanks a collaboration agreement between the Institut Cartogràfic i Geològic de Catalunya (ICGC) and the LIM/UPC, with the aim of analysing the usefulness of LiDAR data for monitoring coastal processes. The authors also express their gratitude to the Generalitat de Catalunya and Puertos del Estado for supplying wave data and to the Ministry of Agriculture, Fish, Food and Environment for the bathymetric data used in this study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Impacts of Sea Level Rise and River Discharge on the Hydrodynamics Characteristics of Jakarta Bay (Indonesia)**

#### **Martin Yahya Surya 1, Zhiguo He 1,2, Yuezhang Xia 1,2,\* and Li Li 1,2**


Received: 8 May 2019; Accepted: 3 July 2019; Published: 5 July 2019

**Abstract:** Jakarta city has been vulnerable to sea level rise and flooding for many years. A Giant Seawall (GSW) was proposed in Jakarta Bay to protect the city. The impacts of sea level rise and river discharge on the tidal dynamics in Jakarta Bay and flooding areas in Jakarta city were investigated using the finite-volume coastal ocean model (FVCOM). Model results showed that the bay is diurnally dominated by the K1 tidal component. The diurnal tides propagate westward, while the semidiurnal tides propagate eastward in the bay. The rise of sea level increases the diurnal tidal component and the inundation areas due to the increased tidal forcing: when considering a sea level rise of 0.6 m, the K1 amplitude increases by ~1% (0.25 cm) near the coastline and the current magnitude increases by 16.6% (0.05 m/s). The inundation area increases with the sea level rise in the low land elevation areas occurring near the coastlines: the inundation area increased by 29.68 km<sup>2</sup> (7.1%) with a sea level rise of 0.6 m. The increase of river discharge amplified the diurnal tidal component as well as the inundation areas at the river mouth due to increased fluvial forcing: if 10 times the mean river discharge occurs, the K1 amplitude increases by ~1% (0.25 cm) and the current magnitude increases by 100% (0.4 m/s), and the inundation areas increase by 26.61 km<sup>2</sup> (6.2%). The K1 tidal phase remains almost unchanged under both the sea level rise and river discharge conditions. The combined increase of sea level rise and the river discharge amplifies the inundation areas and the tidal currents due to increased tidal and fluvial forcing. The construction of GSW would decrease the tidal prism and dissipation effects of the bay, thus slightly increasing the K1 amplitude of the tidal level: by less than 1% (0.2 cm). There would be no significant change of phase lag for the K1 component. Although this study is site specific, the findings could be applied more widely to any open-type bays.

**Keywords:** Jakarta Bay; sea level rise; river discharge; flooding; tides; currents

#### **1. Introduction**

Jakarta Bay is located on the northwest coast of Java Island, extending from 5.8◦ S to 6.2◦ S and 106.6◦ E to 107.1◦ E [1]. It is a shallow bay with an average depth of about 18 m and an area of 662 km<sup>2</sup> [2]. The bay is bordered by Pasir Cape to the west and Karawang Cape to the east, with a total coastline length of about 72 km.

Jakarta Bay is essential for the development of Jakarta city, the capital city of Indonesia. Jakarta city is the economic, cultural and political center of Indonesia. There are 10 million people in Jakarta city and more than 28 million people in the satellite cities of Jakarta [3]. The large population and the rapid development of the city have brought increasing urban development issues. The demographic change influenced by the urbanization causes large-scale rural land conversion and new town development in Jakarta city. This massive urbanization has led to some environmental issues. One of the environmental

problems is related to flooding and sea level rise. Flood, in terms of river flooding and coastal flooding, is the most common natural disaster affecting populations in Asia [4]. For example, extreme rainfall events will increase in frequency and severity within the twenty-first century and lead to extreme river discharges [5].

Jakarta city experiences annual flooding disasters due to heavy rains, high river discharges, and high tides. An extreme flood occurred in 2007, when north Jakarta was hit by unusually high tides and heavy rains. The total inundation area in the extreme flood covered 400 km2. Approximately 70% of Jakarta was flooded with water depth of up to 4 m. This extreme flood cost was around 890 million USD in direct losses and more than 620 million USD in indirect losses [6]. Nearly 600,000 people were displaced and 79 people died. Many roads were closed, including the highway to the international airport, and electricity and telephone lines were cut. Some researchers believed that land subsidence contributed to the heavy flood in 2007. The land subsidence due to groundwater extraction was the main cause for lowering the land topography. At that time, the rate of land subsidence was about 1–10 cm per year [7].

As well as the flooding issue, the problem of rising sea levels has become a major concern in Jakarta Bay. During the 20th century, the mass loss of glaciers and ice caps contributed to sea level rise [5]. These losses were estimated to be 0.5 ± 0.18 mm in sea level equivalent (SLE) per year between 1961 and 2003. Global warming of ocean waters produces water molecule expansions that are responsible for sea level rise. The average contribution of the thermal expansion to sea level rise is 0.42 ± 0.12 mm.

To overcome the flooding and sea level rise, a Giant Seawall (GSW) at Jakarta Bay has been proposed to protect the city from flood events increasing due to subsidence, exceptional high tides, and sea level rise [8]. The GSW project will include sea wall construction, land reclamation, toll roads, and a port extension in Jakarta Bay. The GSW will offer a long-term protection and flooding solution, and at the same time facilitate social-economic development. The GSW project may also have large impacts on tidal level, flow pattern, current and flooding area in Jakarta Bay. There have been many studies around Jakarta Bay, but relatively few focused on river discharge and subsequent sea level rise effect. Most of the applied models of Jakarta Bay have taken into account only tides and winds as the predominant forces to study circulation and hydrodynamics.

Ningsih's model for Java Sea showed that monsoon wind plays important roles in the general circulation of this area [9]. The flow pattern in Jakarta Bay is mainly driven by monsoon winds with current moving easterly during northwest monsoons (NWM), and reversed during southeast monsoons (SEM) [10–12]. In general, the southeast monsoons (SEM) are stronger than northwest monsoons (NWM) [13]. During flood conditions, the westward current induces a clockwise circulation through Jakarta Bay, and an eastward counter-clockwise circulation at ebb conditions [1,14].

Land subsidence in Jakarta is the predominant driver exacerbating coastal floods, followed by the influence of sea level rise [4]. The potential flood area in the 2025–2050 period is estimated to be 3.4 times greater than during the 2000–2025 period [15]. The sea level rise effects on tidal hydrodynamics along the northern Gulf of Mexico have been observed and modelled [16]. Tidal amplitudes within the bays increased by as much as 67% (10 cm) in some areas. Most of the bays would experience faster tidal propagation due to changes in harmonic constituent phases in future scenarios. The tidal inundation increases along the NGOM study area, especially in low-lying marsh areas and barrier islands with low dune elevations.

In recent years, estuaries and coastal system have been modelled successfully using the finite-volume coastal ocean model (FVCOM). Fukuoka Bay has been modelled using the FVCOM and Princeton Ocean Model (POM) [17]. A comparison of these two models showed that FVCOM provides a better result because FVCOM can describe the coastline geometry better. Benoa Bay, which has shallow water depth, has also been modelled with FVCOM [18]. FVCOM can successfully produce the residual current and vertical salinity profile, which is significantly affected by fresh water discharge. In Saemangeum in Korea, FVCOM has been used to successfully explain the mechanism of decreasing tidal current and tidal amplitude after the construction of a sea dike [19]. Thus, FVCOM is a powerful

tool for modelling the hydrodynamics in Jakarta Bay. The unstructured grid of FVCOM would be suitable for the irregular coastline and bathymetry of Jakarta Bay.

In this paper, an FVCOM model for Jakarta Bay was used to construct tidal dynamics and the flooding occurrence in Jakarta city. The impacts of the river discharges and sea level rise on tidal dynamics and flooding areas were numerically examined, and the effects of the GSW on tides and flooding were investigated.

#### **2. Materials and Methods**

The influence of river discharge and sea level rise on the hydrodynamics of the bay was investigated through numerical modelling. The hydrodynamic model (Figure 1) in this paper was based on the model previously established by Rusdiansyah et al. [1], with the model domain being expanded to cover the entire area as described in the Master Plan National Capital Integrated Coastal Development [8]. The impacts of river discharges and sea level rises were analysed using numerical experiments in this study.

#### *2.1. Hydrodynamics of Jakarta Bay*

Jakarta Bay is an open type bay (Figure 1) with a water depth less than 60 m. Surface currents in Jakarta Bay are mainly in an east to west direction during spring tides, with peak current speed occurring at the bay mouth and near the river mouths (Figure 2a,b) [1]. Inside the bay, the current speed is mainly below 0.2 m/s. The ebbing currents are slightly larger than the flood currents due to the river discharge. The current speeds near the bottom level reach 0.1 m/s in Jakarta Bay during spring flood tides (Figure 2c,d). During spring ebbing tides, the bottom currents are small (within 0.07 m/s) in almost the entire bay.

**Figure 1.** Map of the location of Jakarta Bay, and the Giant Seawall (GSW) as described in the NCICD Master Plan (Source: Master Plan NCICD). Numbers 1 to 13 correspond to rivers.

To overcome the flooding issues in Jakarta, Indonesia's government proposed the Giant Seawall (GSW) construction through the Master plan National Capital Integrated Coastal Development [8]. The GSW construction is considered through three phases: phase A, B, and C. In phase A the existing sea wall and river embankments are strengthened (the solid blue line in Figure 1), including land reclamation (the light grey color along the coast in Figure 1). In phase B, the outer sea wall is constructed and land is reclaimed (the dark grey color in Figure 1). This phase B is mainly determined by the required storage capacity of giant water reservoir between the coastline and seawall. There are two reservoirs in the GSW to adjust the water discharges from rivers. Phase C covers the long-term development in the east of Jakarta Bay (the dashed lines).

**Figure 2.** Currents (m/s) during spring tides in Jakarta Bay: (**a**) flooding and ebbing (**b**) ebb currents near the water surface; (**c**,**d**) show flooding and ebbing currents near the bottom level.

#### *2.2. Tides in Jakarta Bay*

The tidal level can be decomposed into different tidal components using harmonic analysis with the least squares method in MATLAB [20]. There are eight astronomical tidal components, named M2, S2, N2, K2, K1, O1, P1, and Q1 tides. When tides propagate into shallow water areas, the shallow water components are generated. For example, the MF, MM, M4, MS4, MN4 tides are all shallow water tides.

The Global Sea Level Observing System (GLOSS) by the Intergovernmental Oceanographic Commission (IOC) provides the sea surface level data for Kolinamil Station, Jakarta. The measurements at Kolinamil Station are from 2 June to 30 June 2015. The maximum tidal range at Kolinamil Station is 0.87 m. The maximum tidal range at Jakarta Bay is about 1 m with the highest level at around 0.6 m and the lowest level at around −0.4 m during spring tides relative to the mean sea level.

Defant (1961) [21] classified tidal type based on the amplitudes of tidal harmonic components called Formzahl number. The Formzahl number (*F*) is obtained from:

$$F = \frac{Amp\_{K1} + Amp\_{O1}}{Amp\_{M2} + Amp\_{S2}} \tag{1}$$

where *AmpK*1, *AmpO*1, *AmpM*<sup>2</sup> and *AmpS*<sup>2</sup> indicate the amplitude of the corresponding tidal component. The classification of the Formzahl number is as follow: *F* < 0.25 is semidiurnal type, 0.25 < *F* < 3 is mixed type, and *F* > 3 is diurnal type.

According to the Formzahl number (*F*) by Defant (1961) [21], the tidal pattern in Jakarta bay is diurnal indicated by the Formzahl number of 3.95 observed at Kolinamil Station. Hence, the tides in Jakarta Bay are diurnal.

#### *2.3. Description of the Flow Model*

The finite-volume coastal ocean model (FVCOM) was used to simulate the hydrodynamics of Jakarta Bay. The finite volume approach and unstructured meshes with 3-D primitive equations make FVCOM ideally suited for coastal areas [22]. The unstructured grid of FVCOM provides topography flexibility that is adequate for the irregular morphology of Jakarta Bay. A σ-coordinate transformation system in a vertical direction is used to get a clear representation of the irregular bottom topography. The σ-coordinate transformation is defined as:

$$
\sigma = \frac{z - \zeta}{D} \tag{2}
$$

where total water column depth is *D* = *H* + ζ, while *H* is bottom depth (relative to *z* = 0) and ζ is the height of the free surface (relative to *z* = 0).

Under the σ-coordinate system, the continuity and momentum equations can be written as Equations (2)–(5), and salinity, temperature, and density as Equations (6)–(8):

$$\frac{\partial D}{\partial t} + \frac{\partial uD}{\partial x} + \frac{\partial vD}{\partial y} + \frac{\partial w}{\partial \sigma} = 0 \tag{3}$$

$$\frac{\partial \underline{u}D}{\partial t} + \frac{\partial \underline{u}^2 D}{\partial x} + \frac{\partial \underline{w}D}{\partial y} + \frac{\partial \underline{u}\underline{v}}{\partial v} - fv\underline{D} = -\underline{g}D\frac{\partial \underline{\zeta}}{\partial x} - \frac{\underline{g}\underline{D}}{\rho\_0} \left[\frac{\partial}{\partial x}\left[D\int\_\sigma^0 \rho^\prime \partial \sigma^\prime\right] + \sigma \rho^\prime \frac{\partial \underline{D}}{\partial x}\right] + \frac{1}{D}\frac{\partial}{\partial v}\left(\mathcal{K}\_m \frac{\partial \underline{u}}{\partial t}\right) + DF\_{\underline{u}}\tag{4}$$

$$\frac{\partial v\mathcal{D}}{\partial t} + \frac{\partial w\mathcal{D}}{\partial x} + \frac{\partial v^2\mathcal{D}}{\partial y} + \frac{\partial w\omega}{\partial v} + fu\mathcal{D} = -g\mathcal{D}\frac{\partial\zeta}{\partial y} - \frac{g\mathcal{D}}{\rho\_0}\left[\frac{\partial}{\partial y}\Big|\int\_{\sigma}\rho'\mathcal{\partial}\sigma'\Big| + \sigma\rho'\frac{\partial\mathcal{D}}{\partial y}\Big| + \frac{1}{D}\frac{\partial}{\partial\sigma}\Big(\mathcal{K}\_{\mathfrak{W}}\frac{\partial\overline{\boldsymbol{w}}}{\partial\overline{\boldsymbol{w}}}\Big) + DF\_{\overline{\boldsymbol{w}}}\tag{5}$$

$$\frac{\partial wD}{\partial t} + \frac{\partial uwD}{\partial \mathbf{x}} + \frac{\partial vwD}{\partial y} + \frac{\partial ww}{\partial \sigma} = -gD\frac{\partial \zeta}{\partial \mathbf{x}} + \frac{gD}{\rho\_0} \Big[ \sigma \rho' \frac{\partial D}{\partial \sigma} \Big] + \frac{1}{D} \frac{\partial}{\partial \sigma} \Big( \mathbf{K}\_m \frac{\partial w}{\partial \sigma} \Big) + DF\_w \tag{6}$$

$$\frac{\partial SD}{\partial t} + \frac{\partial SuD}{\partial x} + \frac{\partial SvD}{\partial y} + \frac{\partial S\omega}{\partial \sigma} = \frac{1}{D} \frac{\partial}{\partial \sigma} \Big( \mathbf{K}\_{\hbar} \frac{\partial \mathcal{S}}{\partial \sigma} \Big) + DF\_s \tag{7}$$

$$\frac{\partial T D}{\partial t} + \frac{\partial T u D}{\partial \mathbf{x}} + \frac{\partial T v D}{\partial y} + \frac{\partial T \omega}{\partial \sigma} = \frac{1}{D} \frac{\partial}{\partial \sigma} \Big(\mathbf{K}\_{\mathbf{h}} \frac{\partial T}{\partial \sigma}\Big) + D\hat{H} + D\mathbf{F}\_T \tag{8}$$

$$
\rho = \rho(S, T) \tag{9}
$$

where *x*, *y*, and σ are the respective east, north, and vertical axes in the σ-coordinate; *u*, *v*, and ω are velocity components; *S* and *T* are salinity and temperature; ρ is total density; ρ is the reference density of ρ0; *f* is the Coriolis parameter; *g* is gravitational acceleration; *Km* is thermal vertical eddy viscosity and *Kh* is thermal vertical eddy diffusion; *Fu*, *Fv*, *Fw*, *FS*, *FT* are horizontal momentum, salt diffusion, and thermal terms.

The finite volume discrete method is applied in the model. The Mellor and Yamada level 2.5 (MY-2.5) turbulent closure scheme [23] was used for vertical mixing and the Smagorinsky scheme for horizontal mixing [24]. More details regarding the model can be found in Chen et al. (2006) [22].

#### *2.4. Flow Model Setup*

The computational domain covers Jakarta Bay, extending from 5.8◦ S to 6.15◦ S and 106.6◦ E to 107.1◦ E (Figure 3). This domain is covered by two grids, one with no GSW (Figure 3a) and one with GSW (Figure 3b). The model grid for Jakarta Bay with no GSW (Figure 3a) consists of 63,932 nodes and 127,273 elements, while the grid for Jakarta Bay with a GSW (Figure 3b) consists of 62,880 nodes and 122,986 elements. The computational grid has a high resolution varying from ~2.3 km at the open ocean boundary to ~25 m near the coastlines. Land areas of ~4.5 km from the nearest coastline have resolutions varying from ~25 m to ~300 m. Ten uniform vertical layers in a sigma coordinate system are set in the model according to sensitivity analysis.

**Figure 3.** Model Grid for Jakarta Bay (**a**) reference scenario 1 and (**b**) Jakarta Bay after the Giant Seawall (GSW) construction (scenario 3). The numbers indicate rivers. The blue dot and triangle indicate field stations.

Coastline data were extracted from the Navy Chart Indonesia and Google Earth (2014). Bathymetrical data were obtained from survey data provided by the Ministry of Public Works and combined with nautical charts of Jakarta Bay. Jakarta Bay is a shallow water bay with maximum depth around 60 m. Land topography was obtained from the Consortium for Spatial Information (CGIAR-CSI) (http://srtm.csi.cgiar.org). The SRTM (Shuttle Radar Topographic Mission) digital elevation data were used for this land topography with a high spatial resolution of 90 m. The CGIAR-CSI SRTM data product applied a hole-filling algorithm to provide continuous elevation values [25].

The observation data of tidal elevation and currents from two field stations (Figure 3) were used to validate the numerical model. The blue dot represents Kayangan Station, which provided tidal elevation and current data, while the blue triangle represents Kolinamil Station, which provided tidal elevation data.

The numerical simulation is driven by ocean tides, winds, and river discharges. The open boundary condition is used for ocean tides as the main forcing, while the closed boundary condition is used for river discharges as the main forcing. The wind data were used as the atmospheric forcing at this research. Other types of external forcing such as precipitation, evaporation, heat flux, and groundwater input were not included in the simulation. Initial conditions for tidal levels and currents were set to zero.

Oceanic tides play an important role when considering the flooding effect. The tide in the Java Sea is a diurnal type with maximum tidal range around 1 m during spring tides. Tidal forcing data were obtained from the TPXO (http://volkov.oce.orst.edu/tides/tpxo8\_atlas.html) [26]. Thirteen tidal constituents included eight basic components (M2, S2, N2, K2, K1, O1, P1, and Q1) and five shallow water components (MF, MM, M4, MS4, MN4).

Jakarta Bay itself is sheltered by many islands and so wind plays a very limited role [27]. However, a monsoon wind still impacts the dynamic circulation in Jakarta Bay. Wind data were derived from the European Center for Medium-Range Weather Forecast (ECMWF). ECMWF can be accessed at ERA-Interim reanalysis (http://apps.ecmwf.int). The spatial and temporal resolutions of the wind data were 0.125◦ × 0.125◦ and 6 h, respectively.

Thirteen rivers discharge into Jakarta Bay. The total mean river discharge from these 13 rivers was around 205 m3/s throughout the year 2012 [28]. Detailed information of the annual distribution of river discharge can be seen in Table 1.


**Table 1.** River discharge in Jakarta Bay based on Wulp et al. (2016) [28]. For rivers location see Figure 3.

Sea surface salinity in the western Java Sea varied from 30.6–32.6 PSU [29], while sea surface temperatures in northern Java ranged from 29–31 ◦C [30]. In this state, the initial salinity and temperature were set as constant values of 31.5 PSU and 31 ◦C, respectively. In this paper, water density was assumed to be homogeneous. The authors assumed that the Manning coefficient for bottom roughness varied from 0.004 at the open ocean to 0.008 near the islands and coastal region. The model ran for 17 days starting from 3 to 19 June 2015. Courant Friedrich Levy (CFL) was used as the numerical stability condition to calculate the external and internal step. The external time step numerical stability is defined as:

$$
\Delta t\_E \le \frac{\Delta L}{\sqrt{gD}} \tag{10}
$$

where Δ*tE* is the time step of the external model and Δ*L* (computational length scale) is the shortest edge from individual triangular grid element, and *D* is the local depth. The internal step is defined as:

$$
\Delta t\_I \le \frac{\Delta L}{C\_I} \tag{11}
$$

where Δ*tI* is the time step of the internal model, and *CI* is the maximum phase speed of the internal gravity waves. Since *CI* is usually smaller than *CE* = *gD*, Δ*tI* could be larger than Δ*tE*. The external time step was 0.15 s and the internal time step was 5 s based on the numerical stability condition.

#### *2.5. Scenarios*

The prevention of future coastal flooding in Jakarta is vital because of the phenomenon of future sea level rises and the frequent occurrence of high river discharge conditions. The Intergovernmental Panel on Climate Change (IPCC) reported that the temperature in Indonesia will increase from 1.3 ◦C to 4.6 ◦C. As the temperature increases, sea surface levels will also increase. The rate of sea level rise accelerated between the mid-19th and the mid-20th centuries based on tidal gauge and geological data [5]. The TOPEX/POSEIDON satellite is joint venture between NASA (National Aeronautics and Space Administration) and CNES (National Centre for Space Studies). The satellite observed that the rate of sea level rises along the coast of Indonesia was approximately 7 mm/year between 1992 and 2015 and that value was used in this paper. There are three scenarios for the sea level; the first is the scenario 1 (reference model), the second is with a sea level increase of 25 cm (scenario 2a), and the third is with a sea level increase of 60 cm (scenario 2b).

On the other hand, high river discharge occurs due to high precipitation and intraseasonal effects such as La Niña. High precipitation usually occurs under the northwest monsoon (NWM) conditions from December to February [13]. There are two high river discharge scenarios, 1.5 times the average discharge (scenario 1a) and 10 times the average discharge (scenario 1b). Scenario 1a is linked to the northwest monsoon (NWM) conditions, while scenario 1b is assumed for extreme river discharge. These scenarios for river discharge are combined with the scenarios for sea level rise to give a view of the areas affected by coastal flooding.

The impact of the GSW on coastal flooding was also included in this model. The effectiveness of the GSW in preventing the city from coastal flooding was examined using river discharges (scenarios 3a and 3b) and the sea level rises (scenarios 4a and 4b), presented in Table 2.


**Table 2.** Descriptions of different scenarios for Jakarta Bay. Scenario 1 is the reference scenario before the construction of the GSW and covers the 3–19 June 2015 period.

The numerical scenarios are listed in Table 2, showing Jakarta Bay (scenario 1) as the reference scenario before the construction of the GSW. This scenario covers a period from 3 to 19 June 2015.

#### **3. Model Validation**

The reference model was validated by comparing with the observation data. Data from two field stations, Kayangan Station (tidal and current) and Kolinamil Station (tidal), were used to validate the model. Data from the Kolinamil Station were obtained from http://www.ioc-sealevelmonitoring.org. Tidal and current data at Kayangan Station were collected by the Ministry of Public Works. Tidal data from the Kayangan Station were collected from 2 to 15 June 2015 and current data from 12 to 14 June 2015. During the June period, Jakarta experienced a dry season because of the southeast monsoons (SEM) conditions. Wind direction moved westward with maximum wind speeds at the simulation period around 4 m/s.

Ma et al. (2011) [31] suggested a statistical method to assess model skill for verification. This skill indicates the corresponding degree of deviation from prediction and observation data. A skill value of 1 means perfect agreement between model and observation.

$$Skill = 1 - \frac{\sum\_{i=1}^{N} \left| X\_{mod} - X\_{obs} \right|^2}{\sum\_{i=1}^{N} \left( \left| X\_{mod} - \overline{X\_{obs}} \right| + \left| X\_{obs} - \overline{X\_{obs}} \right| \right)^2} \tag{12}$$

Tidal verification data at Kolinamil Station started from 6 June 2015 to 19 June 2015 (Figure 4a), while the data at Kayangan Station were from 6 to 15 June 2015 (Figure 4b). The skill parameters for the tidal level calibration at Kolinamil Station and Kayangan Station were 0.97 and 0.96, respectively. The tidal data at Kolinamil Station and Kayangan Station were reproduced by the model with root mean square error (RMSE) of 0.08 m and 0.09 m, respectively. These two tidal verifications show good agreement between the model and observation data.

**Figure 4.** Comparisons of tidal level between the model results and the field data at (**a**) Kolinamil Station and (**b**) Kayangan Station.

Current verification was conducted at Kayangan Station, where the data were collected from 12 to 14 June 2015 (Figure 5). The skill scores between the modelled and the observed current velocity and direction were 0.96 and 0.88, respectively. The velocity and direction at Kayangan Station had RMSE of 0.015 m/s and 50◦, respectively. These skill values indicate that the model has good agreement with the observation data.

**Figure 5.** Comparisons of current speeds and directions between the modelled and observed data at Kayangan Station (12–14 June 2015).

#### **4. Results**

#### *4.1. Tidal Dynamics in Jakarta Bay*

Sea surface level data in the reference model were analyzed using harmonic analysis [20] to examine the main components of tides in Jakarta Bay. The four main astronomical tidal constituents (K1, O1, M2, and S2) explain about 80% of the total variances of water level in the bay. The diurnal components (K1 and O1) account for more than 70% of the tides, with tidal magnitudes of approximately 0.25 m and 0.13 m, respectively. The semidiurnal tides just account for a small part of the tides in the bay, with amplitudes of M2 and S2 tides of ~0.04 m and ~0.07 m, respectively (Figure 6).

Model results showed that diurnal tides and semidiurnal tides propagate differently. Diurnal tides mainly propagate westward, while semidiurnal tides propagate eastward. The different propagation pattern of diurnal and semidiurnal components is caused by the topography of Jakarta Bay, as the bay is located at the Java Sea and near to the Sunda Strait. Ray et al. (2005) [32] reported that the semidiurnal response is clearly dominated by the large tides from the Indian Ocean, while, in contrast, the diurnal tides come from the Pacific Ocean. These diurnal tides propagate westward in the Java Sea. The diurnal tides from the Java Sea propagate westward and meet waves coming from the South China Sea and the Karimata Straits [13,32,33]. The semidiurnal tides penetrate north and merge with the diurnal tides at western Kalimantan. In Jakarta Bay, the diurnal components are mainly driven by tides that come from the Pacific Ocean. The semidiurnal components are driven by tides that come from Indian Ocean. The shallow bathymetry and complicated coastlines generate the distribution pattern of the tides in the bay.

K1 and O1 tides share a similar pattern of propagating westward, from the Pacific Ocean through the Makassar Strait and the Molucca Sea. The semidiurnal components show another behavior. The eastward propagation of an M2 tide comes from the Indian Ocean through the Sunda Strait, while an S2 tide shows a more complicated pattern. An S2 tide comes from the Sunda Strait then propagates northward, and from the Karimata Strait then propagates southward. The complex patterns of tides in the bay are affected by the standing waves which occur in the west part of Java, and the amphidormic point occurring at the Sunda Strait.

**Figure 6.** Co-tidal charts of the main tidal constituents (K1, O1, M2, and S2).

In the reference model (scenario 1), the total inundation areas of Jakarta Bay amount to 32.36 km2, which is 7.7% of the total land area in the study domain (420.13 km2). These inundation areas relate to the land topography around the coastlines of the bay. Some parts of the topography are already under mean sea water level. Most of the inundated parts occur at the eastern part (near the river number 13) and central part of Jakarta (between the river number 6 and 7). In the eastern part, the Citarum River (river number 13) has the highest river discharge of all 13 rivers in Jakarta Bay. The topography of eastern Jakarta also has the same elevation as the coastal areas at around 0 m. In the central part of Jakarta Bay, the elevations of most of the inundated areas are below mean sea water level. The inundated areas are caused by sea water overtopping, as the rivers only have small amounts of river discharge.

Jakarta Bay, without GSW, had the highest water elevation of approximately 0.6 m at high spring tides. In scenario 1a, there is almost no change of water elevation in the entire bay. In scenario 1b, the water elevation along the coastline increased by less than 2.08% (1.25 cm) at the western part of the bay (from river numbers 1 to 6), and 0.9% (0.5 cm) at the center part of bay (from river numbers 6 to 11), and 3.33% (2 cm) at the eastern part of the bay (from the Keramat River to the Citarum River). In scenario 2a, there is a slight change of water elevation increase, 0.33% (0.2 cm) at the center and eastern parts of bay (from river numbers 6 to 13). In scenario 2b, water elevation increases by 4.17% (2.5 cm) at the center part of the bay (near river numbers 2 to 7) and reaches a maximum, 6.67% (4 cm), at the eastern part of the bay (near river number 13).

#### *4.2. Impacts of River Discharge on Hydrodynamics in Jakarta Bay*

The impacts of the river discharge in Jakarta Bay were tested in scenarios 1a and 1b. Wulp et all (2016) [28] reported that in wet seasons mean river discharge in Jakarta increased 1.5 times than the averaged river discharge. If scenario 1a were applied, the total inundation areas of Jakarta Bay would be 33.88 km<sup>2</sup> (8.06% of the total land areas in the study domain), increasing by 1.55 km<sup>2</sup> compared with the reference scenario 1. The inundation areas slightly increase near the eastern part of Jakarta where the Citarum River is located.

Under extreme conditions, when the river discharge is 10 times the normal one (scenario 1b), flooded areas in Jakarta Bay are 58.52 km<sup>2</sup> (13.93% of the total land area), an increase of 26.61 km<sup>2</sup> with respect to scenario 1 (Figure 7). The inundation areas mainly occur in the western and eastern part of Jakarta. A major cause of this extreme condition is the Cisadane and Citarum rivers (river numbers 1 and 13, respectively), which have larger river discharges than the others. The central part of Jakarta is relatively unaffected in this scenario due to the small amount of river discharge nearby. As the river discharge increases, the inundation areas become larger, especially in the eastern part of Jakarta.

To further examine the impacts of river discharge on tidal amplitude, the dominant astronomical tidal component K1 was chosen to represent the changes in the bay. In scenario 1a, the K1 tidal amplitude increases only near the river mouths (river numbers 11 to 13) by 0.8% (0.2 cm). In scenario 1b, the K1 amplitude increases by 1% (0.25 cm) in the western and eastern parts (river mouth numbers 1, 11, and 13) and 0.5% (0.125 cm) in the center part (river number 6) of the bay. There is no significant change for the phase lags in K1 in scenario 1a. In scenario 1b, the phases for K1 advanced, meaning that high tides would occur earlier than in the reference scenario 1.

**Figure 7.** Changes in K1 tidal amplitudes (cm) and inundation areas in Jakarta Bay between scenario 1b and scenario 1. Black dots indicate the inundation area in scenario 1, and red dots indicate the increased inundation area in scenario 1b. Color contours indicate the changes in K1 amplitude (cm).

Surface current directions in Jakarta Bay are largely westward during flood spring tides. The surface current pattern changes slightly due to the variations of river discharge (Figure 8). River discharge strengthens the current magnitude at the eastern part of Jakarta Bay and weakens the current magnitude at the central and western parts of Jakarta Bay. In scenario 1a, current velocities decrease by 33.33% (0.1 m/s) at the center of the bay. In scenario 1b, current velocities increase by around 50% (0.2 m/s) at the center of bay and by around 100% (0.4 m/s) at river mouths (numbers 1, 11, and 13).

**Figure 8.** Differences in surface current velocity (m/s) and direction under flood condition spring tide (**a**) scenario 1a and (**b**) scenario 1b. Black vectors are for scenario 1, and grey vectors are for scenario 2a/2b. Color contours indicate the change in magnitude (m/s) of the surface current velocity.

#### *4.3. Impacts of Sea Level Rise on Hydrodynamics*

The impacts of sea level rise on the hydrodynamics of Jakarta Bay were tested in scenarios 2a and 2b (Figure 9). In scenario 2a (0.25 m sea level rise), the total inundation areas were 43.37 km<sup>2</sup> (10.32% of the total land area), increasing by 11.01 km2 from the reference scenario 1. In scenario 2b (0.60 m sea level rise), total inundation areas were 62.04 km2 (14.77% of the total land area), with an increase of 29.68 km<sup>2</sup> from the reference scenario 1. The flood areas of Jakarta Bay occurred near the coastal region from the western to the eastern part.

**Figure 9.** Changes in K1 tidal amplitudes (cm) and inundation areas in Jakarta Bay between scenario 2b and scenario 1. Black dots indicate the inundation area in scenario 1, and red dots indicate the increased inundation area in scenario 2b. Color contours indicate the changes in K1 amplitude (cm).

When the sea level rise is 0.25 m (scenario 2a), the K1 tidal amplitude increases by around 1% (0.25 cm) at the eastern part of the bay (from river numbers 11 to 13) and decreases by about 2% (0.5 cm) between river numbers 10 and 11. There is a small phase lag change for the K1 component, compared with that in the scenario 1.

When the sea level rise is 0.60 m (scenario 2b), the increase of the K1 amplitude is about 1% (0.25 cm) near the coastlines at the eastern part of the bay, which is almost the same as that in scenario 2a. On the other hand, the K1 amplitude decreases between river numbers 10 and 11 by around 3% (0.75 cm), which is larger than that in scenario 2a. The K1 tidal phase remains almost unchanged. The K1 tidal amplitude increases with the rate of sea level rise.

The pattern of surface currents remains largely the same as the sea level rises (Figure 10). A greater rate of sea level rise strengthens the current magnitude in the entire bay, except in the middle part of the bay. Under a sea level rise of 0.25 m (scenario 2a), current velocities increase by 8.3% (0.025 m/s) in the entire bay and decrease by 13.33% (0.04 m/s) in the middle of the bay. Under conditions of sea level rise of 0.6 m (scenario 2b), current velocities would increase by 16.6% (0.05 m/s) near the coastal region and decrease by 24.9% (0.075 m/s) in the middle of the bay.

**Figure 10.** Differences in surface current velocity (m/s) and direction under flood condition spring tide (**a**) scenario 2a and (**b**) scenario 2b. Black vectors are for scenario 1, and grey vectors are for scenario 2a/2b. Color contours indicate the change in magnitude (m/s) of the surface current velocity.

#### **5. Discussion**

This research improves the understanding of sea level rise and river discharge effect on the Jakarta Bay hydrodynamics. As the Jakarta Bay tidal range is no more than 2 m, it can be classified as a microtidal region. Moreover, it is likely that the other microtidal regions with similar characteristics to Jakarta Bay will experience the same conditions under the effects of sea level rise. These findings can be used to assess future hydrodynamic conditions due to sea level rise effects. The projected change of tidal hydrodynamics and inundation areas in this study will be helpful to adopt sound management strategies and adaptation plans.

The impacts of the GSW on the Jakarta Bay flooding area are shown in Table 3. To overcome flooding disasters, the construction of the GSW has been proposed for Jakarta Bay. The construction of the GSW will slightly increase the maximum tidal range by 0.01 m [1] because of the reduced tidal prism by 20% and increased tidal choking in the bay.


**Table 3.** Flooding areas of Jakarta Bay due to sea level rise and river discharge.

In scenario 3, the total inundation areas of Jakarta Bay are 33.47 km2 (7.97% of the total land area), increasing by 1.11 km<sup>2</sup> from scenario 1. There would be a slight change in flooding area after the GSW had been constructed.

After the construction of the GSW, a slight change would occur to tidal amplitudes during high spring tides near the coastline, compared with scenario 1 (Figure 11). Tidal amplitudes would increase by around 1 cm (2% from scenario 1) at the eastern reservoir. At the western reservoir, tidal amplitudes near the coastal area would increase by around 3 cm (5.88% from scenario 1) and 1 cm (2% from scenario 1) near the channel on the left side of the western reservoir. There would no change in tidal amplitude on the right side of the western reservoir. Tidal amplitudes in the eastern part of the bay (from river numbers 12 to 13) will decrease by around 1 cm (2% from scenario 1).

**Figure 11.** Changes in water elevation (m) between scenario 1 and after the Giant Seawall (GSW) construction (scenario 3).

#### *5.1. E*ff*ects of the Giant Seawall on Sea Level Rise*

The effects of the GSW on sea level rise were tested in scenarios 4a and 4b (Figure 12a). The inundation area in scenario 4a would be 44.05 km2 (10.49% of the total land area in the study domain), increasing by 10.58 km<sup>2</sup> from scenario 3 (the reference model with GSW). In scenario 4b, the total inundation area would be 63.56 km<sup>2</sup> (15.13% of the total land area in the study domain) increasing by 30.09 km2 from scenario 3. The inundation areas in scenarios 4a and 4b would increase by less than 1% compared with scenarios 2a and 2b. The construction of the Giant Seawall in Jakarta Bay would need further design improvement to prevent a flooding hazard due to sea level rise. Furthermore, the implementation of the GSW design should also consider strengthening and increasing the current seawall height on the shore to prevent sea water overtopping the wall.

Compared with scenario 3, when the sea level rises by 0.25 m in scenario 4a, the K1 tidal amplitude slightly increases by around 0.2% (0.05 cm) in the eastern reservoir and on the left side of the western reservoir, and decreases by 0.65% (0.1 cm) near the Cisadane River (river number 1). There are small phase lag changes for the K1 component near the reclamation island areas outside the western reservoir.

When the sea level rises by 0.60 m in scenario 4b, the K1 tidal amplitude increases by around 0.4% (0.1 cm) in the eastern reservoir (from river numbers 9 to 11) and by 1.22% (0.3 cm) near the coastlines inside the western reservoir. The K1 tidal amplitude also decreases by 0.8% (0.2 cm) on the right side of the western reservoir, the eastern reservoir (from river numbers 7 to 8), and the mouth of the Cisadane, Keramat, and Citarum rivers (river numbers 1, 12, and 13). There is no significant change of phase lag for the K1 component.

#### *5.2. E*ff*ects of the Giant Seawall on River Discharge*

The impacts of a Giant Seawall on river discharge were tested in scenarios 3a and 3b (Figure 12b). When river discharge is increased 1.5 times (scenario 3a), inundation areas are 36.7 km2 (8.74% of the total land area in study domain) increasing by 2.82 km2 from scenario 3. The inundation areas in scenario 3a would increase by 0.6% compared with scenarios 1a. The K1 amplitude remains almost the same in the entire bay. There is a small increase of around 0.42% (0.1 cm) at river mouth numbers 11 and 13. There are no significant changes in phase lag for the K1 component in scenario 3a compared to scenario 3.

If river discharge is 10 times higher (scenario 3b), all the areas are flooded due to a blocking of the tidal inflow. Further design improvement is needed for the construction of the GSW to prevent a flooding hazard due to high river discharge, especially under extreme conditions (25-year return period). The designed dike height should be increased along the rivers to prevent overland flow from high river discharge. Protection of the inundation area after construction of the GSW under sea level rise and high river discharge can be considered in mitigation strategies. Overall, this research will help to improve coastal management in Jakarta Bay.

**Figure 12.** Changes in tidal amplitudes (cm) and inundation areas (red dots) in Jakarta Bay: (**a**) K1 tidal component between scenario 4b and scenario 1 and (**b**) K1 tidal component between scenario 3a and scenario 1. Color contours indicate the changes in amplitude (cm), and the black/white dashed lines indicate the phase change.

#### **6. Concluding Remarks**

The impacts of sea level rise and river discharge on the hydrodynamics of Jakarta Bay have been investigated using FVCOM. Sea surface level and current data at two field stations were used to validate the model. The model reproduced tidal elevation and current very well, with statistical skill parameters of 0.96, 0.96 and 0.88 for elevation, current velocity and direction, respectively. Numerical experiments were designed to examine the contribution of sea level rise and river discharge to the changes of hydrodynamics in Jakarta Bay.

The distribution of the amplitudes and phases of diurnal and semidiurnal tides in Jakarta Bay are quite different to other bays as the bay is located in the Java Sea and near the Sunda Strait. The diurnal tidal wave moves westward, while the semidiurnal tidal wave moves eastward. Diurnal tides come from the Pacific Ocean through the Makassar Strait and Mollusca Sea and propagate westward into the Java Sea. In contrast, semidiurnal tides come from the Indian Ocean through the Sunda Strait and propagate northward and from the Karimata Strait, propagating southward. The shallow bathymetry and complicated coastlines contribute to generating the distribution pattern of the tides in the bay.

Future sea level rise and river discharge variations will affect the hydrodynamics and coastal flooding areas in Jakarta Bay. The increase of river discharge amplifies the tidal components as well as the inundation areas near the river mouth due to increased fluvial forcing: the K1 amplitude increases by ~1% (0.25 cm), the current magnitude increases by 100% (0.4 m/s), and the inundation area increases by 26.61 km2 under the 10 times of mean river discharge. The rise of sea level increases the diurnal tidal component due to the increased tidal forcing: when sea level rise is 0.6m, the K1 amplitude increases by ~1% (0.25 cm) near the coastline and current magnitude increases by 16.6 % (0.05 m/s). The increasing sea level rise and river discharge would amplify the inundation areas and the tidal currents due to the increased tidal forcing and the river discharge amounts.

The GSW construction would slightly increase the water level by around 1 cm at the eastern reservoir due to the reduced tidal prism, while in the western reservoir, the water level would increase by up to 3 cm. If sea level rise occurs together with flooding, the K1 amplitude would slightly increase by less than 1% (0.2 cm), due to the reduction of the dissipation effects. The design of the GSW will need further improvement to prevent flooding hazards in the event of sea level rise and the effect of high river discharge. Strengthening the current seawall on the shore and increasing the dike height along the rivers would help to prevent a flooding hazard in Jakarta Bay. The outcome of this research improves our understanding of the sea level rise and river discharge effects in open type bays.

**Author Contributions:** Conceptualization, Y.X. and L.L.; Methodology, M.Y.S.; Software, M.Y.S.; Validation, M.Y.S., Y.X. and L.L.; Formal Analysis, M.Y.S., L.L.; Investigation, M.Y.S., Z.H. and L.L.; Resources, M.Y.S., L.L.; Data Curation, L.L.; Writing-Original Draft Preparation, M.Y.S.; Writing-Review & Editing, L.L. and Y.X.; Visualization, M.Y.S.; Supervision, Y.X. and Z.H.; Project Administration, Z.H.; Funding Acquisition, Z.H. and L.L.

**Funding:** This study was supported by the National Key Research and Development Program of China (2017YFC1405101) and was partially supported by the Natural Science Foundation of Zhejiang Province (LR16E090001) and the Bureau of Science and Technology of Zhoushan (2018C81036). Martin was supported by a scholarship from Zhejiang University for a Master's Degree program at the Ocean College of Zhejiang University, China.

**Acknowledgments:** We thank the Ministry of Public Works, Indonesia, for providing the tide, current, and bathymetry data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
