Coastal Impacts Driven by Sea-Level Rise in Cartagena de Indias

This work analyzes the coastal impacts of the combined effect of extreme waves and sea level extremes, including surges and projected mean sea level rise in Bocagrande, Cartagena (Colombia). Extreme waves are assessed from a wave reanalysis that are propagated from deep waters to the beach considering the hydrodynamic processes and taking into account the interaction between waves and the coastal elevation within the study area. First, we consider present sea level, storm surges and waves affecting the area. Next, we analyze the effect of sea level rise according to a moderate (RCP4.5) climate change scenario for the 21st century (years 2025, 2050, 2075, and 2100). The most pessimistic scenario (year 2100) yields a percentage of ﬂooded area of 97.2%, thus revealing the major threat that represents sea level rise for coastal areas in the Caribbean Sea.


INTRODUCTION
Low crested coasts and beaches are among the most vulnerable ecosystems on the Earth since the mean sea-level rise (SLR), land subsidence, and variations in the frequency and/or strength of marine extreme events have the potential to substantially affect their present morphology, uses, landscapes, and functions (Nicholls and Cazenave, 2010). The availability of these physical drivers at regional and even local spatial scales is crucial to quantify the most probable retreat of coastlines, which is essential to implement appropriate mitigation measures as well as for the development of adaptation plans (Hunt and Watkiss, 2011;Willis and Church, 2012;Lyu et al., 2014;Melet et al., 2018). This kind of information is being increasingly demanded by stakeholders and social actors in order to take the best possible decisions for present and future urban development (Hunt and Watkiss, 2011). Previous works have estimated the impact of SLR on coastal inundation at global scale (e.g., Nicholls et al., 2011;Hinkel et al., 2013;Passeri et al., 2015;Vitousek et al., 2017), and also at regional and local scales. Some examples are found in the United States (Limber et al., 2018), Portugal (Bon de Sousa et al., 2018), Spain (Enríquez et al., 2017(Enríquez et al., , 2019 or in the Netherlands (Verschuur, 2018) Previous results suggest that the retreat rates of coastlines could increase with respect to the present values, notably when the most pessimistic scenarios of greenhouse gases emissions are considered (e.g., Vermeer and Rahmstorf, 2009;Nicholls and Cazenave, 2010;Alexander et al., 2012). The total water level on the coast is the result of the combined effect of mean sea-level (SL) (with the addition of seasonal and inter-annual effects), low-level atmospheric pressures (surges), surface wind fields, tidal dynamics, wave breaking (set-up and run-up) and the topo-bathymetric heights (Nicholls et al., 2014;Guimarães et al., 2015). Mean sea-level variations result from changes in different contributors such as glaciers, ice sheets mass loss, land water storage or the ocean Frontiers in Marine Science | www.frontiersin.org 1 October 2019 | Volume 6 | Article 614 thermal expansion . Tide-gauge based observations corrected by vertical land motions and changes in the geoid show a trend in the global mean SLR of 1.1 ± 0.3 mm y −1 before 1990 for the 20th century, and of 3.1 ± 1.4 mm y −1 from 1993 to 2012 (Dangendorf et al., 2017). According to the last report from the Intergovernmental Panel on Climate Change (IPCC), the global mean sea-level will continue to rise during the 21st century and beyond , but the total amount will depend on the greenhouse-gas concentrations, with a projected global mean SLR for the high-end scenario close to 2 m in the 2100 Le Bars et al., 2017). Moreover, projections at regional scales can differ substantially from those for the global mean SLR , stressing the importance of performing studies that incorporate, where possible, the local SL contributions. The mean SLR increases the exposure of coastal zones to extreme marine events, causing the submersion of low-lying coasts, with obvious socioeconomic implications affecting human activities, coastal defense, development, growth and ecosystems (Nicholls and Cazenave, 2010;Jackson and Jevrejeva, 2016).
As it has been shown in the North Atlantic, the effects of extreme marine events on the coast also may increase due to changes in the storm wave climate (e.g., Masselink et al., 2016).
Based on altimetric data from recent decades Young et al. (2011) found an increment of wind speed and, to a lesser degree, of wave height on a global scale, although, these changes are difficult to predict on a regional scale since storm tracks may shift in the near future (Bengtsson et al., 2006). On a local scale, these changes have been documented by other authors in the north and northeast Atlantic Ocean and in the northeast and northwest Pacific Ocean, where the increase of wave height is associated with the upsurge in the intensity of the storms (e.g., Ruggiero et al., 2010;Masselink et al., 2016;Castelle et al., 2018). Potential changes in wind-waves have also been evaluated taking into account the last IPCC projections (Hemer et al., 2013;Casas-Prat et al., 2018). Apart from wave heights, projections also indicate some variations in the mean wave peak period and in the mean wave direction, which could have an enormous effect on sandy coasts since the beach planform depends on the energy flux direction (Vousdoukas et al., 2018).
Here, we evaluate to local scale the flood hazard in Bocagrande sand spit (Cartagena de Indias, Colombia) including all variables affecting the set-up and the run-up (i.e., the mean SL, tides, storm surges and waves) over the topography of the area. Taking into account all these variables is crucial to determine where and which mitigation actions have to be considered in front the impact of SLR. In order to account for future changes in the wave set-up, we will consider a moderate IPCC scenario corresponding to the Representative Concentration Pathway of 4.5 W/m 2 , hereinafter RCP 4.5 (Moss et al., 2010).

AREA OF STUDY
Cartagena de Indias is located in the northwestern South America within the Colombian Caribbean basin (Figure 1). It is a densely populated area and one of the most visited touristic destinations in the Caribbean region (Asociación Colombiana de Agencias de Viajes y Turismo [ANATO], 2018). Its beaches are mainly located in Bocagrande area (Figure 1, bottom-right panel). They conform a set of dissipative beaches controlled by groins spaced every 250 m approximately. These beaches, totally anthropized and without dunes are around 2.5 km length and between 30 and 70 m wide. The area is mainly exposed to offshore waves from the NE (Figure 2). The astronomical tide in the region is of mixed type primarily diurnal with a range of 0.40-0.60 m (Nicolae-Lerma et al., 2008;Restrepo et al., 2017) and the atmospheric-induced surges can exceed 0.2 m during strong storms . The mean SL for the period 1908-2009 in the Caribbean Sea shows a rise of 5.3 mm y −1 , about 4 times larger than the global mean SLR during the same period, whit around 5.6 mm y −1 in Cartagena for the period between 1950 and 2009 (Torres and Tsimplis, 2012). Moreover, local subsidence measured in the Cartagena region between 2000 and 2008 yielded a rate of 1.4 mm/year (Seemüller et al., 2009;Molares, 2011) thus contributing to exacerbate the mean SLR.
The regional climate is modulated by the location of the Intertropical Convergence Zone (ITCZ) as well as by the American Monsoon System, presenting two marked climatic seasons; the dry (or windy) season from December to April and the wet season from August to November (Poveda et al., 2002). During the dry season, the NE trade winds dominate the area due to the location of the ITCZ at the latitude 0-5 • N. During the wet season the trade winds from the south reach the Colombian basin due to the migration of ITCZ toward latitude 10-12 • N producing intense rain (Andrade, 1993). The rest of the year, from May to June (the transition season) is characterized by the weakening of the trade winds. Wave climate is modulated by these seasons presenting a bimodal distribution; during the dry season the significant wave height (H s ) present larger values and smaller in October (Osorio et al., 2009;Thomas et al., 2011;Andrade et al., 2013). Cold fronts during the dry season are considered the major cause of extreme waves in the central sector of the Caribbean littoral of Colombia (Otero et al., 2016), affecting periodically the Bocagrande area with energetic northwestern waves. From  June to November, the passage of hurricanes near the Colombian Caribbean Sea can be accompanied also by the occurrence of extreme events (e.g., as those caused by hurricanes Joan in 1988 and Lenny in 1999), which have the potential to severely affect the continental Caribbean (Ortiz et al., 2015).

DATA AND METHODS
The mean wave set-up in present day conditions is obtained from a tide-gauge located in Cartagena de Indias (current mean SL), while future conditions under the RCP 4.5 scenario are derived from a global model mean ensemble . ; timeseries of significant wave height (Hs) from the WWIII reanalysis at the closest grid point with respect to Cartagena (B); Kendal correlation between skew-surge and Hs (C) for the period indicated as joint analysis in the top panel (the black cross remarks significant correlation at 95%, see legend). The blue shading depicts the timeseries used to estimate the extreme events in the independent assessment; return level for skew-surge using 2 max. per year (D). The blue shading reflects the 5-95% range of uncertainty.
Other contributions to the wave set-up such as astronomical tides and sporadic extreme surges induced by storms are also computed from the above tide-gauge. We consider that waves propagate over the mean wave set-up, which is the result of adding to the mean sea level the astronomical tides and extreme surges obtained from a nearby tide-gauge located in Cartagena de Indias. Extreme deep water-wave conditions from WaveWatch III (WWIII) are propagated to the coast combing SWAN and SWASH (a non-linear shallow water solver) models to finally obtain the run-up.

Sea Level Observations and Projection
Hourly time series of sea level from tide-gauges located at Cartagena were obtained from http://gesla.org/ (Woodworth et al., 2017). Since tide-gauges records cover different periods (11/1951-4/1993 and 5/1993-12/2013), time series were homogenized, detrended and unified to have a unique long record. Astronomical tides were estimated with the UTide MATLAB software (Codiga, 2011). The analysis of the sea level data confirmed that the amplitude of astronomical tides in the region is around 0.2 m (microtidal). In order to reduce timing errors and non-linear effects we analyzed skew-surges over the sea level time series, which are defined as the difference between the maximum SL and the maximum predicted tidal level over one tidal cycle (Mawdsley et al., 2015). Projections of SLR according to the RCP 4.5 scenario for the Cartagena de Indias area were downloaded from the Integrated Climate Data Center (IDC) of the University of Hamburg 1 . This global database include different geophysical sources that drive long term changes in the relative sea surface height (SSH), such as, ice components, ocean-related components (derived from CMIP5 models), land water storage and glacial isostatic adjustment. From the above data a set of mean wave set-up is constructed by adding to the projected RCP 4.5 sea level scenarios for years 2025, 2050, 2075 and 2100, the spring astronomical tides and extreme storm surge for a return period of 10 years derived from the skew-surges time series (further details on this are provided in section "Selection of Extreme Events of Skew-Surges and Waves").

Wave Data
Extreme wave climate was obtained by performing a statistical analysis of the NOAA hindcast over a period of 38 years 2 . This hindcast is based on the WWIII model forced by the Global Forecast System (GFS) analysis winds for the period of 1979-2016 with a 3-hourly temporal resolution. Here we use the closest model grid point to Bocagrande beaches, located around 15 km off the coast of Cartagena (see Figure 1, top-right panel).

Modeling Set-Up
Based on the statistical analysis of the NOAA hindcast referred above, two energetic sea states were selected in deep waters that were characterized by their H s , D p (peak direction) and T p (peak period). These sea states were propagated using the SWAN model in stationary mode to obtain the corresponding wave conditions at shallow waters. These dispersive waves are transformed up to the beach using SWASH model in non-stationary mode. Present SL as well as those given by the moderate RCP4.5 scenario are analyzed in order to obtain the respective target inundations. The SWAN model solves the equation of action balance for the propagation of the wave spectrum, allowing realistic estimations of wave parameters in oceanic and coastal zones (Booij et al., 1999). To assess the performance of SWAN in the area, the model is first validated in a non-stationary simulation using the Japanese 55-year wind Reanalysis, JRA-55 (Kobayashi et al., 2015). Sea states from a model grid point at the same position of the NOAA buoy # 41194 belonging to the Dirección General Maritima (DIMAR) were used for the validation of the simulated waves. This buoy is 126 km NW off Bocagrande at 11.161 • N, 74.681 • W (see location in Figure 1, left panel). Wave growth by wind was set as exponential following the formulation of Komen et al. (1984), and the deep water nonlinear interaction by using the Webb-Resio-Tracy method. The relevant processes of interest: wave braking, energy dissipation by whitecapping and bottom friction were took into account for the simulation. The time step was set as 30 min and results stored every 3 h. The bias and root mean square error (RMSE) between the buoy and the SWAN model grid point time series are close to zero whereas the correlation coefficient is 0.84 showing a good fit between modeled and observed data. This relationship is statistically significant at a confidence level of 95% according to a two-tailed t-Student distribution (p-value <10 −5 ) (Figure 3).
The SWAN domain had a resolution of 100 × 100 m 2 (471 × 481 grid points) with bottom left UTM coordinates at 1136524.46N, 395970.14E (Figure 1 Swan grid in the right panel). Bathymetry was obtained from the nautical charts COL 042, COL 261, COL 263, COL 615 from the Colombian Hydrographic Institute, with scale ranges of between 1:15000 and 1:250000.
Regarding SWASH, the model solves the non-linear shallow water equations including the non-hydrostatic pressure being suitable for wave propagations up to beaches . SWASH was executed on a fine mesh with a spatial resolution of 3 × 3 m 2 (341 × 551 grid points) and with bottom left UTM coordinates at 836992N, 1641707E. For this domain the bathymetry was obtained from using a single-beam ODOM Hydrotrac 2 echo-sounder merged with a high-resolution LIDAR topography for the emerged part of the domain (Figure 1 Swash grid in the bottom-right panel). All vertical heights in the Digital Terrain Model (DTM) were referenced to the mean low water spring (MLWS) and the horizontal coordinates were referenced to UTM 18N. According to Smit et al. (2014), to accurately resolve waves in a phase-resolving model the mesh resolution has to be at least 10 times smaller than the shortest wavelength to be resolved i.e., L/ x = O(10) [similarly this applies to the time step i.e., T/ t = O(10)], thus satisfying this criteria with Frontiers in Marine Science | www.frontiersin.org

Video Imagery
Video monitoring systems developed after the apparition of digital cameras, have shown to be a powerful and low-cost tool to monitor the coast in a wide range of studies, such as, coastal variability (Nieto et al., 2010;Simarro et al., 2015), intertidal bathymetry (Aarninkhof et al., 2003) or evolution of beach systems (Ojeda and Guillén, 2008). Since video systems have the advantage of providing continuous monitoring of coastal areas they are also a valuable tool to validate numerical models (Salmon et al., 2007;Andrade et al., 2013;Morales-Márquez et al., 2018;Osorio et al., 2019). In this work, we have used video images from the HORUS monitoring system 3 at Bocagrande, to validate the run-up provided by SWASH model. The images passed through a scale rectification process and merged with digital orthophotography for visualization purposes.

SELECTION OF EXTREME EVENTS OF SKEW-SURGES AND WAVES
In this section we evaluate if the extremes of skew-surges and waves are significantly correlated. Our purpose is to estimate return level probabilities of both types of extremes. Based on this result we use the most suitable theoretical framework to estimate the magnitude of the skew-surge and waves that will be employed in section "Results" to model the coastal run-up.

Selection of Independent Events of Surges and Significant Wave Height (H s )
To study the correlation between time series of skew-surges (section "Sea Level Observations and Projection") and H s (section "Wave Data"), we selected the maximum of skew-surges and H s as follows; For each maximum value of skew-surge the associated maximum value of H s within an interval of ± 6 h is selected, thus having a set of pairs of both variables. To be sure that the selected skew-surges corresponded to different events, they had to span at least 3 days between them (see Torres and Tsimplis, 2014 for further details). Since the location of the ITCZ has a large influence on the waves in the Caribbean Sea there are only two seasons to be analyzed; the dry season from December to April and from July to August and the wet season from April to May and from August to November. The common time range between skew-surges and waves extends from July 1979 to June 2001, accounting for a total of 21 years (see Figures 4A,B).

Correlation Between Skew-Surges and H s
If a peak of skew-surge is correlated with a high wave, the probability of having a hazardous total water level will increase (the return period will decrease) respect to the case where both events are independent. In the former case, if the probability is computed independently, the actual risk would be underestimated. To avoid this, first we will assess the correlation between extremes of skew-surges and waves. Then, correlation is assessed with the Kendall rank correlation. Previous studies in the US coast and globally have suggested that an indicative value of 0.2 may be already significant enough (Wahl and Chambers, 2014) to require the computation of the joint probability using a multivariate approach such as the Copula functions (e.g., Sayol and Marcos, 2018;Enríquez et al., 2019;Marcos et al., 2019). Correlations for those independent events are shown in Figure 4C for all percentiles being in all cases very low with Frontiers in Marine Science | www.frontiersin.org 6 October 2019 | Volume 6 | Article 614 the maximum correlation of 0.175 corresponding to the 95percentile. The cross over the circles indicate the correlations that are significant at the 95% level according to a series of permutation distributions (Gibbons and Chakraborti, 2011). In the light of these results, it is reasonable to assume that the extremes of skew-surges and H s are uncorrelated and hence we compute their extremes independently.

Selection of Extreme Skew-Surges
The largest available time series for skew-surges (from July 1952 to June 2001) is used to compute the return levels using the block maxima approach. For each defined year (from July to June) we seek for the highest N independent events of skew-surges (they must be separated at least 3 days) collecting a set of N × M maxima, where M is the number of years. In order to estimate the return levels associated to a given probability, the extremes are fitted separately using the Generalized Extreme Value-GEV distribution (Tsimplis, 1995;Marcos et al., 2009). After checking the return levels for the range N = 1 to N = 5 we did not find a significant difference among them and thus we computed the return level with blocks of N = 2. As a reference we choose a return period of 10 years, which is at the same time extreme and statistically probable. The resulting return level for skew-surges found under this assumption is 0.2 m ( Figure 4D).

Selection of Extreme Waves
Time series of deep-water H s , T p, and D p (section "Wave Data") are here used to characterize wave conditions (Figure 2). Based on this almost 93% of waves come from the NE, 2% from the NNE, 1.5% from ENE and the remaining from other directions (see Table 1). After considering the direction with the highest probability of occurrence and the joint probability analysis of H s and T p , two cases representing extreme events were chosen for the NE direction: the significant wave height that exceeded 12 h per year (H s12 , cases A1-A2 in Table 2) which corresponds to H s = 4.5 m and T p = 11 s; and (ii) the significant wave height with a 10% probability of exceedance (H s90 , cases A3-A4 in Table 2) corresponding to H s = 3.0 m and a T p = 12 s. As the passage of cold fronts is the cause of some of the highest waves in the central of the Colombian Caribbean Sea (Otero et al., 2016), we also study (although they have a lower probability of occurrence) the extreme waves from the NW direction. For this case, considering the joint probability analysis of H s and T p , corresponds to H s = 3.0 m and T p = 14 s (H s12 , cases B1-B2 in Table 2) and to H s = 2.2 m and T p = 10 s (H s90 , cases B3-B4 in Table 2).

RESULTS
Extreme waves at deep waters have been propagated toward the coast following the workflow outlined in Figure 5. Two type of extreme wave conditions have been identified to be the most probable in terms of H s12 and H s90 corresponding to incoming waves from the NE and those resulting from the passage of cold fronts (NW). The remainder of this section describes the validation process followed for the models used to simulate waves in deep and coastal ocean areas. Moreover, coastal flood maps for Bocagrande and surrounding beaches under presentday conditions and for the RCP4.5 sea level projection for the set of time horizons are presented and discussed.

Model Validation
Run-up validation was performed using the images obtained with HORUS video monitoring. As an example, Figure 6 shows the HORUS image corresponding to November 8th of 2010, coinciding with the passage of a cold front in the Colombian Caribbean. Wave conditions for this event were H s = 2.3 m and T p = 9 s from the NE direction. All comparisons were done using ArcGis 10.1, in which the maximum run-up from both data sets were plotted on the coast. Figure 6A shows the maximum run-up manually detected over the georeferenced image (red line) while Figure 6B shows the maximum run-up given by the numerical model (blue line). Along the m-n transect, the longitudinal distance from the coastline to the maximum run-up obtained by the numerical model and the HORUS image was 61 and 72 m, respectively ( Figure 6C). The difference in the elevation between the simulation and the image run-up was 0.3 m, with a higher run-up measured at the image. In other areas of the transect, the differences were lower, in both the horizontal (<3 m) as in the vertical (<0.2 m) ( Figure 6C). As seen in Figure 6C, the run-up from the modeling and from the observations are in very good agreement and thus we used the model approach as a proxy for the projections.

Modeled Coastal Inundation
Based on the analysis of the deep water wave conditions and taking into account the statistical independence of the wave climate and storm surges, we have analyzed 12 cases for the deep water mean direction (NE) as well as for the front cold passages (NW) ( Table 2). Eight cases correspond to present-day conditions (SLR = 0, A1 to B4) and four additional cases under the projected SLR for years 2025, 2050, 2075, 2100 under the RCP 4.5 (C1-C4, Table 2). The flooded areas obtained for all case studies are presented and discussed below.

Flooding Under Present-Day Conditions
Present day conditions correspond to the combination of the deep water wave analysis outlined in section "Selection of Extreme Waves" with the inclusion (or not) the extreme skewsurge (section "Selection of Extreme Skew-Surges") ( Table 2).
Note that for all cases-studies astronomical spring tides of 0.2 m are assumed. The corresponding flooded areas for deep water waves from the NE (A1-A4) and from the NW (B1-B4) are shown in Figures 7, 8 and are summarized in Table 3 (A1-A4 and B1-B4, respectively). The most severe flooding for the present state is given for those cases that account for the combined effect of waves and skew-surges, regardless of the wave direction. This is shown in Figures 7A2,A4, 8B2,B4). H s90 for case A2, which corresponds to NE waves (case A2 in Table 2) correspond to wave conditions of H s = 1.5 m and T p = 12 s over a skew-surge of 0.2 m yields a flooded area of approximately 0.062 km 2 (A2 , Table 3), whilst the H s12 conditions for NE waves (case A3 in Table 2) correspond to wave conditions of H s = 2.3 m and T p = 11 s without skew-surge and results in a flooded area of 0.048 km 2 . The worst scenario for the NE waves in terms of flooded area is given under the case A4 (the most extreme waves and skew-surge) resulting in a flooded area of 0.087 km 2 . In summary, the flooded area for the A1, A2, A3, and A4 cases comprise 9.8, 20.4, 15.7, and 28.4% of the surveyed area, respectively (Figure 7 and Table 3).
For the NW direction (cold front passage), which corresponds to cases B1-B4 in Table 2, the flooded area is smaller compared with NE cases despite that the significant wave height in shallow waters is larger (Figure 8). The H 90 wave conditions (case B1) correspond to H s = 2.3 m and T p = 10 s with no skew-surge, results in a flooded area of approximately 0.014 km 2 . Same wave conditions over a sea level with an extreme skew-surge (case B2) yields a flooded area of approximately 0.052 km 2 . The H 12 conditions correspond to H s = 3.2 m and T p = 10 s (cases B3 and B4 without and with skew-surge, respectively) yield flooded areas of approximately 0.030 and 0.068 km 2 , respectively. For all cases (B1-B4 in Table 3) the flooded area comprises the 4.5, 17.0, 9.6, and 22.3% of the total surveyed area, respectively (Figure 8).

Flooding Under Projected SLR
For the future global warming scenario, we assume the worst present flooding situation given by the case A4 that corresponds to the H 12 from the NE direction (Table 3). Here we consider for all cases a skew-surge of 0.2 m, a high astronomical tide (0.2 m) and the regional SLR from RCP 4.5 for years 2025 (0.11 m), 2050 (0.24 m), 2075 (0.38 m) and 2100 (0.52 m) ( Table 2). We note that the SLR scenario here considered is one of the most conservative. The corresponding flooded area is displayed in Table 3 comprising the 38.7, 53.9, 73.7, and 97.2% of the surveyed area, respectively (Figure 9). It indicates that for the short term simulation (C1, year 2025) the flooded area will increase around 10% with respect to the present-day situation. For the remainder   Table 3.

DISCUSSION AND CONCLUSION
In this work we study the coastal impacts of the combined effect of extreme waves and sea level extremes in Cartagena (Colombia). A summary of the flooded area for the present day conditions and for the projected RCP4.5 scenario is presented in Figure 10. Sea level has been rising since 1870, mainly due to the global warming that melts of land ice and causing ocean expansion (Church and White, 2006). One of the immediate consequences of SLR is land submergence and a higher risk of coastal flooding (Nicholls and Cazenave, 2010). However, the degree of affectation of the coastal areas depends on the regional and local features (profile slope, topographic heights, sediment type, wave conditions, tide conditions, meteorological conditions, etc.). Therefore, to take the most appropriate mitigation and adaptation measures, it is necessary to perform local studies accounting for all physical variables to determine the impact of SLR over the coastal zones. Sayol and Marcos (2018), Their results show that beaches would suffer a coastal retreat between 7 and 49 m under RCP4.5 and RCP 8.5 sea level rise projections, respectively, which is equivalent to half of the present-day aerial beach surface. This last study, present lesser affectation over above mentioned studies because de beach profile has a steeper slope. Flooding in Bocagrande beaches is also related with their low topographic elevation as well as with the discontinuous elevations of the berm with values below the actual sea level in some areas. This fact, together with the coupled to the specific response of the Caribbean Sea under the future scenarios (fast thermos-steric response) makes it a highly risky spot for the SLR. Moreover, the poor correlation between extremes of skew-surges and waves supports the major relevance of swell waves in the generation of flooding in this part of the Caribbean region. Therefore, in order to assess the effect of the projected mean SLR under the RCP 4.5 scenario combined with the tides, high energy waves, surges and topo-bathymetric heights at the Bocagrande beaches we use the coupled models SWAN-SWASH.
Since non-stationary numerical simulations are computationally extremely demanding, we have evaluated the risk of flooding for a set of specific events. In this regard, after a statistical assessment of historical extreme skew-surge and wave conditions, the most representative case-studies have been selected. Besides we have considered the potential effects of SLR in Bocagrande beaches for the most hazardous present-day considered scenario (A4), in an attempt to optimize Frontiers in Marine Science | www.frontiersin.org 13 October 2019 | Volume 6 | Article 614 the computational time. Despite this limitation and in contrast to other methodologies previously used to assessed coastal inundation in Cartagena beaches and its surrounding areas (Nicolae-Lerma et al., 2008;Andrade et al., 2013;Nicolae-Lerma et al., 2013) the coupling of SWAN-SWASH models has allowed us to determine the extreme wave run-up over the Bocagrande beaches in two dimensions, with a good agreement when tested against measured wave data and video images. With this methodology we have taken into account the mean SL, high astronomical tides, and extreme skew-surges and waves, as well as the projected mean SLR under the RCP 4.5 scenario and the actual structural elements (groins, breakwaters, dams, etc.) that interact with the wave field. Therefore, this approach can be used to search for specific solutions (either soft and/or hard engineering approaches) to mitigate the effects of the SLR in the studied area. Indeed, a good representation of waves (H s and direction) is mandatory to adequately study coastal flooding. These processes can only be taken into account with the use of numerical models. From the above results, although the shallow water significant wave height is larger for the coldfront configuration (NW direction at deep waters, B1-B4 in Table 2), the inundation resulting from the NE extreme waves is higher due to waves refraction and location of the coast line. Because the coastline of the study area is parallel to waves from NNW direction (see Figure 1, bottom-right panel), NE waves have a direct impact on the beaches. On the other hand, incoming waves from the NW direction (WNW direction in shallow waters), had impact directly over the groins, thus not affecting beaches. The combined effect of astronomical tides and extreme skew-surge contribute to coastal inundation only between 25 and 45% of the projected SLR under RCP 4.5 (for year 2025 and 2100, respectively). Therefore, SLR appears to be the most dominant factor in this study area. We note that in this work it has been analyzed a conservative climate change scenario, hence even a higher importance of SLR can be expected for more pessimistic projections. This result is inferred from Figure 9 (cases C1-C4), where the inundated area increases dramatically as we move closer to the end of the 21st century. The main reason is the low topographic elevation, which in some cases are directly below the mean sea level.
Finally, we want to remark that among the factors that may influence the flooding in Bocagrande beaches, we have not considered the geodynamic processes, the changes in the beach profile as the sea level rises, the change in beach planform in response to a potential variation in the mean wave energy flux and eventual changes in the wave dynamics related to variations in the atmospheric patterns as a consequence of climate change. These limitations imply that; (i) for the year 2100 the total flooded area could be higher than the one provided in this study. According to the analysis of Andrade-Amaya et al. (2017) from a geodetic station located in Cartagena de Indias, this zone present subsidence processes with values between −1.78 ± 0.4 and −1.88 ± 0.44 mm/year. (ii) The response of the beach profile to sea level rise and wave dynamics also could have an influence on the value of the flooded area provided for each of the scenarios if there were a constant sediment flux and if the beach front would not be totally urbanized. Considering the latter, to limit the number of variables used in this work, we have assumed that the beach profile remains inalterable under the proposed scenarios which, as a first approach, is a reasonably assumption since changes in beach profile are below 20 cm under the worst SLR scenarios (Enríquez et al., 2017). (iii) Flood processes generated by the run -up under wave energetic conditions over Bocagrande beaches could result in an underestimation of the percentage of flooded area obtained in this study, especially for the year 2100 SLR scenario. According to Mentaschi et al. (2017) and Casas-Prat et al. (2018) wave energy conditions will probably change globally in response to a shift in the wind regimes because of the undergoing global warming.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
This work was part of AO-R's Ph.D. thesis. AO-R wrote the manuscript and performed the simulations. AO and JR were the supervisors and contributed to conceiving the topic and the development of the argument, interpretation of the results, and writing. JS, MM, LO, and IH-C contributed to the interpretation of the results and writing.

FUNDING
This work was supported by COLCIENCIAS (Departamento Administrativo de Ciencia, Tecnología e Innovación) through of convocatoria 727 and from the Spanish Government MICINN/FEDER through projects MUSA (CTM2015-66225-C2-2-P) and MOCCA (RTI2018-093941-B-C31). The authors acknowledge the financial support from the Spanish Government MICINN/FEDER through projects MUSA (CTM2015-66225-C2-2-P) and MOCCA (RTI2018-093941-B-C31). AO-R is supported by a grant from Colombian COLCIENCIAS (Convocatoria 727). JS thanks the financial support by the Netherlands Scientific Research Foundation (NWO) through the VIDI grant number 864.13.011 awarded to C. A. Katsman. IH-C acknowledges the Vicenç Mut contract funded by the Government of the Balearic Island. This work was partially performed while AO was a visiting scientist at the Earth, Environmental and Planetary Sciences Department at Brown University through a Ministerio de Ciencia, Innovación y Universidades fellowship (PRX18/00218).