Coastal Risk and Water Flow Analysis in Eastern Algeria (Western Mediterranean)

In this paper, we present the combination of tsunami modelling studies with a GIS approach to analyse the coastal flooding risk in eastern Algeria. The Jijel and Bejaia regions are well-known for their seismic and tsunami risk. The modelling was carried out with the Geoclaw software. It provided the water flow velocities in addition to water heights results for several points of interest located between the Bejaia and the Jijel coast. The simulations were based on a 7.6 magnitude earthquake offshore the Jijel bay. Maps visualizing the susceptibility to flooding were generated in ArcGIS using the GIS integrated weighted overlay tools. By aggregating and weighting morphometric factors (such as the lowest height levels and lowest slope gradients) influencing the susceptibility to flooding the coastal areas prone to inundations in case of high energetic flood waves (tsunami, storm surge, meteo .tsunami) from the sea side can be determined. The entrances of river mouths are the locations marked by the highest susceptibility to potential


Introduction
The Geodynamic setting in the western Mediterranean often leads to moderate to destructive earthquakes [1]. Sometimes, these earthquakes trigger slumps, slides and tsunami waves. Located at the boundary between the African and the Eurasian plates, the Bejaia-Jijel region is subjected to a convergent motion with a slow convergence rate that ranges about 4 to 6 mm/years. In 1856, on the 22 nd and 23 rd of August, a tsunami triggered by a destructive earthquake flooded the eastern coast of Algeria. In particular, 5 meters for Bejaia and 2 to 3 meters for the Jijel harbour have been reported in historical documents [2]. Inland, the Kherrata region, south of Bejaia is frequently the location for moderate seismic sequences with sometimes landslides that cause damage in rural areas [3,4]. Geological field investigations have been carried out in the eastern coast and highlighted evidences for past destructive events (tsunami deposits, boulders, seismites) [4,5].
Most often, the understanding of the tsunami generation mechanism at the origin of a coastal disaster is complex because of (1) the behaviour of the fault rupture and the identification of the source parameters and (2) the combination of the slides evidenced by the rupture of undersea cables well imaged by swath bathymetry and seismic data (Fig. 1). The 2011 Tohoku tsunami is a good example that shows the complexity to assess a tsunami source with accuracy in real time. It is quite common to constrain the sources parameters from a mechanism entirely due to slip on an earthquake fault. When the earthquake generated the tsunami, scientists quickly collected data and worked on seismic inversions to assess the tsunami source model. Nevertheless, Tappin et al. [6] revealed that the tsunami from the earthquake-only source cannot generate the high-frequency components seen in the observed data. Moreover, the tsunami inversion models based on tsunami waveforms recorded at nearshore and offshore GPS buoys failed to reproduce the run-ups data and the tsunami inundation penetration along the Sanriku coast. Tappin et al. [6] concluded that an additional tsunami source had to be considered as the high-frequency components could be reproduced by a submarine mass failure (SMF) that has been further identified at the northern boundary of the earthquake rupture zone. Finally, Maercklin et al. [7] developed a method to track the rupture process of the Tohoku earthquake from a back projection procedure. They reconstructed the slip distribution in space and time and the high resolution images revealed three localized slip patches and evidenced that a twin ruptures grew to build up the giant Tohoku earthquake with distinct seismological parameters (depth, seismic energy radiated and magnitude). Consequently, the modeling of the tsunami triggered at the onset of the first shock is very challenging as for the behavior of the rupture area (earthquakes and associated slides). The Algerian coast (North Africa) hosts about ten harbours dedicated to the maritime traffic (containers, oil terminal, storage tanks, vessels, ships, etc.). Urban nearby cities and marine and terrestrial ecosystem and livelihoods are exposed to an additional environmental risk potentially triggered accidentally by the subsequent earthquake or slides related damage. Evaluations of high resolution satellite imagery can be used for a contribution to a better understanding of the interaction of incoming waves with the coastal morphology and its influence on streaming dynamics. The data mining is aiming at visualizing critical areas and providing information about damage in case of emergency due to flooding hazards as fast as possible, as the civil protection units need this information for their management. The almost actual inventory of land use and infrastructure of settlements and cities is an important issue for the hazard assessment and damage loss estimation. Near-shore bathymetry and the coastal morphology might be leading to local resonance phenomena (see Fig. 1), for example as they have a great influence on the inundation extent of tsunami waves. In 2003, the damage reported in the Palma de Majorca bay after the 2003 Boumerdes earthquake were related to high frequency oscillations of relatively large amplitude waves and combined effect of the tsunami and local resonance [8].
The Bejaia harbor is mostly dedicated to the oil commercial activities and maritime transport ( Figure 2). Coastal risk is mainly related to the consequences of the harbor potential damage in case of an earthquake, a storm or tsunami waves for the region. Within the last centuries, several reports and studies showed the impact of tsunami and earthquakes events offshore on the vessels. This has been the case for the 1790 Tsunami and Earthquake in Oran in western Algeria (references). The moderate earthquake that triggered a local tsunami in Mostaganem probably because of related slides due to submarine canyons is as well a threat for industrial harbors located nearby ecosystems in Arzew [9]. The 1856 Jijel earthquake was largely described in historical documents with comments on the flooding of the Jijel harbor [2]. Today, nearby Jijel, a harbor is as well in expansion, in Djendjen, to provide more space for the trade market in the west Mediterranean. The growth of harbors to support the economy supply is in constant progress [10]. During the last decades, the Bejaia Mediterranean Terminal has been developed to handle the worldwide current maritime traffic (amount, size of the containers, etc. ) by Protek, the Singapore port operator (https://www.portek.com/port-operations/algeria-bejaia-me diterranean-terminal/). Closed to Jijel, the djendjen harbor is as well in development to handle large containers dedicated to the trade market and facilitate the Euro-Mediterranean trade market (Dubai Port World's Group). Coastal damages or impacts due to water waves are not only related to the height reached by the waves but also to the velocity and the energy that characterize the currents. The knowledge of the tsunami currents flow and velocities provides insights as for the depth at which tsunamis erode and deposit sediments and whether tsunamis can leave submarine records [11,12]. In this paper, our goal is (1) to test and propose a tsunamigenic source that can potentially affect the eastern coast and (2) examine the flow velocities of the induced water waves and currents from the numerical modeling results and the GIS -approach.

The Tsunami Numerical Modelling
Various types of data indicate that tsunamis are made up of a very long dispersive gravity wave train [13]. Nonlinear dispersive theory is necessary only when examining steep gravity waves, which is not the case in deep water. If the tsunami is generated by a large and shallow earthquake, its initial wavelength and period will be greater. On the other hand, if the tsunami is caused by a landslide, both its initial wavelength and period will be shorter [13]. Moreover, the importance of dispersive effects strongly depends on the extent of the source area (the smaller the source, the stronger the dispersive effects and the ocean depth in the source area. The wavelength of tsunamis and consequently, their period depend essentially on the source mechanism [13]. These water gravity waves are governed by the set of hydrodynamic laws deduced from the theory of the fluid dynamics. The equations of fluid dynamics are representations of the laws of conservation of mass, momentum, and energy as applied to a fluid. In this work, the tsunami modeling has been carried out by using the CLAWPACK package that solves hydrodynamic problems [14]. This package is an open source project and the code is available at http://www.clawpack.org/. The series of algorithms aim at solving linear and nonlinear hyperbolic systems conservation laws (reference). The numerical approach is based on the Godunov's method and the Riemann solver to fix stability problems (nonphysical oscillations) at the interface between grid cells (Finite Volume Method) [15]. Here, the spatially varying flux functions are cell-centered and the method can be applied to gravity waves problems. This package has been adapted to geophysical problems like solving the nonlinear shallow water equations using a depth-averaged of partial differential equations over bottom topography [16]. In addition to classical tsunami modeling where the purpose was to examine the inundation limits and the water heights in the coast, the Geoclaw software was as well applied to examine the velocity flows that resulted from the March 2011 Tohoku Tsunami [12]. Based on this modeling approach and comparing their results with data observations in the Pacific (gauge, video footage), Arcos and Leveque [12] suggested the tsunami current velocity was a more sensitive variable for model validation. The wave propagation analysis and the depth -averaged system of equations are developed as follows [12,15,16]: Where q is the flux (vector of conserved quantities) and = (ℎ, ℎ , 1 2 ℎ 2 + ℎ 2 , ) ; A is the jacobian matrix of the flux function and ∈ ℝ 4 * 4 , h is the water depth, u is the horizontal velocity, g is the gravitational constant and b is the bottom surface elevation.

The GIS Approach
In the scope of this study, open-source tools as OpenStreetMap or Google Earth were used for gaining the necessary information, as well as evaluations of satellite imageries and further Web-tools, in order to create a most current, high spatial resolution GIS integrated, reference database. Satellite data such as Landsat, Sentinel and ASTER images and Digital Elevation Model (DEM) data were downloaded from the USGS EarthExplorer of the USGS / USA and the Global Land Cover Facility, University of Maryland. ArcGIS / ESRI and ENVI / Harris digital image processing software were used for the digital image processing.
Responsible authorities related to disaster management and mitigation are aware that not all potential flooding events caused by earthquake triggered tsunamis, landslide triggered tsunamis or meteo-tsunamis can be covered by precise models due to the complexity of their causes. More frequent extreme weather events due to climate change might increase meteo-tsunamis and storm surge occurrence. Further on, the different causes and their interactions with the coastal morphology and bathymetric conditions have to be considered, depending among other factors on the high energetic wave directions.
Whatever might be the reason of coastal flooding events, it can be stated that generally the flattest and lowest areas are flooded first. Morphometric properties and disposition of coastal areas can influence to a great extent the susceptibility to flooding due to flash floods, storm surge, meteo-or tsunami waves. Thus, the systematic inventory of morphometric properties according to a standardized GIS-approach based on digital elevation model (DEM) data and long-term evaluations of satellite imageries from flooding prone areas contribute considerably to the detection of coastal areas, that are more susceptible to flooding due to their geomorphologic disposition. Based on DEM data morphometric maps were created and terrain parameters were extracted such as shaded relief, aspect and slope degree, minimum and maximum curvature or plan convexity maps using ENVI and ArcGIS software. The integration of different factors in a GIS environment using weighting procedures represented one objective in the GIS application in the frame of this study. The basic pre-requisite for the use of weighting tools of GIS is the determination of weights and rating values representing the relative importance of factors and their categories. The Weighted Overlay Tool in ArcGIS for the flooding susceptibility-approach is aggregating and summarizing causal / preparatory factors influencing the morphologic / morphometric disposition to flooding. Those causal factors were extracted from morphometric maps. Such causal factors are  relatively lowest local height levels, in this case height levels below 10 m,  slope gradients < 10°,  minimum terrain curvature (values = 0 corresponding to flat terrain),  drop raster < 100.000 (The drop raster is calculated using the Hydrology tools of ArcGIS. It is calculated as the difference in z-value divided by the path length between the cell centers, expressed in percentages. For adjacent cells, this is analogous to the percent slope between cells. The result is a map of percent rise in the path of steepest descent from each cell.),  flat areas derived from aspect maps (extraction of values [-1]) and  high flow accumulation values.
These causal factors were aggregated and weighted, in this case starting with equal influence. The resulting weighted overlay map is divided into susceptibility classes. The susceptibility to flooding is classified by values from 0 to 7, whereby 7 is the highest, assumed susceptibility to flooding due to the aggregation of causal morphometric factors. This weighted overlay approach in ArcGIS has been applied in various tsunami prone areas and the results were compared then with tsunami flooding events in the past such as during the 11.March 2011 tsunami in NE Japan and recent catastrophic tsunami events in Indonesia, Morocco and in Chile. A clearly coincidence of the tsunami flooded areas with areas showing the highest susceptibility values on the susceptibility map of each investigation area could be observed in all cases. Therefore, this approach was used as well in the coastal areas of Algeria in order to derive the position of areas prone to flooding due to their geomorphologic properties.

The State of Research
Since 1980, several authors published methods and magnitude scale to emphasize (1) the role of the energy budget from the source mechanism and (2) the total tsunami energy from tsunami waves as a key component for coastal damage. In particular, a tsunami magnitude only related to the energy was proposed [17]. Flooding, water height and run ups observed correspond in fact to the energy budget that is the sum of the potential energy and the kinetic energy. In fact, the common parameters when quantifying what could be the source mechanism at the origin of tsunami waves are the potential energy and its variation and transmission to the water body. Dutykh and Dias [13] proposed a set of equations for tsunami wave energy that are in fact an extension of the simple formula given by Okal and Synolakis [18] to more spatially realistic sea-floor deformations. Dutykh and Dias [13] reported and plotted the time evolution of the kinetic and potential energies for three mechanisms of tsunami generation. For an instantaneous bottom motion (dislocation model) the energy is higher than for an exponential bottom motion (slumping event).
The physics of the earthquakes and the understanding of the rupture processes with measured physical parameters are now investigated since at least the mid-nineteenth century. Thanks to the seismological observatories dataset, several empirical relationships have been proposed to define a rupture zone. Nevertheless, Kanamori and Brodsky [19] well highlighted that despite the approaches widely used by the seismologists are favoured by the development of the seismic inversion data methods, the model that physically represents the failure at the origin of an earthquake is the spontaneous failure.
Until now, it is quite difficult to evaluate the geometry and the location for the rupture area that played during the 21 st -22 nd of August 1856 in the eastern coast of Algeria. The re-appraisal research compiled from historical documents revealed that from a seismological view, the scenario can be approximated but not defined [2]. The macroseismic rupture zone and affected area proposed by the authors raise a critical issue as for (1) the onset of the first main shock, (2) the seismic energy versus the gravitational energy (slide, slumps) radiated during the whole seismic period of the seismic sequence (foreshock, mainshock, aftershocks) and (3) the series of water waves that triggered the flooding of the coast (harbours, inland, and the Soumam valley role on the coastal flooding in Bejaia). Here, the modelling performed aimed at studying the water flow (currents, velocities) that could bring an insight into the coastal risk (seismic, tsunami, environmental).
The sea bottom deformation is simulated based on the theory that estimates the internal displacements and strains due to shear and tensile faults [20]. The source considered is a finite rectangular source. Length and Width of the fault plane are related to the energy released by the rupture area induced by the seismic scenario. The quantification of the size of the rupture area from the seismic moment helps to propose the sudden sea bottom uplift. Geological and Geophysical offshore surveys and investigations in eastern Algeria revealed several features related to destructive past events and active faults offshore and inland in the Algerian margin [21,22]. The deformation identified in soft sediments (quaternary deposits) located in the coast of Jijel are attributed to strong earthquakes (sismoslumps, seismites) [22]. Nevertheless, Benhamouche et al. [22] pointed out that the soft deformation features may have different origins such as seismic shocks, strong turbulence, and high tractive current as well as a collapse of a granular material in slope morphology. Here, the authors highlighted as well that "the action of strong currents and collapse of material may be ruled out because of the absence of slope in the deposit environment of the studied quaternary deposits". The observations gathered by Benhamouche et al. [22] correspond to a spatial distribution of the features that fit the area of maximum intensity of the historical major earthquake of August 1856 that triggered a destructive tsunami. In the coast of Jijel (between 3°75E to 6°10E), "the observed soil settlements (paleo-landslides and liquefaction) is recorded in the landscapes as a 10-m long and 4-m wide paleo-valley where the liquefied sand and the sliding traces are still visibles".
In previous modelling studies, Roger and Hebert [23] considered three subfaults with a strike varying between 80° to 60°N. In fact, four segments were revealed from seismic surveys and swath bathymetry [24,25]. However, only the western segment (5.4764°N, 36.95°E), the center segment (5.736°N, 37.0791°E) and the eastern segment (6.15°N, 37.1784°E) are considered as active structure possibly at the origin of the 1856 events with dips ranging between 75 to 85° [24,25]. The results obtained for six sites along the Algerian coast by Yelles-Chaouche et al. [25] do not exceed 1.5 meters. In particular, the maximum sea wave amplitude between 5°42E and 6°E. Yelles-Chaouche et al. [25] also indicate the height the computed in Bejaia city do not exceed one meter. Their model showed that the maximum energy of the tsunami dissipates more easily towards the North than laterally, indicating that the majority of exposed areas area located southward and northward of the seismic sources. The authors acknowledge that the use of their model do not allow to reach run-up values higher than two meters [25]. For both of the modelling studies [23,25], all considered an echelon fault pattern (NE-SW or E-W directions)at the foot of the djidjelli margin that express the deformation inland and offshore [25].

The Scenario
Marine surveys along the Algerian margin succeeded to map submarine slides from west to east [24]. These campaigns helped to highlight seafloor instabilities related to seismic and neo tectonic activities. In the Bejaia-Jijel region, the continental slope morphology is dominated by widely spaced, rectilinear submarine canyons (see Fig. 3) [26]. Deposit areas and locations where sediment accumulations have been identified and examined to understand the triggering mechanisms for gravity events in the Eastern Algeria margin and off Bejaia in particular [27]. They concluded the sediment accumulations were related to (1) major floods of the Soummam Oued and sedimentary discharge (prodeltaic deposition) controlled by the steepness of the catchment hillsides, (2) past climatic changes and palaeoenvironmental parameters, (3) episodes of strong instabilities of the slopes connected to a clustered-type period of seismic activity [27].
In this work, simplified assumptions are suggested so as to include within the seabottom displacement, all factors that generate water waves. We do not consider a series of subfaults but rather a single segment that covers the whole extent of the affected area during the 1856 seismic sequences in Djijelli-Bejaia. When estimating a co-seismic slip, the procedure uses the Okada theory [20]. The Okada model calculates the deformation due to a rectangular source. It provides analytical solutions for finite faults (shear and tensile) in a half space. Therefore, we suggest in this study to consider a unique fault segment that would fit with the historical observations (earthquakes and sea wide tsunami) but include as well the potential slides' sources for seafloor instabilities. Another criteria for such an approximation is the storm reported by Roger and Hebert [23] in the Western Mediterranean in Mahon Harbour.
The table 1 lists the points of interest selected for the tsunami modelling and the Figure 3 shows the slope bathymetry for the western Mediterranean. Nearby the coast and the four overlapping structures mapped offshore Djidjelli by Yelles et al. [25], the slope varies between 5 to 70° either parallel to the coastline or angles that reflect the location of distinct geomorphological pattern undersea. Hence the slope gradient for the region closed to the points of interest is another argument for slide sources triggering locally tsunami waves.  The geometry of the rupture area tested here presented consists of a rectangular plane of 94.4 km in length by 24.8 km in width and dipping 40° in the SE direction. The plane is oriented with a strike of 75 degree north. In regards with the tectonic setting and previous studies and the regular records of seismic events in the region, the rake angle proposed is 90° degree (pure reverse motion) and the rupture point or focal depth is at 7 km. The epicentre suggested is located at 37.10°N, 5.7°E.
The calculations with the geoclaw software are computed using a regional topographic grid from the NOAA database with a resolution of 1minute (source: http://www.ngdc.gov). It covers an area that extends from 36°N to 42° N and 0°E to 8° E.

The Tsunami Modelling
The bay of Bejaia-Jijel trapped the water waves induced from the seafloor deformation introduced in the modelling (Figure 4). It takes then about 30 minutes for the water flow to propagate from the epicentre area to the Balearic Islands. The tsunami water heights and water waves' velocities are represented in the figures 5 and 6.   The maximum water wave's height computed in this work is obtained for the POI 4 (Point of Interest) with a value of 1 to 1.5 meters in height (Fig.5). Later, the water wave profile is dominated by lower water height waves with larger periods. The POI 3 and POI 5 are marked by maximum wave amplitudes of about 1 meter in height. The lowest waves are computed for POI 2 and POI 6 with a water height of 0.5 meter.
The velocities estimated in this work are oriented eastward for the u component and northward for the v component. The results show the values obtained are higher for the v component (Fig. 6). The lowest velocities are estimated for POI 2 and POI 4 with a maximum value of 8 cm/s for both locations. The maximum of the u component is reached for the second crest for POI2 while it reaches POI 4 for the first crest. The highest velocities are calculated at POI 5 and 6. The u component is about 35 cm/s for POI 5 and 57 cm/s for POI 6. POI 5 is marked by the highest velocity in the southward direction (-82 cm/s).

The Susceptibility Map Flooding
The coastal topography of the Jijel area influences the water current development forming circular patterns. Created by wind blowing across the surface, wave patterns are complex and highly varied. Superimposed wave sets that have been reflected and bent by the coastal features can be observed on satellite imageries. Assuming high energetic waves coming from eastern direction towards Jijel water currents will form that -combined with return flow -will probably have a larger impact than waves coming from western direction due to the coastal morphology.
River mouths are forming an entrance for intruding water waves from the sea. Therefore, even in longer distances from the seaside flooding along the coast-near rivers might be possible.
The flooding of this coastal area, however, is not only related to the high energetic waves such as incident tsunami waves but also to the return flow. In general, the return current flows into the low-elevated areas or into river beds and channels, sometimes partly remaining there for days before discharging into the sea. The weighted overlay approach in ArcGIS can contribute to the detection of those lowlands, that after the return-flow could remain longer time flooded than the environment. Fig.7 presents the results of the weighted overlay based on morphometric factors derived from ASTER digital elevation model data (DEM). The ASTER DEM derived drainage pattern and the weighted overlay results, presented together with the sealed areas of Jijel contribute to a better visualization of urban areas susceptible to flooding (Fig.8) because of their morphologic disposition. By combining these basic susceptibility maps with results of tsunami modelling of specific cases early warning and disaster management could be improved, provided a quick data flow in case of emergency. Figure 7. Weighted overlay of causal, morphometric factors such as the lowest height level (< 10 m), slope degree < 10°, curvature = 0, aspect = (-1), drop raster < 100.000. The dark-blue areas are related to those areas with the highest susceptibility to flooding due to their morphometric disposition

The Flooding Risk
The tsunami modelled water waves all correspond to water flow simulated nearby the entrance of rivers, at river mouths. Hence additional considerations regarding the sediment discharge and the canyons geomorphologic features must be as well examined as part of the flooding risk in case of a potential destructive earthquake in the region. The results obtained from the earthquake scenario offshore the Jijel bay indicate the waves' energy mostly propagates in the north and south direction. Our study is in good agreement with the work carried out by Yelles et al. [25].
Nearby the harbour of Jijel (POI 4), the water height simulated is the highest while the velocity is the lowest. As this location is the epicentral area, the ground shaking is as well the highest. Waves trapped in the harbour are sources not only for the coastal flooding but for collisions between the vessels anchored. In the Djendjen harbour (POI 5), the water waves' velocities simulated are the highest and are oriented in the southward direction. Meanwhile, the remote sensing approach revealed the region is marked by a high susceptibility flooding. Finally, the Soummam valley constitutes an important asset for the region as for the economy of the Bejaia region (agricultural, etc…). The results obtained near Bejaia (POI 2 and POI 3) emphasize the necessity to add local tsunami waves from nearby canyons. The energy transported from the main source (Jijel epicentre) mainly dissipates in the south direction. The oued Agrium would be largely flooded if we add a local slumping event from the Bejaia Canyon.

Conclusions
The combination of the remote sensing and GIS approach and the tsunami modelling helped to define the potential flooded and/or affected locations in case of a destructive earthquake in Eastern Algeria. Satellite images contribute to a better understanding of the interactions between the coastal morphology and the development of water currents, of course depending on meteorological and tidal conditions at the acquisition time of the images. The weighted overlay approach in ArcGIS contributes to the detection of those coastal areas in the Jijel and Bejaja regions that are prone to flooding due to their morphometric disposition such as river mouths and flat lowlands, thus, providing a basic flooding susceptibility map. As this GIS approach has shown promising results in comparable areas [28][29][30][31], it might be useful here as well.
Local tsunamis must be considered when studying main destructive events like the1856 damaging earthquake in Jijel. The occurrences of slides inland (Kherrata region) and the canyons offshore near the Soumam valley and the Agrium river are sources for slumping events. Deep canyons might focus incoming high energetic waves from the sea-side resulting in potential higher run-up, inundation and damages in the affected coast lines. The water flow along the coast might induce material damage and debris floating nearby the harbours. The water waves' energy and the areas marked by a high susceptibility flooding mostly at rivers mouth are sources for pollution and related-environmental risk issues. Today, the trade market consists mainly on industrial activities with oil tankers and refineries in Bejaia and in Djendjen regions.