**Geodiversity Evaluation and Water Resources in the Sesia Val Grande UNESCO Geopark (Italy)**

#### **Luigi Perotti \*, Gilda Carraro, Marco Giardino, Domenico Antonio De Luca and Manuela Lasagna**

Earth Sciences Department, University of Torino, 10125 Turin, Italy; gilda.carraro@edu.unito.it (G.C.); marco.giardino@unito.it (M.G.); domenico.deluca@unito.it (D.A.D.L.); manuela.lasagna@unito.it (M.L.)

**\*** Correspondence: luigi.perotti@unito.it

Received: 16 September 2019; Accepted: 4 October 2019; Published: 9 October 2019

**Abstract:** This paper aims at systemizing knowledge related to geodiversity assessment for water resources and its evaluation. The novel aspect connected to geodiversity of this paper is the analysis of the components of hydrological system, both at the superficial and underground level, in the territory of the Sesia Val Grande United Nations educational, scientific, and cultural organization (UNESCO) Global Geopark (Northwest Italy). More specifically, the research establishes a conceptual model and a specific procedure for the evaluation of geodiversity connected to water resources on a regional scale, by means of a qualitative-quantitative geographic information system (GIS) process, renamed here as hydro-geodiversity assessment. For these purposes, a targeted ecosystem approach is applied to consider the assets of the Geopark territory that has been derived from the interaction between water and other components of geodiversity, i.e., the hydro-geosystemic services. The element selection and processing operations led to the identification of areas characterized by greater values of hydrological geodiversity, in which the link between surface and underground hydrodynamics became closer and intense. The single geodiversity factor maps that were obtained from partial data aggregations were added together in map algebra operations, then subjected to weighing to formulate the hydro-geodiversity map of the Sesia Val Grande UNESCO Global Geopark. The results of the present study strengthen the strategic management of geological, geomorphological, and hydrological heritages of the study area by identifying different landscapes and local peculiarities determined by mutual influences between geology and hydrological dynamics.

**Keywords:** water resources; geodiversity assessment; geosystem services; geoheritage; hydro-geodiversity; Sesia Val Grande UNESCO Global Geopark

#### **1. Introduction**

The term "geodiversity" has no intrinsic value; its importance relies on the quality of the relationships built between the systems or spheres of which it is composed, i.e., the Earth system science (atmosphere, lithosphere, hydrosphere, and biosphere; Figure 1) [1] and those specifically addressed to describe surface processes, landforms, and materials, such as pedosphere and anthroposphere. Interactions between these different spheres constitute the variety of geological and geomorphological phenomena and landscapes [2] to which human beings attribute several values. Therefore, we agree with Sharples [3] in considering geodiversity as "the quality we are trying to conserve," and geoconservation as "the endeavor to conserve it".

**Figure 1.** Conceptual scheme of the relationships between the main spheres considered by the Earth systems science [1]. Asterisk indicates the area of relationships with pedosphere and anthroposphere.

According to Gray's (2013) definition, geodiversity is not just a matter of different Earth features [1] but also of their assemblages, structures, systems, and contribution to landscapes. The complexity of geodiversity is a challenge for its assessment. A system of norms and modes of action, as well as a model of interaction of all the related variables, has to be created in order to achieve significant geodiversity assessment, in term of both biotic and abiotic ecosystem services. Such a comprehensive assessment can offer relevant contributions to both geoconservation and sustainable use of georesources [4].

This work is focused on systematizing the relationships between geodiversity and water resources by analyzing parts of the hydrological cycle that interact with geological features, satisfy essential needs, and for allow biological and human evolution. The first research question we want to address is: What value should be attributed to water resources in the definition of geodiversity?

Despite the broad interests and various focuses of contemporary geodiversity methods and tools, the related literature considers hydrology as a fundamental component of geodiversity assessment [1,5]. However, relationships between geology and water resources have been frequently analyzed in a sectorial way:


Although the methods used in the two aforementioned studies have opposite aims and results (quantitative versus qualitative), their main tendency is to perform partial analyses of the interaction between water and geodiversity, i.e., by focusing on specific geological or geomorphological phenomena (in these last two cases: landforms and river dynamics), thus neglecting other relevant elements of the hydrological system.

On the other hand, some authors [9,10] posed a greater effort in a systematic attempt to define the structure of the hydrogeological interactions, including factors and variables useful for creating a map of water resources. In the study by Arajuo and Pereira, such a factor map represents a fundamental part of the final geodiversity map of the State of Cearà, in Brazil [9].

Despite the undoubted inclusion of hydrological processes in the geodiversity equation adopted by Brazilian and Polish assessments [9,10], the real crux of the matter is that, to date, there is no uniform choice of essential variables to be considered for classification of hydrological elements relevant for geodiversity assessment. This is particularly relevant for regional geodiversity studies, such as large-scale surveys and assessments, where a theoretical framework has to be carefully discussed for achieving successful applications [11].

To overcome the problem, we sought to define and test a set of relevant variables for a qualitative-quantitative assessment of hydrological-geodiversity in the Alpine territory of the Sesia Val Grande UNESCO Global Geopark (UGG). In a preliminary phase, it was necessary to devise a relational model that could provide support for the collection and management of information to be processed. This important and preliminary research operation included the definition of a conceptual model for geodiversity assessment in relation to water resources, based on the specific territorial context, and the adopted regional scale of analysis.

#### **2. Study Area**

The Sesia Val Grande UGG [12] study area covers about 2000 km<sup>2</sup> and is located in the Piemonte Region (NW Italy) (Figure 2). The area includes 106 municipalities across 4 provinces: Verbania, Vercelli, Novara, and Biella. Elevation of the territory varies from 190 m a.s.l. at the lower alpine piedmont area to 4554 m a.s.l. at the top of Monte Rosa (Pennine Alps), the second highest massif of the European Alps. Indeed, the area is mainly mountainous, including high plains and large floodplains, as well as a portion of the Maggiore Lake.

**Figure 2.** The study area of the Sesia Val Grande United Nations educational, scientific, and cultural organization (UNESCO) Global Geopark (UGG).

The study area comprises three important hydrographic sub-basins of the Po drainage basin (Figure 2): those of the Sesia river (3079 sqkm area, 138 km length), Toce river (1785 sqkm area, 57 km length), and Ticino river (6033 sqkm area, 284 km length) [13]. Many factors and competing morphodynamic processes contributed to the shaping of the mountain relief and the valleys system, i.e., the litho-structural and tectonic conditions posed by alpine orogeny and the morphoclimatic variations, such as the Pleistocene glacial/interglacial phases and later Holocene stages. These "regional" factors are related to long-term processes, giving the basic shape of the Alpine mountains and valleys, then followed by important Holocene "local" morphogenetic processes of fluvial/torrential, as well as of gravitational origin. Currently, the dominant geomorphic agency in the valleys of Sesia Val Grande Geopark is the fluvial-torrential one, which is accompanied by consistent gravitational instabilities, where the slopes are steeper. [14]. A large portion of the territory shows fluvial/torrential landforms along steep slopes (16% or more), such as deep river incision, mainly in bedrock, with abundant debris deposits, often forming debris fans up to km-size. On the other hand, in valley floor and high plain areas, rivers created terraces at the valley sides, and multiple to single channels had a tendency to meander at the valley mouth [15].

At higher elevations, glaciers have been the most important morphogenetic agent during Holocene. Indeed, the upper areas of the high valleys are dominated by glacial and periglacial processes [16]. Despite the ongoing climate warming and the predominant southern exposure of the slope of the Monte Rosa Massif, seven glaciers are still present [17], whose hydro-geosystem value is undeniable: they constitute the sources of the Sesia river, as well as of a beautiful landscape, even if they are extremely sensitive and endangered by climate change.

Concerning the tectonic setting, Figure 3 shows the main structural units and geological complexes, in which the geopark is located. Units of the Southern Alps are aligned along a northeast-southwest direction, juxtaposed to the Austroalpine units along the Insubric Line (towards W) and to the lower Pennidic units along the Sempione-Centovalli Line (towards N). In turn, the Austroalpine is in tectonic contact by means of complex polyphasic deformation with the Pennidic domain, herein represented mainly by oceanic units [18].

The lithological geodiversity constrains the Geopark's hydrological features. The whole area (Figure 3) is mainly characterized by crystalline basement units, i.e., lithologies whose permeability depends on fractures density. A little portion of conglomerates, limestones, and marbles outcrops in the Lower Sesia Valley, while fluvial and fluvioglacial deposits characterize the valley floors and the piedmont. These last units represent the reloading areas of the high productivity aquifers present in the territory.

The various geo-lithological complexes, combined with their structural settings, contribute to the large geological diversity of the area and, at the same time, define highly diversified hydrographic network and rich hydrogeological structures [19], as follows:


6. Alluvial deposits: these sediments are located in the valleys bottom and in the piedmont. Due to their porosity they show prevalent permeability from high to medium and contain a shallow unconfined aquifer in connection with surface rivers.

**Figure 3.** Map of the geological and structural units of the Sesia Val Grande UNESCO Geopark. (modified from Brak et al. 2010) [18].

#### **3. Materials and Methods**

In the following paragraph, we describe the adopted methodological approach, the selected parameters, and the methodology chosen for carrying out the evaluation process that led to the final map of hydro-geodiversity of the Sesia Val Grande UGG.

#### *3.1. Methodology*

Before defining the input data for hydro-geodiversity assessment, it is necessary to define the conceptual structure of the assessment. In this specific evaluation, the intention is to consider values and services that the territory and community derived from the abiotic components in connection with the dynamics of the water and the formation of aquifers.

Gray's (2013) model of geosystemic services identifies five different categories of geoservices [2] and was adapted in a hydrogeological overview and applied to the local context of Sesia Val Grande UNESCO Global Geoparks area (UGGp). This kind of approach, usually part of qualitative methods to assess geodiversity [20], is human-centered. Our analyses identified hydro-geosystem services, which represented hydrogeological features capable of offering a range of specific services and goods.

The conceptual process guided the assessment of hydro-geosystem services and is described in Figure 4. Starting from the analyses of relationships between geodiversity and water, we proposed a framework for hydro-geodiversity. Since it represents the part of geodiversity concerning the hydrosphere, it includes hydrogeological phenomena that interacts with geolithological features, the component of geomorphological landscape, and the way in which human societies manage them.

According to the conceptual scheme of geosystem services by Gray (2013), interlinked categories have been found and the intensity and types of relationships are described in Figure 4:


These interactions determine a range of hydro-geoservices that forge the structure of the hydro-geosystem, which corresponds to a certain degree of hydro-geodiversity. The conceptual definition of the hydro-geosystem services in the territory under study was essential to understand and define the input data to consider the hydro-geodiversity map of Sesia Val Grande UNESCO Geopark.

**Figure 4.** Hydrogeoservices from Gray (2013) [2].

#### *3.2. Hydro-Geodiversity Assessment*

Once analyzed, the characteristics of the territorial context defined the conceptual setup of the research, it was possible to proceed with the definition of the parameters of the hydro-geodiversity assessment. The specific parameters are described in detail in Table 1. The operational purpose is to identify areas characterized by high hydro-geodiversity using a qualitative-quantitative evaluation technique.

**Table 1.** Chosen parameters for hydro-geodiversity assessment in the case study area (modified from Zwolinski, Najwer, and Giardino, 2018) [20].


The hydro-geodiversity assessment procedure is typically quantitative, based on the pioneering work of Serrano and Ruiz [11]. It is the kind of approach based on the construction of map algebra indexes and techniques, using geographic information system (GIS) software to process information. To achieve this practical purpose, we performed a GIS analysis by using ArcGis 10.5 software (developed by ESRI Redlands, USA) on a complete set of georeferenced spatial data. On the basis of the available data and the survey scale, a hydro-geodiversity equation was established, whose variables corresponded to the factor maps that were added together using weighing techniques in the map algebra phase.

Due to the large study area, the chosen scale for the evaluation was 1:100,000. Indeed, several relevant features for geodiversity assessment are at a nominal scale of 1:10,000. In order to obtain a final representation that was compatible with the finale factor maps and consistent with the chosen scale of representation (1:100,000), we did a semi-automatic data generalization.

Based on the selected scale of analyses, a geodatabase was constructed by collecting the public data provided by regional and territorial agencies, as described in the data source field described in Table 1.

The main steps for the hydrological geodiversity assessment were:


Data selected in the initial phase have undergone changes due to scale compatibility or type of input data. An iterative process was adopted [27]. Available data were collected and analyzed, and we observed the results and determined which data could be used and how.

Four main factors were chosen for the evaluation of hydro-geodiversity:


These factors represent the variables of the hydro-geodiversity (HGD) equation, which can be summarized as:

$$\text{HGD} = \text{tP} + \text{tLU} + \text{SWD} + \text{MR} \tag{1}$$

From the vector data and the expert classification of the elements, the data has been transformed into a raster format. This allows us to assign a value to each identified class and to add the obtained images with a final resolution of 25 m.

Figure 5 briefly describes the methodology adopted in the evaluation assessment, as well as highlighting the relationships between the fur factor maps in creating the hydro-geodiversity map.

**Figure 5.** Flow diagram for the creation of the final map of hydro-geodiversity in the Sesia Val Grande UNESCO Global Geopark.

Once the four factor maps have been obtained, the next step includes the processing of partial maps via attributing weights through map algebra. By varying the weight of the individual partial maps in the map algebra phase, it was possible to evaluate and compare the results of the weights assigned to the individual factor maps. To find the best weight proportions used to obtain a result that identifies sufficiently homogeneous areas of hydrological geodiversity, the AHP (Analytic Hierarchy Process) method [28] was used, which is a multi-criteria decision support technique.

Finally, seven hypotheses of hydro-geodiversity assessment were formulated. Only one was chosen as a representative for the Sesia Val Grande UGG hydro-geodiversity. It was then reclassified into three distinct classes using the natural break method.

Consequently, an interpretation of the results were made, leading to the definition, within the geopark, of a certain number of landscapes characterized by high hydro-geodiversity.

#### **4. Results**

#### *4.1. Factor Map of Total Permeability*

The geological lithology [21] was divided by type and degree of permeability. More specifically, the basement rocks, quaternary, and pre-quaternary deposits were classified with values from 1 to 5 based on the hypothetic degree of permeability, which is directly related to the constitution and the productivity of aquifers that hide from lithological formations (Table 2). Both the deposits and rock basements were classified by hypothetical permeability, which underpins their predisposition to constitute aquifers.

**Table 2.** Lithology of rock basement and pre-quaternary deposits in the Sesia Val Grande UGG.


The highest values were assigned to quaternary and pre-quaternary gravelly and sandy deposits. Intermediate values were instead assigned to lithologies such as marble, limestone, and sandstone with mixed deposits or debris flow. Ultimately, very low values were assigned to most coherent lithologies or glacial deposits (Table 3).


**Table 3.** Classification of quaternary deposits in the Sesia ValGrande UGG.

A particular procedure was followed with regard to lacustrine and marsh deposits, since these deposits contain clay and are characterized by low permeability. However, theses deposits constitute layers of protection for aquifers and represent superficial environments dominated by the water dynamics that constitute extremely precious biotopes (e.g., peat bogs, marshes, ponds). Thus, classifying these deposits on the basis of permeability means giving them low values; this is not in line with the objective of the present study, which seeks to enhance the centrality of the water element in its interaction with geodiversity. This led to the choice of not considering these deposits in the present classification, but to evaluate them with high value in the factor map of total land use. To complete the factor map of total permeability, the state of fracturing of the substrate and influencing the degree of

permeability of the rock was considered. Indeed, the lithologies of the basement, such as crystalline rocks, have a wide permeability range depending on the level of incidence of tectonic structures. In the study area faults, fault systems and ductile shear systems were distinguished [21]. Moreover, structures responsible for ductile, ductile-fragile, or brittle-ductile deformation were distinguished, because in ductile-fragile areas in the tectonic action produce a greater incidence of fracturing.

In order to create a map of the fracturing index, several hypotheses of classification were advanced. The final decision was to create manually areas of fracturing relevance, attributing more value to fragile deformations as compared to ductile ones. These areas were added to the map of deposits and rock basement permeability, in order to obtain a map of total permeability (Figure 6).

In the map, it is possible to observe how rocks with low degree of lithological permeability assumed maximum values, i.e., the Insubric Line, the Cossato-Mergozzo-Brissago, and the Pogallo Lines.

#### *4.2. Factor Map of Land Use*

Land use is an important factor to consider in the equation of hydro-geodiversity because it explains the human impact on natural environments. In the hydro-geodiversity assessment, it is important to highlight all variables that seal or pollute the ground. Thus, the factor map of land use collects all types of land uses identified by the corine land cover satellite tracking system, which integrates them with wetlands, marsh areas, and lakes (Figure 7).

The identified elements were classified based on their possible effect on ground permeability, as well as the possibility to create underground reserves and water resource pollution.

Because of this, the factor map of land use collects all types of land uses identified by the corine land cover satellite tracking system. The elements identified are then classified based on the possible conditioning of the ground or riverbed permeability, as well as the quality of the underground reserves.

Table 4 summarizes the considered variables, classifying them from the lowest hydro-geodiversity value (1) to the highest (5). Regarding the lake data, provided by the regional landscape plan (PPR), only the water elements with a surface greater than 100 m × 100 m were selected. This measure corresponds to the minimum "cartographic" resolution at the 1:100,000 scale.


**Table 4.** Expert classification of Land Use types in Sesia Val Grande UGGp.

Considering the number of areal landslides, the map of land use has been integrated with the map of slope instability. Landslide phenomena were considered on the basis of their density; for reasons of scale adaptation, only landslides with a surface area greater than 100 m × 100 m were selected. The slope instability index was obtained by analyzing the density of the area landslides converted to point format. The result of the kernel density analysis was then reclassified into three classes using the natural breaks method (values from 0 to 2).

#### *4.3. Factor Map of Springs and Wells Density*

The location of springs and wells were mapped separately. Then, in order to obtain a final representation compatible with the other factor maps and consistent with the chosen scale of representation, areas with a higher density of springs and wells were identified (Figure 8).

Natural springs and wells were subjected to a kernel density with a radius of 1000 m [29]. The raster file obtained was then reclassified into four classes (0, 3, 4, 5) with manual classification, turning the areas characterized by low density into a value of 0. For this classification, high values were used to stress the importance of these factors.

**Figure 6.** Factor map of lithological permeability integrated with fracturing index. (Esri ArcGis 10.5, [30].

**Figure 7.** Factor map of land use integrated with landslides density index [30].

**Figure 8.** Factor map of springs and wells density [30].

#### *4.4. Factor Map of Morphogenetic Relevance*

The factor map of morphogenetic relevance (Figure 9) is used to consider the predominant geomorphological factors that characterize the study area, as well as the dynamics and genetic processes that are the basis of morphological conformation. The territory of the Sesia Val Grande UGG was therefore divided into three areas of morphogenetic relevance: glacial, fluvial, and gravitational. Geomorphological elements taken into consideration are glaciers and glacial cirques for glacial relevance, hydrographic network, alluvial fan, lakes for fluvial relevance, and areal landslides for gravitational relevance. These areas were expertly classified with values from 3 to 5. For the areas dominated by the glacial processes was given the maximum value (5). Glacial modelling is indeed a central factor in the geomorphology of alpine areas.

**Figure 9.** Factor map of morphogenetic relevance [30].

Glacial cirques often host lakes and mountain pastures, which constitute rare habitats and areas used for anthropic purposes for grazing. At the same time, glaciers constitute a reserve for drinkable water and important element of river flow regulation.

The river and lake elements constitute fluvial relevance, which represents a high value of hydrological geodiversity. Lakes and rivers are reservoirs of water. Moreover, rivers can be used for energy production (e.g., dams, hydroelectric power plants) and provide aggregates (e.g., gravel, sand, silt, peat) for various uses. Lastly, the lowest value (3) was attributed to the areas of gravitational relevance. This value has a moderate to high estimation since the landslide processes are firmly interrelated with the water dynamics. They can, in fact, activate and be activated by water processes.

Once the final factor maps were obtained (Figures 6–9), they were added together in a map algebra operation, their sums weighted with GIS, and put through AHP.

This criterion allowed us to create measures that judged consistency, derived priorities between criteria that allowed for comparisons, and established a hierarchy of priorities among the elements. The weighing method adopted made it possible to elaborate many hydro-geo-assessment solutions

(from HG\_A\_1 to HG\_A\_7 in Table 5). Each time, a greater weight was assigned to one of the four factors, as illustrated in Table 5.


**Table 5.** Hydro-geodiversity assessment solutions. The underline values indicate the group of factors that has more weight in the map algebra process.

The last two solutions, HG\_A\_6 and HG\_A\_7, are the result of a reasoning that considers a more rigorous approach, which was adopted in the present work and the objectives set. A greater weight was used in the land use factor map, containing the elements interacting with the human dimension. While the HG\_A\_6 shows an increasing value of factors, the HG\_A\_7 shows that lower weights are equal. In the HG\_A\_7 solution, land use weight is equal to 40% of the total weight. Moreover, the springs and wells density is equal to 25%. In fact, since these factors are connected to human activities, these results indicates systems of provisioning and pumping of water resources, as well as a strong point for monitoring the quality and quantity of water in deep and shallow aquifers. To the natural abiotic factors, like permeability and areas of morphogenetic relevance, a weight of 17% was assigned to both levels. It was therefore considered that the HG\_A\_7 solution was the best solution for the hydro-geodiversity assessment.

In order to obtain more homogeneous areas, the raster file was reclassified into three classes using the natural breaks method (presented in Figure 10 with specific areas and hydrogeosites).

#### **5. Discussions and Conclusions**

In this study, nine peculiar areas in the Sesia Val Grande Geopark were identified on the base of the prevailing landscape and its propensity to develop a sustainable relationship between man and the hydro-geosystemic services (Figure 10):


Once identified, hydro-geodiversity areas were compared with the geosite location in order to validate the correspondence between them and the areas of hydro-geodiversity, as well as to verify their representativeness and to test the functioning of the qualitative-quantitative procedure previously applied.

**Figure 10.** The final solution of the hydro-geodiversity assessment reclassified into three classes with natural breaks in the area definition and hydrogeosites [30].

Geosites from the list of that identified for candidacy to the UNESCO program of the Sesia Val Grande Geopark [31], as well as those extrapolated from the ISPRA national inventory of geosites [32] were selected.

Only geosites with a significance in terms of hydro-geodiversity were chosen. In particular, 25 geosites in this selection were plotted on the final map. All geosites (with the exception of 1) fell into areas with high hydro-geodiversity.

Despite this good correspondence, it should be noted that in some areas there are more than one geosite (areas 4, 8) and in other areas, geosites are classified with a high hydro-geodiversity. Occasionally, we noted the total lack of geosites (areas 1, 2, 6,7). If we were analyze the features of geosites, we would note that the categories of representativeness expressed are fluvial, glacial, gravitational, karst, and lacustrine.

Geosites that represent and test the aquifer dynamics and the relationship between the geological structure and the concentration of springs are missing. However, this aspect is considered fundamental in the hydro-geodiversity classification. This is not surprising and is in line with the tendency to underestimate hidden geosystem features, such as underground processes. This is also demonstrated by recent results of a systematic literature review [33] that show how goods and services derived from the subsurface are underrepresented in the contemporary literature on ecosystem services.

Based on the results and the comparison with the current state of geoconservation of the study area, area 7 (landscape of springs) and area 9 (landscape of the deep aquifers of the Upper Po Plain) are the most important areas in terms of hydro-geosystemic services, as they are directly related to the withdrawal and consumption of water (e.g. drinking water, for agriculture, for breeding). They are also the areas in which human impact is deeper and where there are no instances of hydrogeological protection sufficient for a good preservation. Therefore, more studies and insights about these issues is needed.

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

**Funding:** This study was supported and funded by the Compagnia San Paolo bank foudation in partnership with University of Torino, under the funding program 2016; the project "GeoDIVE—From rocks to stones, from landforms to landscapes" was funded with grant #CSTO169034.

**Conflicts of Interest:** The authors declare no conflicts 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* **Impact of the October 2018 Storm Vaia on Coastal Boulders in the Northern Adriatic Sea**

#### **Sara Biolchi 1, Cléa Denamiel 2, Stefano Devoto 1,\*, Tvrtko Korbar 3, Vanja Macovaz 4, Giovanni Scicchitano 5, Ivica Vilibi´c <sup>2</sup> and Stefano Furlani <sup>1</sup>**


Received: 15 August 2019; Accepted: 22 October 2019; Published: 25 October 2019

**Abstract:** Boulder detachment from the seafloor and subsequent transport and accumulation along rocky coasts is a complex geomorphological process that requires a deep understanding of submarine and onshore environments. This process is especially interesting in semi-enclosed shallow basins characterized by extreme storms, but without a significant tsunami record. Moreover, the response of boulder deposits located close to the coast to severe storms remains, in terms of accurate displacement measurement, limited due to the need to acquire long-term data such as ongoing monitoring datasets and repeated field surveys. We present a multidisciplinary study that includes inland and submarine surveys carried out to monitor and accurately quantify the recent displacement of coastal boulders accumulated on the southernmost coast of the Premantura (Kamenjak) Promontory (Croatia, northern Adriatic Sea). We identified recent boulder movements using unmanned aerial vehicle digital photogrammetry (UAV-DP). Fourteen boulders were moved by the waves generated by a severe storm, named Vaia, which occurred on 29 October 2018. This storm struck Northeast Italy and the Istrian coasts with its full force. We have reproduced the storm-generated waves using unstructured wave model Simulating WAves Nearshore (SWAN), with a significant wave height of 6.2 m in front of the boulder deposit area. These simulated waves are considered to have a return period of 20 to 30 years. In addition to the aerial survey, an underwater photogrammetric survey was carried out in order to create a three-dimensional (3D) model of the seabed and identify the submarine landforms associated with boulder detachment. The survey highlighted that most of the holes can be considered potholes, while only one detachment shape was identified. The latter is not related to storm Vaia, but to a previous storm. Two boulders are lying on the seabed and the underwater surveys highlighted that these boulders may be beached during future storms. Thus, this is an interesting example of active erosion of the rocky coast in a geologically, geomorphologically, and oceanologically predisposed locality.

**Keywords:** rocky coast; extreme waves; active erosion; geohazard; Croatia

#### **1. Introduction**

This study investigates the movement of the boulders located in the southern sector of the Premantura (Kamenjak) Promontory (Croatia, southern Istria Peninsula) after an extremely severe storm named "Vaia" which struck during 29 October 2018 [1]. Terrestrial, aerial and underwater surveys were performed to analyze the possible geomorphological effects of the storm on the boulder field in terms of movements of boulders, as well as their disappearance and the appearance of new ones.

In 2018, between Saturday, October 27 and the early hours of Tuesday, October 30, large sectors of the Italian and the northern Adriatic coasts were hit by one of the most intense, complex, and damaging storms in recent decades as the result of the passage of an exceptionally strong cyclone named Vaia. Monday, October 29, was characterized by violent gusts of Sirocco, (a southerly wind blowing across the Adriatic Sea), storm surges, and extraordinarily high tides in the northernmost Adriatic Sea combined with severe rainfall, particularly in the northeastern Alps. The storm caused 16 casualties in an area from Trentino (northern Italy) to Campania (southern Italy) and severe damage totaling more than two billion euros. Tens of thousands of customers were still without electricity two days after the event, especially in the Trentino, Veneto, and Friuli areas (Northeast Italy).

We investigated the relationship between this extreme Sirocco storm and an extensive boulder accumulation located at the southernmost tip of the Istrian Peninsula, in Premantura, Croatia [2]. This study provides the results of two years of monitoring using unmanned aerial vehicle (UAV) based multitemporal photogrammetric image comparison and the outputs of underwater environmental monitoring using a three-dimensional (3D) model obtained through photogrammetric analysis of underwater pictures. Field activities were carried out, assisted by a UAV survey, approximately two weeks after storm Vaia. These surveys revealed that the boulder positions had changed with respect to the previous field investigations carried out in November 2017. For example, the isolated boulder, K8, widely described by Biolchi et al. [2] had been moved towards the northwest and rotated.

The transport of rocky boulders in the Mediterranean and more generally along oceanic coasts has been widely studied. It has been documented that waves associated with stormy events and major tsunami events can generate a force strong enough to detach boulders from the ground and transport them along the shore. Despite boulder motion being mostly attributed to tsunamogenic activities in the literature, several authors have recently proven that the impact of storm waves is, in fact, the primary mechanism for boulder detachment and transport [3–9]. Moreover, at higher latitudes, Orviku et al. [10] explained how the decomposition of sea ice and its drifting onshore can transport and accumulate boulders over 1 m in diameter, especially during the spring after ice melting.

Marriner et al. [9] analyzed tsunami and storm data, contained in the EM-DAT (Emergency Events Database) for the 1900 to 2015 period, and found that up to 90% of tsunami attributions of high-energy events in the Mediterranean coastal records should be reconsidered. Cox et al. [11,12] and Cox [13] reported on strong evidence that even some mega boulders in the open ocean setting (Northern Atlantic) were displaced along with the cliff platforms by storm waves.

This study presents a multidisciplinary approach based on the hypothesis that it was storm Vaia that impacted the coastal boulders. This approach, which integrates inland and submarine surveys, as well as numerical modeling aims to:


The focus of this work is to study the relationship between a previously documented storm and coastal boulder movements, highlighting that a severe storm is sufficient to move or initiate boulder transport both from the sea bottom and in a subaerial environment. The work can be included in the scientific debate concerning the explanation of coastal boulder detachment and movement in terms of "storm vs. tsunami". A clear relationship between a severe storm that occurred during 2014 and the emplacement of a large boulder in the southern Istria peninsula in Premantura was provided by [2]. This study, performed in the same area after a different storm, confirms that sea wave energies are capable of moving or detaching new boulders from the submarine floor.

#### *1.1. Study Area*

The studied boulder accumulation is located on the rocky coast of the Premantura Promontory (Kamenjak Nature Park) in the southernmost sector of the Istrian peninsula (Croatia, northern Adriatic Sea). The promontory is formed by a succession of stratified Late Cretaceous carbonate rocks, dipping gently towards the east [2] (Figure 1).

Its southeastern-most tip is named "Jugo Promontory", after the powerful southeasterly wind (jugo in Croatian) and consists of a stratified, shallow-marine Cretaceous limestone succession. The succession is characterized by alternations of thin-bedded (10–30 cm), fine-grained peloidal packstones, and thick-bedded (50-150 cm) mudstones and wackestones containing algal oncoids and rare rudist shells, typical of the lower part of the Gornji Humac Formation, which dates back to the Turonian [14]. The Kamenjak Nature Park is characterized by elevations that are up to 50 m above sea level, rounded bays, pocket beaches, and small islands. The south and east coasts of the park are characterized by gentle slopes, which generally follow the dip angles of the bedrock strata.

At the studied site (southern part of Jugo Promontory), the washed ("white") coastal zone is up to 70 m wide. The zone is the narrowest (40 m) at the highest elevation region, in the central part of the study site. The bed dip direction and the dipping angle on the southwestern part is 88/12, and boulders are grouped along the washed zone, situated 50 to 70 m from the sea, along the border with the vegetated zone, figured in Biolchi et al. [2]. There are a few solitary boulders that are considered the youngest. More than 950 mostly meter-sized boulders (volumes from 0.385 to 11.440 m3) were recognized on the promontory during a former survey [2], with an estimated total volume of ~2000 m3.

The submerged part of the coast is of similar morphology, and gently dipping limestone ramps continue under the sea. Cliffs with heights of a few meters developed locally (out of the studied localities) and are related to massive rock masses within the succession of the limestone rocks and to small local faults. Thus, the maximal runup during extreme storms is defined by the height of the breaking waves and the gentle dip angle of the coastal ramp and is marked by the boundary of the vegetated zone along the gentle coast that is around 70 m away from the sea [2].

The Adriatic Sea is a semi-enclosed basin in the northern Mediterranean, with the following dominant winds: Sirocco (jugo in Croatian) blowing from the southeast, Bora (bura in Croatian) blowing from the northeast, and Libeccio blowing from the southwest (Figure 1c) [15]. The Sirocco has the longest fetch, and therefore generates the highest waves in the northern Adriatic Sea [16]. Significant and maximum wave heights of 5.3 m and 10.8 m, respectively, were measured over a roughly 10-year interval (1978–1986) about 50 km southwest of the Premantura Promontory. The other two winds generate lower waves in the area under investigation. The Sirocco may also generate extreme storm surges and sea levels in the northern Adriatic Sea [17], particularly between November and February, which can lead to the flooding of coastal regions. Together with the Adriatic seiche [18] and tides, extreme sea levels can reach up to 1.5 m above mean sea level [19,20]. On top of the tides and the storm surges, tsunami-like waves of meteorological origin, meteotsunamis, may further raise the sea level along the open coastline of Istria by a few tens of centimeters [21,22].

Although they have been found to impact coastal boulders in various parts of the world's oceans, no significant seismic tsunamis have been reported in the northern Adriatic Sea over the last two hundred years [23,24]. The worst-case hazard scenarios provide for a maximum tsunami height no higher than 20 cm [25].

**Figure 1.** (**a**) Location of the Istrian peninsula, (**b**) location of the Premantura study area, (**c**) Windrose, (**d**) geological map of the South part of the Premantura (Kamenjak) Promontory, and (**e**) oblique view of the study area (Jugo promontory and Cape Kamenjak) where the K8 boulder and clusters of boulders spread along rocky coast are clearly visible.

#### *1.2. Storm Vaia*

At the time of this study, storm Vaia had not yet been described in the literature and the only available information came from reports provided by the Italian Meteorological Society [1] and the numerical results (Figure 2) obtained during this study with the Adriatic Sea and Coast (AdriSC) modeling suite presented in more details in Section 3.3.

Atmospherically speaking, the Vaia depression developed on Saturday, 27 October 2018, within an extensive low-pressure trough stretching from the Baltic to the western sector of the Mediterranean Sea.

This depression held its position at sea between the Gulf of Lion, the Balearics, and Sardinia until the morning of Monday, October 29. This first phase of the storm was marked by humid Libeccio air currents coming from between the south and southwest, accompanied by intense rainfall events in the northern Apennines and mountainous areas from northwestern Italy to the Carnic Alps (northeastern Italy). By midnight on October 28, more than 300 mm of rain had already fallen in many areas from the Prealps above the city of Brescia to the mountains of Friuli, with a noticeable increase in river flows. After a pause of a few hours, the storm's second phase, fueled by the arrival of the first major cold front of the 2018 season in the Alps, developed on the morning of Monday, October 29. During the course of the day storm Vaia underwent rapid deepening (about 17 hPa in 18 h), almost classifiable as "explosive cyclogenesis" (the threshold for which is considered a pressure decrease of 24 hPa in 24 h), with the center of the depression moving from the west of Corsica to the northeast during the afternoon and then to northwestern Italy in the evening (red circles in Figure 2). This deepening of the storm was accompanied by the following: (1) a brutal reinforcement of the winds across Italy, (2) a major strengthening of the Sirocco, and (3) the development of violent self-regenerating storm cells in the area of Sardinia, and the Tyrrhenian and Ligurian Seas. In the afternoon and evening of the October 29th a southeasterly windstorm reached the Adriatic basin while a heavy rain began to fall again on the already saturated soils of the Alps and Prealps but also in Northwest Italy due to its more southeasterly airflow.

The Vaia depression will be remembered not for the rainfall, but for the violence of the Sirocco that blew between morning and afternoon of Monday, October 29, when it swung round to the Libeccio in the evening, starting with Italy's western sea areas (i.e., measured wind speeds reached 102 km/h at Rome Ciampino airport, 119 km/h at Genoa airport; 128 km/h in Lugano in southern Switzerland, 128 km/h on the Valles Pass in the Dolomites, 148 km/h at Capo Carbonara in southeastern Sardinia, 155 km/h at the Colle di Cadibona near Savona in NW Italy, 171 km/h at both La Spezia and Follonica in Tuscany, and 200 km/h on Monte Rest in the Carnic Prealps in northeastern Italy). The violent southerly gusts of wind, which persisted for many hours along considerable lengths of coastline, raised devastating storm surges, particularly in Liguria, with severe damage to the coastal roads and railways, buildings, and tourist facilities with dozens of boats destroyed in ports. On the evening of Monday, October 29, a buoy belonging to Liguria's Environmental Protection Agency lying offshore from Capo Mele recorded a maximum wave height of 10.3 m with a very long period of 11 s, an indicator of highly destructive wave power along the coast. Also noteworthy was the episode of high water in Venice and along the northern Adriatic coast, including the Istrian peninsula. In particular, the tide gauge at Punta della Salute (at one end of Venice's Grand Canal) measured a maximum of 156 cm at 14:10 on 29 October, a value exceeded by only three other events in the historical series since 1872: 1 February 1986 (158 cm), 22 December 1979 (166 cm), and the infamous 4 November 1966 (194 cm).

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

The investigations on the boulder accumulation located along the southern coast of the Premantura (Kamenjak) Promontory have been developed through a multidisciplinary approach, which envisaged inland surveys, submarine investigations, and wave modeling.

#### *2.1. Inland Surveys*

Field activities have been performed since 2017 and include traditional geomorphological and geological investigations and UAV surveys. Geomorphological and geological activities have included extensive fieldwork and are reported in [2]. The above-cited paper includes a list of boulders and their axis lengths.

UAVs are widely used in various fields of geoscience [26–28] and offer major benefits through their ability to provide high-resolution photographic images from reduced flight times [29]. The digital photogrammetry (DP) technique enabled the reconstruction of high-resolution orthophotographs and a digital elevation model (DEM) of a large sector of the Premantura area. This reconstruction was done through the processing of 135 images collected by a DJI Phantom drone TM (DJI, Nanshan District, Schenzen, China) in 2017 (flight altitudes between 20 and 60 m).

These images were processed using Agisoft Photoscan Professional software version 1.4.0 (Agisoft, St. Petersburg, Russia) and the results were processed in 2018 before storm Vaia in a GIS environment using QGIS version 2.18 Las Palmas (QGIS Development Team, QGIS Geographic Information System. Open Source Geospatial Foundation Project).

One of the main outputs of the photogrammetric processing was a detailed map of the position of 950 boulders spread along the Jugo Promontory, as shown in [2].

In order to carry out a multitemporal photogrammetric survey devoted to recognizing possible movements of blocks triggered by storm Vaia and other future storm events, we performed three drone campaigns in November 2018 (two weeks after storm Vaia) and 2019, using a quadcopter DJI Spark droneTM (DJI, Nanshan District, Schenzen, China), equipped with a 12MP camera. DJIFlightplanner softwareTM (AeroScientific, Blackwood, Australia) and Litchi softwareTM (VC Technology, London, UK) assisted in the choice of 2019 flight plans and permitted the operator to easily set the altitude, radius, number of waypoints, speed, and directions.

We performed two flights on 15 November 2018, and six flights at different altitudes on 30 April 2019, and 14 June 2019. The DJI Spark droneTM executed automatically and autonomously the 2019 flight plans using the Litchi softwareTM installed on the remote controller device.

Table 1 lists which sector of the study area was investigated in each flight and the main characteristics of the UAV surveys.


**Table 1.** Main characteristics of unmanned aerial vehicle (UAV) flights carried out in 2018 and 2019.

In order to define the ground control points (GCPs) for the production of georeferenced aerial images, the flights carried out in 2019 were accompanied by a differential global navigation satellite system (DGNSS) survey [30]. The GCPs were mainly located on the west side of the boulder accumulation where we had noticed the recent movement of the isolated boulder named K8 in [2], and, we assumed, further possible movements of a few solitary boulders because their distribution was closer to the coastline. The master GCP was located in an area that may not be reached by any boulders, whereas 6 DGNSS bolts accompanied by black-yellow photogrammetric targets were spread in positions clearly identifiable from aerial photos taken in previous UAV campaigns. These DGNNS points were crucial for the quantification of possible boulder motion, with a maximum error of 15 cm calculated on the maximum discrepancy between the reconstructed coordinates and the one measured at the edge of the area.

The images were processed using 3DF Zephyr Aerial softwareTM version 4.007 (3Dflow, Verona, Italy) and detailed orthophotographs of 2018 and 2019 were computed. The best quality photos were taken on the later flights of 2018 and 2019 and for this reason we devoted our attention to these orthophotographs to analyze the western and central sector of the study area. These photos were imported into the QGIS software, which permitted the evaluation of possible boulder rotation and translation transport.

The distance between two known and measurable points on the ground, such as the pothole described in [2] and the old position of boulder K8 (both clearly visible on aerial photos of 2018 and 2019), was measured using a QGIS tool. This measurement was compared with the real distance measured directly on the ground on 30 April 2018. The difference of approximately 10 cm validates the 15 cm accuracy of the model computed by the 3DF Zephyr Aerial softwareTM.

#### *2.2. Submarine Investigations*

In order to reconstruct the 3D model of the seabed, underwater images were collected via a snorkel survey in June 2019. Photogrammetric procedures are a common and rapid technique to acquire metric measurements in a submarine environment using the combination of field operations and post-processing [31,32]. Submerged photogrammetry requires a precise preparation and setting of the images to reduce water refraction and distortion. Underwater photography is vulnerable to the different refraction coefficient of the marine water, which produces a reduction of the field of view (FOV) and a complete change in the optical parameter of the lens [33]. To reduce the changes in focal length, we used a dome glass. Conversely, this increases the distortion of the system significantly. In this work, a "quick and dirty" procedure was applied to evaluate the limitations of a low-cost approach using GoPro action cameras together with the capability of acquiring good quality data snorkeling from the surface. Action cameras are usually equipped with wide-angle lenses and the option of selecting a narrow FOV via software adjustment, losing a large percentage of the sensor resolution and resulting in a low-quality set of photos. We operated in full resolution to manage the wide FOV resulting in strong distortion and a huge color cast at the picture boundaries.

The 3D model was assembled in Agisoft Metashape TM (Version 1.5.0) starting from a selection of 916 of the 1200 images produced by a GoPro Hero 6 BlackTM (GoPro, San Mateo, California, USA) set in a six inch waterproof plastic case.

In order to remove the color cast and the large amount of suspension present at the picture corners, a mask was applied to each photo. This strategy removed almost one-third of the pixels belonging to the area with the strongest distortion and permitted a precise alignment.

To better cover the submarine area, an S-shape route was followed with a time-lapse camera set to 0.5 s. Therefore, the overlap of the underwater images was higher than in the aerial photogrammetry. This procedure avoided inaccurate orientation of the cameras at different depths of the seabed [34].

#### *2.3. Wave Modeling*

The numerical modeling was carried out using the Adriatic Sea and Coast (AdriSC) modeling suite, that was jointly developed within the ADIOS [35] and MESSI [36] projects with the aim of accurately representing the processes driving the Adriatic's atmospheric and oceanic circulation at different temporal and spatial scales. This modeling suite consists of: (1) a basic module that deals with the coupled ocean and atmospheric general circulation and (2) a nearshore module that provides high-resolution fields during extreme events.

The basic module is based on the Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) modeling system [37]. The module is built around the Model Coupling Toolkit (MCT), which exchanges data fields and dynamically couples the Weather Research and Forecasting (WRF) atmospheric model, the Regional Ocean Modeling System (ROMS), and the Simulating WAves Nearshore (SWAN) model. The basic module was set up with (1) two different nested grids of 15 km and 3 km resolution used in the WRF model and covering the central Mediterranean area and the Adriatic-Ionian region, respectively, and (2) two different nested grids of 3 km and 1 km resolution used in the ROMS and SWAN models and covering the Adriatic-Ionian region (the same grid as the WRF model) and the Adriatic Sea, respectively.

The nearshore module additionally accounts for the nearshore bathymetry changes and combines the WRF model with the fully coupled ADCIRC-SWAN unstructured off-line model [38]. In this module, the hourly results from the WRF 3 km grid obtained with the basic module were downscaled to a 1.5 km

grid covering the Adriatic Sea, while the hourly results from the 1 km ROMS grid and the 10 min results from the 1 km SWAN grid are used to force the unstructured mesh of the ADCIRC-SWAN model.

The AdriSC modeling suite was installed and fully tested on the European Centre for Medium-Range Weather Forecast (ECMWF) high-performance computing facilities. More details on the AdriSC modeling suite setup can be found in [39].

In this study, in order to reproduce the storm which took place in the Adriatic Sea on 29 October 2018, the SWAN model was set up in both modules to be coupled with the oceanic and atmospheric models (i.e., WRF, ROMS, and ADCIRC). In addition, the computation of the bottom stress of the ocean models (respectively, ROMS and ADCIRC) was updated in order to consider the spatial distribution of the sediment grain size at the bottom of the Adriatic Sea, extracted from the Adriatic Seabed database [40], and the wave effects. To reproduce the storm as accurately as possible, the basic module was set up to run for three days between 27 October and 30 October 2018, with initial conditions and boundary forcing provided by (1) the 6 hourly ERA-Interim reanalysis fields [41], (2) the daily analysis MEDSEA-Ocean fields [42], and (3) the hourly MEDSEA-Wave fields [43]. The nearshore module, forced by the results of the basic module, was set up to run between midday on 28 and 30 October 2018.

In this study, only the wave results from the unstructured SWAN model of the nearshore module (hereafter referred to as AdriSC unSWAN) were analyzed.

Due to the lack of precise bathymetry data and resolution in the nearshore area where the waves are breaking (i.e., surf zone) and the known limitations of the physics used in the SWAN model (e.g., the parametrization of the wave breaking), the wave heights modeled by the AdriSC unSWAN module during the event were extracted off the surf zone (ideally in deep water) and the Sunamura and Horikawa equation [44] was applied to evaluate the wave height at the breaking point, Equation (1):

$$\frac{H\_b}{H\_0} = \left(\tan\beta\right)^{0.2} \* \left(\frac{H\_0}{L\_0}\right)^{-0.25} \tag{1}$$

where *Hb* is the breaking wave height; *H*<sup>0</sup> and *L*<sup>0</sup> are the wave height and the wave length off the surf zone (ideally in deep water), respectively; and β is the slope of the seabed near the coast (i.e., in the area where the waves are expecting to break). Finally, the breaking wave height obtained was compared with the hydrodynamic model proposed by Nandasena et al. [45] which is extensively used in the literature to define the storm wave heights capable of detaching and transporting coastal boulders.

#### **3. Results**

#### *3.1. Onshore Surveys*

The area where the limestone boulders are distributed is located along the southern coastal section of the Jugo Promontory, mostly between the sea and the vegetated zone. The southwestern coast of Cape Kamenjak is directly exposed to the Sirocco-induced waves. The dip direction and dip angles of the limestone beds are 88◦ and 12◦, respectively. An indistinct joint system has developed along the bed strike (i.e., generally running North-South), whereas a distinct open fractures system generally strikes east and west (with a mean dip direction and dipping angle measure of 350◦/85◦) as per meter-scale distances. Thus, quadrangular limestone fragments are formed by the fracture network together with layer planes [2]. The seven tonne boulder (K8), here renamed boulder #1, is characterized by its unusual orange surface coloring due to karst weathering. Its elevation is 2 m above sea level. Before the storm, boulder #1 was located at 27 m distance from the coastline and was oriented with the longer axis facing the main wave direction. Given its isolated position and intense orange color, it is clearly distinguishable in aerial and terrestrial images. The occurrence of very fresh subrecent biogenic carbonate encrustations on its southern and upper surfaces, mainly produced by coralline algae and serpulids, as well as by more fragile barnacle shells, attests to its marine origin. Its deposition has been reconstructed by [2] and is the result of a severe Sirocco storm between January and February 2014.

During the field activities carried out on 15 November 2018, just two weeks after storm Vaia storm, the shift of boulder #1 was the most evident and visible impact of the storm (Figure 3).

**Figure 3.** (**a**) Oblique view of boulder #1 and (**b**) the red rectangle indicates its past position on the limestone pavement.

As a result of the print that was left by boulder #1 after four years lying on the limestone pavement and the multitemporal comparison of UAV-derived orthophotos, its displacement was measured as being approximately 3 m towards the northwest and with a counterclockwise rotation of 18◦.

The UAV surveys in 2018 and 2019 were carried out to identify and quantify any other boulder movements. The comparison between the 2017 and 2018 UAV orthophotos showed that 14 of them had changed their position or had been rotated and revealed the appearance of a new small boulder. Figure 4 highlights the boulders that had moved, plotted on the 2018 UAV orthomosaic, and compared to their 2017 positions.

**Figure 4.** Spatial distribution of 14 boulders affected by the 29 October storm Vaia. The base image was taken during the November 2018 UAV flight. The star indicates the presence onshore of a new limestone boulder detached and moved from the seafloor during the storm. No movements were recorded after 15 November 2018, and 30 April 2019.

All the boulder that had moved were measured using field measurements (sizes, position, and distance from the coast), UAV-obtained images, and 3D models. These 3D models are crucial for a desk-based geomorphological analysis, such as boulder recognition, and size measurement. Figure 5 displays a view from Flight #2, performed on 30 April 2019.

**Figure 5.** A view of the three-dimensional (3D) model obtained during Flight #2 performed on 30 April 2019. The red arrows show the shift of boulders #1 and #9 and the position of the new boulder #10.

Table 2 lists the boulders that were affected by the storm including: size, quantification of translational movement (with an accuracy of 1 m), and rotation characterization. Moreover, an estimation of the wave height required to move a boulder in a subaerial scenario according to the Nandasena model is reported.

Boulders #2, #3, #4, #7, and #8 are trapped in two different clusters populated by tens of blocks and have mainly been rotated. Boulders #5 and #6 are close to the coastline and their motions were interrupted by a persistent east–west oriented joint that crosses the limestone promontory. The discontinuity is clearly visible in Figure 4. This movement caused a rotation of Boulder #5 and its fragmentation. Boulder #11 was deeply fragmented by severe waves and probably by collisions with other fragments, whereas boulders #12, #13, and #14 have mainly been rotated and fragmented. This behavior is related to the presence of other boulders and the abovementioned persistent discontinuity that inhibited their movements. The same scenario took place with boulders #5 and #6. The western part of Figure 4 shows well the severe impact of the last storm. The yellow star indicates the presence of a new limestone boulder (boulder #10), detached from the seafloor and transported by the waves. North of this new boulder (Figure 5), Boulder #9 was moved about 5 m and is now embedded within a niche.


**Table 2.** List of boulders moved or rotated by the storm.

Conversely, no movement or rotation of the boulders was detected between mid-November 2018 and 30 April 2019, either in the field or after the comparison of the orthophotographs. During this period, no severe storm affected the Istrian coasts.

#### *3.2. Submerged Landforms from Photogrammetry*

The surface of the seabed exhibits a tabular shape, due to the limestone bedding (Figure 6a), and is pitted by several holes, ranging from a few centimeters to about 2 m in size. Most of them are potholes (Figure 6c), occasionally filled by rounded cobbles.

**Figure 6.** (**a**) Limestone beds below the mean sea level, (**b**) a detachment scar, (**c**) the largest pothole in the submerged area, and (**d**) isolated boulder.

An oblong-shaped hole, about 2 m in length and 1 m in width, is less colonized by marine life than the surrounding seabed and lies offshore from the limestone bed hosting boulder #1.

The 3D model shows a sloping seabed, reaching its maximum depth at about 50 m from the coastline. Landward, the seabed is cut by a persistent joint parallel to the coastline (the same as described in the previous paragraph), while a limestone bed borders the eastern side of the study area (Figure 7).

**Figure 7.** (**a**) A 3D model of the seafloor at the Premantura site, (**b**) the roughly S-shape route followed during the snorkel survey, (**c**) a 3D sketch of a sector of the submerged area, and (**d**) a submerged profile. The submerged landforms on the seafloor are clearly visible in the model, particularly the limestone beds, potholes, and other rounded abrasional landforms.

#### *3.3. Wave Modeling of the 29 October 2018, Adriatic Sea Storm*

In order to evaluate the skills of the AdriSC modeling suite to reproduce storm Vaia, the unSWAN wave model results were, first, compared with both state-of-the-art wave model running on the entire Mediterranean Sea and wave measurements from buoys located along the Croatian coastline.

The qualitative comparison across the entire Adriatic Sea of the AdriSC unSWAN significant wave height and peak period with the operational MEDSEA-Wave results (Figure 8) during the peak of the storm (i.e., 29 October 2018, at around 20:00) shows that the AdriSC unSWAN model is capable of capturing the wave conditions during the storm. However, it is fundamental to notice that the peak period in the northern Adriatic Sea is higher for the AdriSC unSWAN model than for the MEDSEA-Wave results. This might be due to the difference in wind force, the physics of the models or the difference in resolution, which leads to a better representation of the nearshore geomorphology of the Adriatic Sea by the AdriSC unSWAN model than the 4 km MEDSEA-Wave model.

**Figure 8.** (**a**) Significant wave height using MEDSEA-Wave, (**b**) significant wave height using AdriSC unSWAN, (**c**) peak period using MEDSEA-Wave, and (**d**) peak period using AdriSC unSWAN. All four models are related to 29 October 2018, at 20:00.

In Table 3, the AdriSC unSWAN maximum significant wave height and its associated peak period during the 29 October 2018 storm are compared at four different nearshore locations (Rovinj, Split, Dubrovnik, and Ploˇce in Figure 8) with the measurements obtained by the Croatian Hydrographic Institute (Hrvatski hidrografski institut, HHI) and reported in [46].

**Table 3.** Comparison of the wave parameters (significant wave, height, and peak period) measured by the Croatian Hydrographic Institute at four different locations (Rovinj, Split, Dubrovnik and Ploˇce) together with the AdriSC unSWAN model results during the peak of the storm on 29 October 2018.


In general, the AdriSC unSWAN model shows some skill in reproducing the measurements at the four nearshore locations along the Croatian coastline, however, the unSWAN model overestimates the peak period in Rovinj by nearly 2 s and underestimates the significant wave height in Dubrovnik by about 0.5 m (Table 3). The overestimation of the peak period of the AdriSC unSWAN model in Rovinj is consistent with the results obtained in Figure 8, and thus the model is likely to have generally overestimated the peak period in the northern Adriatic Sea during storm Vaia on 29 October 2018.

Given the qualitative and quantitative comparisons performed, the AdriSC unSWAN model is thought to have reasonably reproduced the storm on 29 October 2018, and can be used to assess the impact of the waves on the boulders of the Premantura area. Near Premantura (Figure 9a), the AdriSC unSWAN model has a resolution ranging from 100 m at the coastline to 1 km further offshore, and a bathymetry that captures the main geomorphological features of the seabed, however, the islands of Fenera, Ceja, and Bodulas are too small to be included in the mesh and are each represented with one point with a depth of 0 (zero) m (Figure 9b).

**Figure 9.** AdriSC unSWAN mesh structure (panel **a**) and bathymetry (panel **b**) in the vicinity of the Premantura area.

Results from the AdriSC ADCIRC and unSWAN models during the storm on 29 October 2018, for the study area, are presented in Figure 10. The modeled storm produces sea surface elevations up to 0.5 m in most of the area (Figure 10a), except in Medulin Bay, where it reaches 1 m due to seiche activity. Currents (Figure 10b) reach their maximum strength in coastal regions. They are particularly strong at the tip of the Premantura Promontory, where they were modeled with velocities up to 1.7 m/s. Such strong currents are presumably induced by strong waves inclined to the coastline. Significant wave heights reached 7 m (Figure 10c), with a peak period of 12 s (Figure 10d) and a wavelength of 105 m (Figure 10e) in the wider Premantura area. These wave parameters, particularly the significant wave height, are greater than any recorded measurement in the northern Adriatic Sea and, following an analysis performed 50 km west from the Premantura Promontory [47], may be considered as occurring once every 20 to 30 years.

**Figure 10.** Spatial variability of (1) ADCIRC maximum sea surface height (panel **a**) and maximum current speed (panel **b**) and (2) unSWAN maximum significant wave height (panel **c**), maximum peak period (panel **d**), and maximum mean wave length (panel **e**) obtained in the vicinity of the Premantura Promontory during the storm of 29 October 2018.

To investigate the boulder motions on the Premantura Promontory (as the nearshore bathymetry of the model may not be accurate due to a lack of both resolution and accuracy of the bathymetry data), the wave parameters were extracted using a point off Cape Kamenjak, at a depth of about 29 m (Figure 10) and the wave height at the breaking point was calculated using the Sunamura and Horikawa equation [44]. On 29 October at 20:34 the maximum significant wave height reached 6.2 m with an associated maximum wave height of 10.8 m, peak period of 11.1 s, and mean wave length of 103.3 m. Assuming the same bottom slope as in [2], the wave height at the breaking point for the storm on 29 October 2018, is 13.6 m, which is higher than the theoretical wave heights needed to move or

transport the boulders presented in Table 2, except for Boulder #7 that requires wave heights of more than 19 m according to Nandasena [45].

#### **4. Discussion**

The comparison of orthomosaics obtained from UAV images, repeated annually, can be a very useful tool for monitoring the variation of the boulder positions. This methodology has been successfully applied in both North Atlantic [11–13,48] and Mediterranean boulder sites [49]. UAV-derived orthophotographs and digital surface models (DSMs) can provide excellent data and information on coastal boulder patterns. Orthophotographs allow for the mapping of a axes and b axes, including their orientation, whereas precise values for c axes and boulder volume can be taken from the DSM [50]. In recent years, after an initial decade when the analyses were mostly dedicated to explaining boulder detachment and transport mechanisms and distinguishing between those of tsunamigenic and those of storm origin, boulder studies are now more oriented towards the observation of boulder movements and reassemblage after exceptional storm events [13,50]. In addition, Hastewell et al. [51] have proposed an innovative technique for boulder movement reconstruction, using radio frequency Identification (RFID) tagging and DGNNS technology. A single RFID tag is inserted inside a boulder and can be activated with an electromagnetic signal emitted by a pole antenna moved by an operator also equipped with a backpack reader and a handheld computer.

On the Premantura Promontory, where the origin of the boulders has been ascribed to storm events, totally excluding those of tsunami origin [2], and where boulders analysis began only recently (end of 2016), a proper monitoring network including orthophotograph comparison, both aerial and submarine, geomorphological observation, instrumental monitoring techniques, and numerical modeling, is most definitely needed as severe storms are forecast to increase in the near future [52].

Following the long-term monitoring experience of coastal landslides affecting the northwestern part of Malta [53–55], a network of GNSS benchmarks has also been installed in Premantura. At the same time, the submarine 3D model, presented here, could be the starting point for monitoring the effects of future storms in the underwater environment, that represents the source area for the quarrying and detachment of future boulders. Finally, every day since early 2019, the AdriSC operational component has been providing a 48 h forecast of high-resolution wave parameters, which can be used to detect in advance whether or not a potential storm can move boulders.

By comparing UAV-derived orthophotographs, we detected the movement of 13 boulders and the emplacement of a new one after a severe storm. The western part of the boulder field suffered the intensity of the last storm most, whereas the central and eastern part did not show any evidence of boulder movements. This lack of movement is related to the structural setting of the promontory, where, in the central and eastern sectors the boulders (probably older) are located at higher elevations. Thus, active erosion occurs in the western part of the Jugo Promontory (Cape Kamenjak).

By observing boulder movements (Figure 4), the maximum inundation flow of the storm waves capable of transporting boulders was estimated to about 50 m of distance from the coast (that has been measured considering wave directions and orientation of the limestone ramps). The displacement followed the azimuth of the limestone beds strikes that, in turn, have conditioned the detachment of rocky material forming scarps throughout the past. These scarps, together with flat, gently inclined limestone bedding planes (pavements), have acted as ramps for boulder transport. Isolated boulders, such as boulder #1 (K8 in [2]), are those that have been subject to the largest movements, while the boulders gathered in clusters were only rotated or toppled, but remained in position, trapped by other boulders or the abovementioned scarps, fractures, and faults. The latter, in particular, have been enlarged over time by coastal karst and marine weathering, as well as by the removal of small portions of rocky material. Some boulders exhibited signs of fragmentation, due to collisions with other boulders or against the rocky scarps.

The submerged scenario, accurately obtained using digital photogrammetric reconstructions, is a high energy environment, with fresh detachment scarps and rounded or sub-rounded potholes, where decimetric cobbles and boulders lie and move, causing the abrasion of the rocky seafloor. The topographical setting of the seafloor reflects the coastal geomorphology and topography, where limestone beds alternate with scarps and faults. With respect to the previous underwater geomorphological surveying, any recent detachment scar has been noted, while a clear detachment niche has been reconstructed in detail. The 3D submarine model obtained, reconstructed using digital photogrammetric analysis of underwater pictures, will be the starting point for future monitoring of the submerged environment, both before and after exceptional storms, that may affect the Adriatic Sea even more given the future climate predictions [56].

After three years of studies and surveys, we can state that storm wave erosion at the Premantura site is active. The appearance of new boulders (especially the seven tonne boulder #1 during the 2014 storm) attests to a dynamic process of sea erosion that is continuously changing the coastal setting. As already suggested by [2], the erosion process was more intense, in the past, when the coast was initially brought into contact with the sea, and the removal of rocky material due to the repetitive action of waves began. Radiocarbon dates provided by [2] suggests that this process has been active over several centuries.

Regarding the numerical modeling, the AdriSC modeling suite presented some skills to reproduce storm Vaia and the unSWAN wave model provided high-resolution wave parameters in the vicinity of Premantura with an unstructured mesh capturing the global geomorphology of the area. The models also allowed for the quantification of the coastal dynamics, particularly the total elevation, including the wave set up and the wave-induced currents, which may significantly impact the wave field when reaching values as high as those modeled off the Premantura Promontory. However, the wave results indicate that the storm produced an estimated maximum wave height of 13.6 m at the breaking point near the area of interest, which was not enough to move Boulder #7. The hydrodynamic equations were used to have an idea about the wave height. The model of Nandasena is mostly used and generally presents a good fit with the measured waves, such as in our study. Of course, the model has its limits and does not offer a precise wave height. The parameters of the equations, such as the lift coefficient, which was calculated using laboratory experiments by [57] for certain particle size and density conditions, are probably not applicable to the local geological, topographical, and climatic conditions where boulders move.

This study shows the following: (1) the need to accurately survey the nearshore bathymetry in the vicinity of the Premantura Promontory (e.g., via shipborne single-beam, multi-beam, side-scan sonar sensors or airborne laser scanning bathymetric surveys) and (2) the limitation of using the SWAN model which, contrarily to the XBeach model [58] — which is developed for wave propagation, sediment transport and morphological changes of the nearshore area, beaches, dunes and backbarriers during storms, cannot correctly reproduce the dynamics in the swash zone (i.e., the land-ocean boundary) where some of the boulders are located. This study provides the first attempt to model the wave conditions responsible for the boulder motion during storm events in the Premantura Promontory, but more modeling efforts and better bathymetry data will be required in the future to truly quantify the wave effects on the boulder dynamics, particularly in the swash zone as described in [59,60]. In addition, the return period of the wave height modeled during storm Vaia, in the northern Adriatic Sea, is about 20 to 30 years. Given that the most extreme storms could, although not significantly, tend to increase slightly in the northern Adriatic Sea in future climate scenarios [61], this might, together with the increase in mean sea level, increase the erosion of rocky coasts. Additionally, other areas, such as the lowlands along the coasts of the northern Adriatic Sea (some of which are subsiding) [56], or the coastal cities with substantial cultural heritage [62], may even be more endangered by the combination of mean sea level and wave activity. Therefore, both these factors should be included in any assessment of the vulnerability of the northern Adriatic Sea to climate change.

Finally, taking into account all the geomorphological observations, the aerial orthophotographic comparison, the wave data in a time range of 1 year (from November 2017 to November 2018), and the wave results from the AdriSC model during the storm, the recent boulder movements and rearrangements can be definitively linked to the 29 October, storm Vaia.

#### **5. Conclusions**

The Premantura Promontory represents a unique example of an extensive coastal boulder accumulation in the northern Adriatic Sea triggered and rearranged by storm waves, as reported by [2] and confirmed in the present study.

The occurrence of this kind of boulder deposit depends on various factors (i.e., the discontinuity network, marine wave direction, and coastal exposure) and for these reasons the integration of underwater and onshore surveys is crucial in the understanding of the processes involved in boulder transport.

The integration of geomorphological surveys and multitemporal UAV-DP permitted the identification of 13 boulders that were moved or rotated and the emplacement of a new boulder during the severe storm that hit the northern Adriatic coasts on 29 October 2018, i.e., storm Vaia.

The western part of the promontory suffered the greatest impacts of the waves as a result of the structural settings of the limestone layers that acted as favorable pavements for the movements of the boulders.

The wave model confirmed that the storm waves had the potential to move these boulders during the storm Vaia event.

Boulder movement on the Premantura Promontory is, thus, linked to the frequency of exceptional storms capable of generating extreme wave heights. Following previous studies carried out in the Adriatic Sea, the return period of the waves modeled during storm Vaia is roughly 20 to 30 years, however, in the last four years, at least two major extreme storms have affected the northern Adriatic Sea (2014 and 2018), causing breaking waves exceeding 10 m height and capable of causing boulder movements.

We demonstrated that the tsunamogenic origin for coastal boulders movements is not required, especially where and when severe waves repeatedly hit the coast. In fact, in the study area the topographical and the geomorphological setting, together with the exposition of southern Istria towards the most severe winds and waves of the Adriatic Sea, favoured boulder detachment and accumulation.

Ongoing research should mainly focus on the improvements of the DP procedure in submarine environments and on the analysis of outputs of DGNNS surveys. Survey zero was carried out on 30 April 2019, when six GNSS benchmarks were installed on six boulders already displaced by storm Vaia. Moreover, the underwater 3D model could be the basis for future comparisons of submerged environments after future storms.

**Author Contributions:** S.D. and G.S. operated the UAV surveys; S.B., S.F., and V.M. carried out diving operations; the geomorphological and geological surveys were carried out by S.B., S.D., S.F., and T.K.; C.D. and I.V. performed the wave modeling. Original idea, S.F., S.D. and S.D.; Conceptualization, all authors; Writing—original draft, all authors; Supervision, S.F.; Writing-review & editing, S.D. and S.B.

**Funding:** This work was carried out within the framework of the Geoswim Project (Resp. Stefano Furlani, University of Trieste), the MOPP–Medflood INQUA-CMP 1603 Project, and the Croatian Science Foundation projects GEOSEKVA (Grant IP-2016-06-1854) and ADIOS (Grant IP-2016-06-1955).

**Acknowledgments:** The authors are grateful to Falck family for their financial support and thank Matteo Mantovani, Linley Hastewell, Andrea Morosetti, and Valeria Vaccher for field surveys and monitoring support. Acknowledgement is also made for the support of the ECMWF staff, as well as for the ECMWF's computing and archive facilities used in this research. The ECMWF ERA Interim dataset is available at [63].

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **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* **Coastal Vulnerability Assessment along the North-Eastern Sector of Gozo Island (Malta, Mediterranean Sea)**

#### **Angela Rizzo 1, Vittoria Vandelli 2,\*, George Buhagiar 3, Anton S. Micallef 4,5 and Mauro Soldati <sup>2</sup>**


Received: 2 April 2020; Accepted: 12 May 2020; Published: 15 May 2020

**Abstract:** The coastal landscape of theMaltese Islands is the result of long-term evolution, influenced by tectonics, geomorphological processes, and sea level oscillations. Due to their geological setting, the islands are particularly prone to marine-related and gravity-induced processes, exacerbated by climate change. This study aligns different concepts into a relatively concise and expedient methodology for overall coastal vulnerability assessment, taking the NE sector of Gozo Island as a test case. Geomorphological investigation, integrated with analysis of marine geophysical data, enabled characterization of coastal dynamics, identifying this stretch of coast as being potentially hazardous. The study area features a high economic value derived from tourist and mining activities and natural protected areas, that altogether not only make coastal vulnerability a major concern but also the task of assessing it complex. Before introducing the methodology proposed for overall vulnerability assessment, an in-depth revision of the vulnerability concept is provided. The evaluation was carried out by using a set of key indicators related to local land use, anthropic and natural assets, economic activities, and social issues. Results show that the most critical areas are located east of Marsalforn including Ramla Bay, an important tourist attraction hosting the largest sandy beach in Gozo. The method combines physical exposure and social vulnerability into an overall index. It proves to be cost effective in data management and processing and is suitable for the identification and assessment of overall vulnerability of coastal areas to consequences of climate- and marine-related processes, such as coastal erosion, landslides and sea level rise.

**Keywords:** coastal morphodynamics; climate change; vulnerability index; Gozo; Malta

#### **1. Introduction**

Coastal zones are host to very dynamic and complex environmental systems, subject to the direct and indirect influence of a number of factors that have contributed to their evolution over time. The present-day landscape of coastal areas is the combined result of interactions between natural agents that include physical factors inherent in the system and external climatic and marine forces [1,2]. Human activity also plays an important role in shaping coastal dynamics, often exerting additional pressures that may dominate over natural processes [3]. During the past decades, the rate of human occupation along littoral areas has risen significantly [4] and currently, in Europe, about 86 million

people are estimated to live within 10 km from the coastline [5]. Furthermore, approximately one-third of the Mediterranean population is concentrated along its coastal regions and around 120 million inhabitants are concentrated in coastal hydrological basins located in the southern region of the Mediterranean Sea [6]. Coastal areas are the transitional zone between the aquatic and the terrestrial ecosystems and they have an environmental intrinsic value on account of their high level of biological diversity, which supports the provision of several ecosystem services essential for human well-being [7,8]. In view of these considerations, the sustainable conservation of coastal areas is a worldwide issue and coastal vulnerability evaluation and risk assessment are of paramount importance for integrated coastal management [9].

Coastal hazards, including marine-related and gravity-induced processes such as landslides, coastal erosion, storm water runoff and coastal flooding, have different impacts on coastal values, due to differing local geomorphological and anthropogenic settings. The impact of these processes is expected to increase with climate change. The assessment reports of the Intergovernmental Panel on Climate Change [10–13] emphatically forecast increase in both the frequency and the severity of extreme weather/climate events, combined with sea level rise, which altogether will undoubtedly impact coastal areas more adversely.

The Sendai Framework for Disaster Risk Reduction 2015–2030 [14], whose implementation is overseen by the United Nations Office for Disaster Risk Reduction (UNDRR), has recognized risk assessment as an important preliminary action on the basis of its relevance internationally. This points towards the prioritisation of the role of "understanding disaster risk" and "enhancing disaster preparedness" and promotes the use of "foresight" and "scenarios" for improving the level of preparedness for existing, emerging and new types of risk. At the same time, the identification of suitable climate adaptation strategies, which also account for future scenarios, is essential for increasing the resilience of coastal areas and ensuring the long-term conservation of coastal natural services and anthropic activities [15].

At the European level, risk assessment and mapping have been included in a number of legislative instruments, such as the Floods Directive [16], which requires member states to provide maps and indications for risk management, prescribing specific requirements for climate change impact evaluation.

As highlighted in the Special Report on "managing the risks of extreme events and disasters to advance climate change adaptation" published by the Intergovernmental Panel on Climate Change (IPCC) [17], vulnerability has been a pivotal concept for disaster risk reduction since the 1970s. Initial studies in this field, carried out by Baird et al. [18], O'Keefe et al. [19], Lewis [20], and Hewitt [21], introduced the concept of disaster risk as the combination of the probability of occurrence of a hazardous event with its negative consequences. Over time, the concept of vulnerability has been interpreted in different ways [22–25] reflecting how the approaches and conceptual models for its assessment differ among the scientific research communities [23,25].

This paper proposes and applies a methodological approach that integrates physical exposure and social vulnerability into an "overall vulnerability index", to identify areas that can be negatively affected by climate- and marine-related processes, such as coastal erosion, landslides and sea level rise.

#### **2. Conceptual Framework**

Given the spectrum of possible conceptual interpretations of "vulnerability" and other related terms such as "exposure", it is deemed essential that, as of primacy, and prior to proposing the methodology adopted in this research for the overall assessment, the range of most-commonly used concepts of vulnerability and exposure coined in the contexts of both climate change and disaster risk management are first introduced here, before clarifying which concept and definition are adopted by this study. The disaster risk community identifies the notion of vulnerability as one of the elements required for risk assessment [26]. In this scientific context, "risk" is defined by exposure to a hazard, which is a potentially damaging physical event or phenomenon, and vulnerability, which denotes the relationship between the severity of the hazard and the degree of damage caused to the exposed element [27]. Specifically, according to the United Nations International Strategy for Disaster Reduction [26], vulnerability is defined in a qualitative way as "the characteristics and circumstances of a community, system or asset that makes it susceptible to the damaging effects of a hazard". In the same document, exposure is also defined qualitatively as "people, property, systems, or other elements present in hazard zones that are thereby subject to potential losses".

In earlier definitions, vulnerability was considered as the degree of loss of an element at risk after the occurrence of a natural process [28]. On one hand, the above-mentioned and more recent (but not the latest) definitions consider the notion of vulnerability as denoting a pre-existing condition related to the characteristics of the elements at risk and gives less emphasis to the process [29]. On the other hand, in the Climate Change Adaptation (CCA) context, more focus is placed on the concept of vulnerability defined as a function of exposure, sensitivity and adaptive capacity [10,30]. Therefore, in the risk context, the definition of vulnerability is quite similar to that describing the sensitivity of the system's components in the climate approach [10], while in the climate change community it is defined in a similar way to the concept of risk used in the Disaster Risk Reduction (DRR) context [31].

Only in more recent years have the United Nation Institutions (UNDRR, IPCC) had a key role in converging on a common and shared definition of the vulnerability concept in CCA and DRR fields. In fact, the IPCC integrates the different conceptualizations of vulnerability and provides its upgraded and clearer definition in the glossary of the Fifth Assessment Report (AR-5, [11]), in which vulnerability expresses "the propensity or predisposition to be adversely affected. Such predisposition constitutes an internal characteristic of the affected elements (or societies) and includes the concepts to cope with, resist, and recover and the lack of capacity to cope and adapt to the adverse effects of a physical event". Indeed, in this latter AR, the definition of risk (as a whole) is expressed as being the result of the interaction between vulnerability, exposure and hazard, whereby the distinction between the contribution of vulnerability and exposure to risk has been made clearer, which is a view also shared with the DRR community. This paper adopts a method of combining "physical exposure" information and "social vulnerability" data, which is in line with this evolution and convergence of concepts of vulnerability and exposure, across the above discourses.

In parallel to the vulnerability concept development, several actions taken in the last couple of decades have also focused on the development of methods and tools for supporting decision-makers in the reduction of coastal hazards' impacts. Particularly noteworthy is the index-based approach, which was introduced by Gornitz at the beginning of the 1990s [32,33]. This is mainly based on indirect analysis supported by photointerpretation and topographic maps and is still the most commonly used method for coastal vulnerability assessment at regional levels [34–40].

At the European level, a number of transnational research projects have been funded with the aim of improving scientific knowledge and increasing the awareness of decision-makers on coastal vulnerability, the potential damages of marine processes and the effectiveness of coastal adaptation strategies (such as EUROSION, MICORE, RISC-KIT, ANYWHERE, OPERANDUM). In particular, the main outputs of the RISC-KIT Project are represented by a set of tools that allow identification of the most vulnerable areas (hot-spots) by means of a coastal index [41,42] and then selection of suitable measures for increasing coastal adaptation and therefore favouring risk reduction [43].

The evolution and emergence of concepts of vulnerability and exposure as being different features is warranted in certain situations, to provide a better understanding of "why" and "how" societal assets and values may be under threat or at risk under different scenarios of exposure. This can provide guidance to planners and policy makers in identifying suitable adaptation measures and actions, but it also makes risk assessment far more complex.

Some factors that contribute to the vulnerability of society can be expressed in terms of the dependence of societal wellbeing on certain exposed physical assets or other values that are at risk. Examples of this are infrastructure and utilities, which underpin societal wellbeing to different degrees, depending on their nature. This component of vulnerability is distinct from vulnerability associated

with the inherent characteristics of society. A community relying on physical elements for its wellbeing can be deemed as being vulnerable to the extent (or level) of its dependence on those specific physical elements located in a hazard-prone area. This creates different levels of "physical exposure", in direct relation to their level of importance and to the degree to which societal wellbeing is dependent on and would be affected by losses to them. Thus, "physical vulnerability" is considered as one part of "overall vulnerability", while "social vulnerability" (attributed to inherent societal factors) is considered as another component of "overall vulnerability".

In this context, the research presented here aims to evaluate the overall vulnerability, and to refine the analysis and understanding of two distinct components, physical vulnerability (representing exposure from dependence on physical elements) and social vulnerability to a given set of external climate- and marine-related processes. Mirroring physical exposure as physical vulnerability and combining it with social vulnerability into an overall vulnerability index, an expedient and cost-effective method for the identification (and assessment) of the overall vulnerability of coastal areas at risk is provided.

The selected coastal study area, located along the NE sector of the Island of Gozo (Maltese archipelago), is one for which considerable applicative research has already been undertaken in order to showcase the high geological and geomorphological significance of its coastal landscapes [44]. Although the area is known to be particularly susceptible to several coastal hazards (erosion, storm surges, floods, sea level rise, landslides), it still lacks a detailed and refined appraisal of its natural and anthropic coastal elements, based on an analysis of assets and values that could potentially be at risk. Approaching the assessment of vulnerability of this coastal area from the perspective of the two components of social vulnerability and physical vulnerability provides a sharper analytical insight into overall vulnerability. Furthermore, the analytical methods and tools used in this study provide a methodological template for overall vulnerability assessment in the form of an Overall Vulnerability Index (OVI) that can be computed and used, with relative ease at different scales.

The scientific literature review carried out as part of this study reveals a gap in the area of assessment of vulnerability. Studies carried out to date have focused on the overall evaluation of hazard and susceptibility along the Maltese coastal sectors [45], rather than on the methodological assessment of vulnerability per se. Integrated approaches have been applied along the NW coast of Malta for landslide hazard assessment [46–49]. Landslide susceptibility assessment assisted by Persistent Scatterers Interferometry (PSI) was carried along the same stretch of the NW coast of Malta [49,50]. Based on the assumption that the identification and analytical mapping of the exposed elements represents a key step for the evaluation of the coastal risk, this study tries to bridge this gap by proposing a relatively simple yet reliable, cost-effective and easily replicable procedure for the assessment of coastal vulnerability. Nonetheless, while not being overly complex, the proposed approach provides an analytical discernment of different components of overall vulnerability that, when mapped geographically, provide sufficient guidance as to which areas are most vulnerable and on what account.

Following on from the methodological approaches proposed in previous studies [42,51–54], this study formulates and applies a research method for overall vulnerability assessment, which is based on the following steps: (i) identification of the main exposed elements (natural and anthropic) located in the investigated area; (ii) definition of their relative exposure level, in economic and ecological terms; (iii) assessment of the social vulnerability of the population living in the investigated area; (iv) calculation of the overall vulnerability by means of a combined index. As highlighted in [55] the use of an index as an evaluation tool requires the definition, weighting and aggregation of a number of indicators, which are defined as variables that are "an operational representation of a characteristic or quality of a system" [56,57].

In detail, the method applied here foresees a set of indicators for the evaluation of the exposure level of local land uses (considering both the presence of economic activities as well as natural protected areas) and anthropic assets (such as transport networks and main utilities). Furthermore, the social context that characterizes the investigated area is also taken into account, in order to include social vulnerability in terms of population capacity to respond to and cope with a hazardous event in the overall evaluation of vulnerability. The combination of physical and social vulnerability provides the overall index, which expresses the level of overall vulnerability.

The use of distinct physical and social vulnerability indicators, and the definition of an Overall Vulnerability Index (OVI), has the main advantage of summarizing complex issues making them more easily interpretable, facilitating decision making and communication among stakeholders [58]. Meanwhile, the distinction between the two components at analysis stage provides a more in-depth understanding of specific anthropic elements that contribute to different risk levels.

#### **3. Study Area**

#### *3.1. Geological and Geomorphological Setting*

The Maltese archipelago is located in the central Mediterranean Sea and comprises the main islands of Malta, Gozo and Comino (Figure 1). The coastal landscape (Figure 2) is the result of long-term evolution under the influence of tectonic activity, geomorphological processes and sea level oscillations. Due to their geological and geomorphological setting, these islands are particularly prone to different marine-related and gravity-induced processes such as landslides, coastal erosion, storm water runoff and sea level rise, enhanced by ongoing climate change. Multidisciplinary research and integrated investigations have been carried out to better understand the evolution of the geomorphological processes within the archipelago in the wider dynamic scenario of ongoing climate change as a way towards the reduction of risks associated with these processes [47,59–61].

**Figure 1.** Geological setting of the Maltese archipelago (cf. [62]) and location of the study area, bounded by the NE coastline of Gozo and the red dashed line parallel to it on the inland side (after [44], modified).

**Figure 2.** Main geomorphic features of the study area: (**a**) block slides (west of Da*h*¯let Qorrot Bay); (**b**) rock fall at the bottom of a limestone plateau and earth flow/slide affecting the underlying clayey terrain (between San Blas Bay and Da*h*¯let Qorrot Bay); (**c**) plunging cliff between Da*h*¯let Qorrot Bay and Ras il-Qala; (**d**) sloping coast between Da*h*¯let Qorrot Bay and Ras il-Qala; (**e**) shore platform east of Marsalforn Bay; (**f**) Ramla Bay pocket beach; (**g**) cliff shaped in Blue Clay east of Marsalforn Bay; (**h**) built-up coast of Marsalforn Bay. After [44], modified.

Previous studies of the Island of Gozo have focused on general aspects of coastal features [60,63–66]. Only a few papers deal with the specific geomorphological aspects of the island [67–69], some of them referring to its rich geoheritage [70–73]. A detailed geomorphological map of the investigated stretch of coast (NE Gozo), based on geomorphological and geological field surveys integrated with the analysis of marine geophysical data, has been published recently [44]. The geomorphological map is accompanied by two other maps that show land use and the distribution of coastal geomorphotypes (Figure 3). Such documents constituted the basis for carrying out the vulnerability analysis presented below.

**Figure 3.** Distribution of coastal morphotypes in the study area (after [44], modified).

From a geomorphological perspective, the investigated coastal stretch is characterized by limestone plateaus bounded by steep structural scarps which are progressively reshaped by gravitational and/or degradation processes. Clayey slopes, located at the foot of the limestone plateau, accommodate terraced fields of actively used or abandoned agricultural land. Numerous blocks of rock are strewn over the clayey terrain (a unique landscape known as *rdum* in Maltese) that slopes more gently away from the plateau edge.

The investigated coastline is characterized by the alternation of inlets and promontories. The accumulation of sand and mixed grainsize deposits results in the formation of pocket beaches where this corresponds with bays and coves. The large sandy beach of Ramla il-Ħamra Bay ('red sandy beach') is a particularly noteworthy example. It is also partly bounded by dunes on the landward side and is protected as a Special Area of Conservation (SAC) within the Natura 2000 network. The main inhabited centre of the investigated stretch of coast is Marsalforn, with its homonym bay, the latter being intensely urbanised and anthropized, comprising a predominantly built-up coastline.

The landscape of the investigated stretch of coast, and in general of the entire Maltese archipelago, is largely controlled by the different erodibility of the exposed lithostratigraphic units constituted by marine limestones and marls of the late Oligocene-Miocene [74,75]. The outcropping geological formations comprise (from the oldest to the youngest): Lower Coralline Limestone Formation (late Oligocene, Chattian), Globigerina Limestone Formation (late Oligocene—middle Miocene, late Chattian—Langhian), Blue Clay Formation (middle-late Miocene, Serravallian—Tortonian) and Upper Coralline Limestone Formation (late Miocene, Tortonian—early Messinian).

The Lower Coralline Limestone Formation consists of pale grey, hard, shallow marine biomicrites and biosparites [76], and outcrops in a restricted coastal stretch forming subvertical cliffs. Few Lower Coralline plunging cliffs (*sensu* [77]) are found in the investigated stretch of coast. They host roof notches with an asymmetric shape, recognized by Furlani et al. [78]. Lower Coralline Limestone is more commonly found in sloping coast formations, as typified in the stretch between Dahlet Qorrot Bay and Ras il-Qala.

The Lower Coralline Limestone Formation underlies the Globigerina Limestone Formation, the latter being younger and more erodible with respect to the former. The Globigerina Limestone

Formation consists of flattened areas along the coast, and includes yellowish, fine-grained, planktonic foraminiferal limestones. This formation features prominently as shore platforms along the investigated stretch of coast, where it outcrops above sea level and where the overlying softer blue clay layer has been eroded away, with examples at the Rdum il-Kbir promontory and on the eastern side of Marsalforn Bay. The Blue Clay Formation, overlying the Globigerina Limestone Formation, consists of grey, soft marls, clays and silty sands, forming gentle slopes. Sea cliffs shaped in Blue Clay can be also found within the study area. Finally, the youngest unit outcropping in the study area is the Upper Coralline Limestone Formation. This layer forms the plateaus at the top and is frequently weathered into steep cliffs and well-developed karst topography. It is also the source of the blocks of rock strewn onto the clayey slopes that slope away from the plateau edge more gently. A schematic stratigraphic coloumn is shown in Figure 1.

The NE coast of Gozo is particularly prone to gravity-induced processes such as rock spreads, block slides, rock falls and earth flows/slides. This is mainly due to the tectonic and geological features of the area which is characterized by a dense network of joints and fractures [79] and by the superposition of terrains with different geomechanical behaviour (cf. [80]). In fact, the brittle limestone plateaus overlying clayey terrains enhance the fracturing of the plateaus and the development of lateral spreading locally evolving into block sliding [44,46,81–85].

Moreover, clayey slopes are more prone to shallow earth flows and earth slides [86], while limestone cliffs are affected by rock falls, which has caused scarp retrogression over time.

Geomorphological investigations, relevant to the submerged area along the investigated stretch of coast [44], have revealed that block slides and earth flow/slide runout continue locally below the sea level, reaching ca. −20 m of depth. This is also in agreement with evidence from other submerged areas along the Maltese Islands [44,87]. Evidence from landslide dating in the NW coast of Malta infer that these deposits were emplaced in a subaerial environment during sea level lowstands and subsequently became submerged during the post-glacial marine transgression [61].

#### *3.2. Social, Economic and Tourist Setting*

Studies that include geoheritage assessment and geosite inventory highlight that the Maltese archipelago is considered as an attractive geotourist destination due to the strong interaction between natural and cultural aspects (cf. [71,88]). Data from the World Travel and Tourism Council (WTTC) for the year 2017, retrieved from Selmi et al. [88], show that 27.1% of Malta's GDP and 28.3% of total Maltese employment (corresponding to 55,000 jobs), were accounted for by activities directly related to and induced by travel and tourism [89,90]. Taking the type of tourism into consideration, data referring to 2017 show that the majority of tourists (almost 85%) visited Malta for holidays, while a very low percentage visited the country for business and other purposes (8% and 7%, respectively). More recent data show that the number of tourists is constantly increasing (standing at 2.6 million in 2018). The WTTC forecast that the travel and tourism contribution to national GDP will rise to 34.6% by 2027.

The tourism on the Island of Gozo is highly dependent on tourism activity of the main island, Malta [91], given that all of Gozo's tourist traffic necessarily passes through mainland Malta. Tourism is a significant source of income and employment and it is one of the primary contributors to the Gozo economy. Gozo's tourism also relies heavily on domestic tourists (with 400,000 domestic tourists or visitors per year coming from Malta) and on one-day trips of international tourists [92]. In 2018, almost 100,000 guests including resident and non-resident spent an average of three to four nights in one of the accommodation facilities of Gozo [93].

The Gozo Island attracts many tourists, especially during the summer period, for its environmental, cultural and geological heritage. Moreover, senior and middle-aged foreign residents more commonly chose the Island of Gozo as a place for their retirement, mainly for its relative peacefulness and quiet, combined with its array of scenic features (Gozo is colloquially known as "the place where time stood still").

The coastal sector investigated by this study includes four administrative districts: Zebbu ˙ g˙ (7.6 km2), Xag*h*¯ra (7.6 km2), Nadur (7.2 km2) and Qala (5.9 km2), which are characterized by extensive urban development. The study area is host to two renowned tourist destinations, Marsalforn Bay and Ramla Bay, both of high importance to Gozo tourism. Marsalforn Bay is one of the main inhabited centres located on the Gozo coast. The availability of a significant number of accommodation facilities, shops, restaurants and several diving centres contribute to a dense tourist population and to lucrative businesses in Marsalforn, especially in summer.

Total guests in the Gozo and Comino region increased by 12.6% to 97,781 in 2017, while total nights spent went up by 11.1% to 347,943 when compared to the previous year [94]. Within a continuous upward trend in the last five years or so, the Gozo and Comino region recorded a strong growth in terms of inbound tourist arrivals in 2017 [94]. It is worth noting that among the top five localities where inbound tourists to Gozo stayed longest [94], two are within the study area, with important tourist destinations such as Marsalform and Ramla Bay.

Ramla Bay is the largest sandy beach in Gozo. It is characterized by golden-reddish sand, which lends it its name in Maltese (Ramla il-Ħamra) and makes this beach peculiar and unique to the Maltese Islands. The area around the beach is also of archaeological and historical interest, hosting Roman remains lying beneath the sand, a submerged seawall, fugass and battery defensive structures constructed by the Knights of St. John of Jerusalem in the mid-18th century, and the famous Calypso Cave looking over the western side of the beach.

Furthermore, the study area as a whole includes areas of high natural and ecological importance, hosting two wide natural protected zones, G*h*¯ajn Barrani (located west of Marsalforn, including Ramla Bay) and Il-Qortin tal-Magun u l-Qortin il-Kbir (located close to Dahlet Qorrot Bay), both included in the Natura 2000 network as SAC. Finally, it should be noted that the high economic value of the study area is also derived from specific land uses in certain areas, such as quarrying, which occupy a surface of 0.13 km2.

#### **4. Materials and Methods**

The method proposed and applied in this research for the assessment of coastal vulnerability is in line with the most recent index-based approaches generally used for this kind of analysis. In detail, the method relies on the outcomes of other research in the field [41,95] and is based on a new approach for the evaluation of indicators referring to: (i) the potentially exposed categories of assets (defined as "physical indicators"); and (ii) a number of parameters related to the social context (defined as "social indicators").

The indicators are interpreted here as operational representations (cf. [56,57]) of the physical and social characteristics of the area. Each physical indicator comprises different categories of assets, to which a score representing an increasing level of exposure was attributed. In order to tailor the approach to the local settings of the study area, the exposure level assigned to each category of assets was defined based on expert judgment. It is a widely shared opinion that the scoring and aggregation of indicators into indices may have a large impact on the resulting rankings and, consequently, on decision-making [55].

The spatial overlay of the physical indicators provides an estimate of the physical vulnerability level, which ranges from very low to very high vulnerability. Meanwhile, the various social indicators refer to vulnerability of the socio-economic aspects. These are also classified into five levels ranging from very low to very high, to provide a measure of the social vulnerability level for the investigated area. The overlay of the two sets of aggregated physical and social indicators enables an overall zonation that shows grades of the combined physical-socioeconomic vulnerability, defined as "overall vulnerability", thus identifying which are the most vulnerable stretches of the investigated coastal area.

Specifically, the method applied here comprises the following four main steps:


The analyses are supported by GIS tools, which allow identification of the exposed coastal assets (i.e., natural and semi-natural environments, buildings, infrastructure, and agriculture), their combination with social data, and the calculation and mapping out of the overall vulnerability levels.

#### *4.1. Definition of the Landward Limit of the Coastal Area to be Investigated*

The area investigated was delimited according to the definition of a RICE area (Radius of Influence of Coastal Erosion and Flooding) proposed in the framework of the EUROSION project [95] for the identification of the coastal areas potentially impacted by marine-related process. The limit was set at a maximum distance of 500 m inland, or reaching a maximum of 5 m a.s.l. from the coastline. Furthermore, a buffer area of 100 m inland from the edge of the scarps formed in the Upper Coralline Limestone—which are extensively affected by landslide processes—was also added to the area defined in the above manner, in order to specifically account for landslide hazard. Information about landslide distribution provided by the geomorphological map of the study area (cf. [44]) was used specifically for this purpose, so as to identify the areas prone to slope instability.

The resulting landward limit was simplified and shown as a dashed red line in Figure 3, in order to provide a more linear indication of the investigated coastal sector.

#### *4.2. Classification of Physical and Social Indicators*

This step comprises the identification of data required for the evaluation and classification into five different levels, from very low to very high, for each of the two sets of indicators (physical and social indicators).

The evaluation of the anthropic and natural assets potentially exposed to coastal hazards is based on the elaboration of the following indicators: land use, transport network, and utilities. A GIS layer has to be created for each indicator in order to identify polygons representing the spatial unit to which the related exposure level is assigned, indicating its specific numerical value, ranging from 1 (very low exposure) to 5 (very high exposure). These layers are first converted to a raster format (5 × 5 m) and then overlaid to estimate the physical vulnerability level, by considering the maximum exposure value assigned to each cell.

At the detailed level, the land use information was collected from the land use map available in Prampolini et al. [44], in which land use categories were classified as: (i) natural and semi-natural areas, (ii) agricultural areas, and (iii) artificial surfaces. This land use map was supplemented with additional detailed information in order to include strategic elements, such as civil protection posts, police stations, emergency posts, fire corps posts, port authorities' posts, hospitals, and schools. The identification of these elements was supported by photointerpretation of the most accurate maps available for the study area and verified in the field during the investigation.

An exposure level (ranging from 1, very low exposure, to 5, very high exposure) was assigned to each land use category based on expert judgement. The highest exposure level was attributed to the areas occupied by settlements, constructions, and human activities while the lowest level was assigned to abandoned agricultural zones. Details concerning the land use indicator and the exposure levels related to each category are shown in Table 1.


**Table 1.** Physical indicators and expert-based exposure level (ranging from value 1 to value 5) assigned to each category of assets.

The transport network information was retrieved from Open Street Map downloaded from [96] as a shapefile. In this case, the roads were classified as follows: primary road, secondary road, tertiary road, inhabited street, residential, services, footway, path, track, and steps; the greater exposure level was attributed to roads of national or international importance. For each polyline element, a buffer distance was built (up to 20 m of total width for the primary road) in order to take into account a proper area of pertinence. Details concerning the exposure levels defined for the transport indicator are shown in Table 1.

The electricity network was taken into account for the utilities. The spatial distribution of these elements is available on Malta Inspire Geoportal [97]. In this case, a buffer zone up to 10 m width (10 m for substations, feeder pillars, and overhead lines, and 5 m for street lighting) was considered for each element. Details concerning the exposure levels associated with each electricity element are indicated in Table 1.

The social indicators allow characterization of the districts located in the investigated coastal area by evaluating, directly and indirectly, the social characteristics of the population living in the zone and therefore prone to be affected by coastal hazards. Generally, the social vulnerability information is obtained from census data. In this study, social data was downloaded from the Inspire Geoportal [97] and integrated with information available in the special report recently published by the Malta National Statistics Office (NSO) [98]. The social indicators are gauged from the economic budget allocated by the Government for supporting the population. The social vulnerability analysis carried out here is based on the assumption that higher allocation of social schemes corresponds to higher vulnerability of the population living in the investigated districts.

Specifically, the following indicators were taken into account: health care, disability, old age, children, and unemployment. These aspects represent some of the higher-risk groups of society, which social services and budget allocations aim to protect. The NSO provides a classification of these indicators in five classes that were here converted in a numerical value ranging from 1 to 5. The available data is provided by the NSO at district level and thus a numerical value was assigned to each district as a whole. The social vulnerability was then calculated as geometric mean of the values attributed to each indicator and then classified into five levels (from 1, for very low vulnerability, to 5, for very high vulnerability) by means of an equal interval classification method. A detailed description of the social indicators used in this study is reported in Appendix A.

The physical and social indicators identified are weighted equally, as in most of the indicator-based studies [99], meaning that each individual indicator has the same influence on the final calculation of the overall vulnerability.

#### *4.3. Data Overlay and Computation of an Overall Vulnerability Index*

Once all the required data related to the physical and social vulnerability are collected and expressed in five levels, they are overlaid by means of a specific GIS tool. The overall vulnerability calculation is therefore the result of the aggregation of the physical vulnerability and social vulnerability data, the combination of which provides the definition of the Overall Vulnerability Index (OVI) as follows:

$$\text{Overall Valnerability Index} = \text{(Physical vulnerability} \times \text{Social velnerability)} \tag{1}$$

Finally, the Overall Vulnerability Index is classified into five levels of increasing vulnerability (from very low to very high vulnerability) by means of the equal interval classification method.

#### **5. Results**

The results of this study are represented by: (i) GIS-based maps that show the spatial distribution of the natural and anthropic exposed elements and their classification in terms of physical vulnerability; (ii) social vulnerability level estimated for each investigated coastal district; (iii) overall vulnerability map of the NE part of the Gozo Island; (iv) areal extent and relative percentage of each vulnerability level. The land use classification (Figure 4) enables the evaluation of the surface occupied by different categories of natural and anthropic elements (cf. Table 1 in Section 4.2).

**Figure 4.** Detailed land use map. The red dashed line identifies the inland boundary of the study area. The districts' boundaries are also indicated.

In detail, 3.6 km<sup>2</sup> are occupied by agricultural areas, 1.2 km2 by abandoned agricultural areas, 1 km<sup>2</sup> by terraced agricultural field, 1.2 km<sup>2</sup> by natural protected areas, 0.8 km2 by bare rock, 0.7 km2 by natural and vegetated areas, and, 0.5 km<sup>2</sup> by residential areas. Finally, quarrying and land occupied mainly by agriculture with significant areas of natural vegetation occupy 0.1 km2. The remaining surface area (corresponding to 602,268.2 m2) is occupied by the following categories: strategic elements (including a Police station), historical/archaeological sites (The Tower of Ta Sopu, west of Da*h*¯let Qorrot Bay, and Saint Anthony's Battery at Ras il-Qala), green urban areas, beaches, dunes, sand, abandoned terraced fields, entertainment zones, and a cemetery. These values are indicated as percentages in Figure 5.

**Figure 5.** Percentage of the area occupied by each land use category according to the classification shown in Figure 4.

The interpretation and spatial representation of the overlaid physical indicators (cf. Section 4.2) (cf. Table 1) is representative of the different levels of physical vulnerability, given that the physical exposure levels are correlated directly to the level of value of the assets at risk.

This enables the zonation and calculation of what are, in effect, the extent of the areas with different levels of physical vulnerability, expressed in square kilometres and in percentages of the total surface investigated. The spatial interpretation of physical vulnerability is mapped in Figure 6 while the corresponding numerical values are reported in Table 2.

The numerical values assigned to each of the social vulnerability indicators for the investigated coastal districts are summarised in Table 3.

**Figure 6.** Physical vulnerability map resulting from data overlay of the physical indicators (land use, transport and utilities). The red dashed line identifies the inland boundary of the study area.



**Table 3.** Social indicators evaluated for the four coastal districts of the study area. The method for attributing the numerical value assigned to each indicator is indicated in the methodological paragraph (cf. Section 4.3).


The aggregation of the social vulnerability indicators (cf. Table 3) enables the evaluation of the social vulnerability level for each of the investigated districts (Figure 7). Specifically, two districts (Nadur and Xag*h*¯ra) are characterized by medium vulnerability and two (Qala and Zebbu ˙ g) by ˙ high vulnerability.

**Figure 7.** Social vulnerability map of the four administrative districts located in NE Gozo. The red dashed line indicates the inland boundary of the study area.

Finally, the overlay of the physical vulnerability levels with the social vulnerability levels enables a computation (and mapping) of the Overall Vulnerability Index for different polygons in the study area (cf. Section 4.3). It provides the basis for evaluating the variation in the overall vulnerability levels across the investigated coastal sector and enables this to be represented spatially (Figure 8).

**Figure 8.** Overall vulnerability map resulting from the spatial aggregation of the physical vulnerability levels and the social vulnerability levels over the area. The red dashed line indicates the inland boundary of the study area.

In sum, 1.9 km2 are occupied by areas showing a low overall vulnerability level, 5.6 km2 are occupied by areas showing a medium overall vulnerability level, 0.7 km2 are occupied by areas with a high overall vulnerability level and, finally, areas characterized by a very high overall vulnerability level cover 1 km<sup>2</sup> (Table 4).

**Table 4.** Areal distribution of overall vulnerability levels that are the result of the combination and overlay of spatial distributions of the physical and social vulnerability factors.


#### **6. Discussion**

The proposed index-based method allows zoning of the investigated coastal stretch into different levels of vulnerability to climate- and marine-related processes. Results show that the study area is divided in four zones only (from low to very high vulnerability), since there are no areas with a very low level of vulnerability.

The coastal sectors located east of Marsalforn Bay, which include Ramla Bay and the area on the western side of Da*h*¯let Qorrot Bay, show the highest overall vulnerability level. This is due to the combination of very high physical vulnerability levels pertaining to the presence of two Natura 2000 sites (cf. Sections 3.1 and 5) and high social vulnerability levels for the Xag*h*¯ra and Nadur districts. The high social vulnerability of Xag*h*¯ra and Nadur is explained by the fact that each of the latter districts accounts for more than 4000 inhabitants, which is twice the number of inhabitants of the two other districts considered, Zebbu ˙ g and Qala, for which a medium social vulnerability is assigned. ˙

The coastal sector surrounding Marsalforn Bay is mainly characterized by a medium level of overall vulnerability as a result of the combination of a medium physical vulnerability level, explainable by the presence of residential areas, and, a medium social vulnerability level obtained for the Zebbu ˙ g district. ˙ Furthermore, the innermost sector of the study area shows a prevailing medium overall vulnerability level, resulting from the combination of a low to very low physical exposure level, owing to the presence of cultivated and abandoned agricultural fields respectively, and a high social vulnerability level resulting for the districts of Xag*h*¯ra and Nadur, as discussed above. The eastern part of the investigated sector, with bare outcrops of rock, hosts areas characterized by a low overall vulnerability.

The index-based method here proposed can be considered a suitable approach for the identification of coastal areas that are most vulnerable to consequences of different climate- and marine-related processes. Since the method adopted here combines vulnerability and exposure, the results represent an overall vulnerability assessment. It relies on the evaluation of the two main components that are used to define the overall vulnerability: physical and the social vulnerability, thus accounting for both the physical assets potentially exposed to climate and marine processes and the social aspects of the local population. This approach is in line with the convergence of different terms by the IPCC in its Fifth Assessment Report (AR-5, [11]).

The application of the proposed OVI method has allowed identification of a number of operational advantages. First of all, as the analysis is based on a sequence of steps, the method is relatively simple and easy to apply. It has been found to be cost-effective on account of the possibility of managing and processing the acquired data in a GIS environment with relative ease. Furthermore, it does not require intense field work, since it relies on indirect analysis, supported by existing geological-geomorphological information that is either readily available in the literature or easily collected. It can therefore be applied to wider coastal stretches, as in the case of other index-based approaches [40,100–103]. It should be underlined that the applied method can take advantage of data

sets and information concerning the natural and anthropic assets which are generally freely accessible and downloadable, e.g., from national geoportals available in most of the countries (at least across Europe), and from open databases, such as OpenStreetMap. Therefore, the method may have a wide range of applicability at different scales. Large scale analyses are based on detailed information about exposed/vulnerable elements and population (such as the land use and vulnerability maps proposed in this study). Meanwhile, small scale vulnerability maps concern wide areas, taking into account only the most prominent and spatially persistent exposed/vulnerable elements and regional data about the exposed/vulnerable population (as in the case of EUROSION Project), which accounted for all the European countries and provided a comprehensive European-level data repository at scale 1:100,000 [95].

The OVI method provides a sound approach suitable for the identification (and assessment) of vulnerable areas and sectors, even as an expedient and cost-effective scoping stage to identify which area may be analysed more specifically and at greater expense. It provides a useful tool for an overall time-efficient and cost-effective approach to focus limited research resources onto where they are needed most.

The presented research is a pilot-study and a first-ever combined overall vulnerability assessment carried out for any area in the Maltese archipelago, and it is promises to provide new contributions at different levels and in different ways. At the localised level, the proposed method, and the results obtained from it, are promising as they reveal the potential applicability of this method for other (and potentially wider) coastal areas of the Maltese archipelago. Moreover, even at the local level, the overall vulnerability assessment could be easily updated by using any newly available or more detailed data regarding the social vulnerability of the investigated districts, as and when this becomes available. The relatively simple combined approach makes it more possible to keep the OVI up to date, even at the localised level. By way of example, it is relevant that during 2016, Malta's social protection outlay rose by €70.3 million in comparison to 2015 and that Old Age and Sickness/Health Care witnessed the biggest increases in social outlay [98]. The NSO [94] published statistical data on the basis of six different regions, from which localities can be studied specifically as socio-economic parameters change. This means that the social vulnerability level of the investigated district could be expected to rise in the next decades, influencing the overall vulnerability assessment. At the national level, the method could be applied to provide a first-round assessment of overall coastal vulnerability, to identify the most vulnerable stretches of coastal areas and values, and then to follow this up with a more detailed qualitative and quantitative analysis and comparison of risk assessments for specific areas and sectors, for different hazard types.

Running parallel to other research at the European level, further studies could include the spatial susceptibility analysis, including each of the coastal hazardous processes already identified here as affecting the investigated stretch of coast (e.g., erosion, sea level rise, landslides). As already done for other coastal zones in Europe [39,40,54], the investigated area should be classified into zones with different susceptibility, accounting for their proneness to be potentially affected by the specific impacts from extreme events and related processes pertaining to the particular hazard types (e.g., erosion, sea level rise, landslides). The susceptibility to different climate- and marine-related hazards could then be coupled with the vulnerability data derived by this study to perform a complete risk analysis (as done for example in [25,53]).

In this context, it is worth noting that the Maltese archipelago is situated centrally in the Mediterranean Sea, which has been classified as one of the regions most sensitive to climate change and, therefore, it is considered as a hot-spot area [104]. Climate change projections for the Mediterranean region [105] show that the area is experiencing an increasing temperature with consequent change in spatial and temporal distribution of weather/climate extremes [106]. More intense events, with alternating and more severe drought and precipitation, are expected to trigger and exacerbate erosive and mass movement processes [11,13], affecting the spatial distribution of the susceptibility to these types of events and hazards.

A direct consequence of global temperature increase is the rise in sea level, the direct impacts of which on coastal systems become manifest with larger rates of erosion, increase in flooding events, wetland loss, and saltwater intrusion [107]. Studies regarding the reconstruction of sea level changes in the Mediterranean Sea during the last 2000 years have shown that, in tectonically stable Mediterranean areas, the sea rose about 1.1 m [108,109] while for the last two decade the estimated rise was of about 3 cm/decade [110]. However, differential vertical land movements, including uplift and subsidence, characterize the Mediterranean coasts and, for this reason, the trends of sea level rise in the Mediterranean Sea have a large spatial variability [111]. Taking into account the role of terrestrial ice melt, steric effects and glacial isostatic adjustment, the future total Mediterranean averaged sea level rise has been estimated to be between 9.8 and 25.6 cm by 2040–2050, depending on the Representative Concentration Pathways scenario [112]. Sea level projections, obtained by coupling modelled eustatic trends with local ground movements, are in general used for supporting the assessment of the future coastal risk to sea level rise, which is aimed at reducing its impacts on natural coastal environments and human economic activities [3].

Based on the assumption that the investigated coastal area could potentially be affected by all the above mentioned climate-related hazards, further research activities should focus on risk analysis and mapping, as already done with reference to hazards related to sea level rise in several parts of Europe [113,114] and worldwide ([115–119] and reference therein). In this context, the results obtained in this study can prepare the ground for a comprehensive risk assessment, combining the identification of the exposed assets and the evaluation of vulnerability levels with a quantitative assessment of the hazardous processes affecting the area.

Finally, it is worth noting that the study is in line with the methodological approach proposed by the European Environment Agency (EEA) for the identification and implementation of climate change adaptation strategies (European Climate Adaptation Platform Climate-ADAPT). Specifically, in the EEA approach, six steps of analysis are indicated: (i) preparing the ground for adaptation; (ii) assessing risk and vulnerability; (iii) identifying adaptation options; (iv) assessing adaptation options; (v) implementation; (vi) monitoring and evaluation. The vulnerability analysis shown here represents a useful tool for addressing Step 1 and Step 2, while further activities should be planned with the aim of developing and tailoring the most suitable adaptation actions for the protection of the natural ecosystem and the maintenance of the anthropic activities.

#### **7. Conclusions**

This research represents a first attempt at the evaluation of coastal overall vulnerability in the Island of Gozo (Maltese archipelago) to the main climate- and marine-related processes, such as coastal erosion, landslides and sea level rise, to which the island is particularly prone. The analysis method developed and tested in this study is based on a conceptual interpretation of overall vulnerability, defined as the combination of two distinct components, namely, physical vulnerability and social vulnerability.

The study departs from a discussion of different notions of risk, including commonly used notions of exposure and vulnerability, which are defined qualitatively in the most authoritative sources (cf. Section 2). Given the identified, often conflicting understanding of the concept of vulnerability by various sectors of researchers and policy makers, it was determined that an in-depth consideration of the vulnerability concept was required prior to determining the vulnerability assessment methodology.

The study proposes a method for the assimilation, analysis and computation of overall vulnerability, and for its graphical spatial representation. The method manages to converge different qualitative notions of physical exposure and social vulnerability, and also makes it possible to derive a spatial quantitative distribution of the Overall Vulnerability Index. The latter can be represented in a map showing different coastal vulnerability levels that are easily readable spatially, making it simple to communicate comprehensibly to decision makers.

The Overall Vulnerability Index is calculated by means of an index-based method that is proposed along the lines of approaches developed in the framework of previous research projects and adapted specifically for the context of the study area and to the available information.

The study area was identified as possessing features of high economic value derived from tourist and mining activities and natural protected areas, altogether making coastal vulnerability a major concern. The final results of the analysis reveal that most of the investigated area (61.3%) is characterized by a medium level of overall vulnerability. A very high overall vulnerability level (11%) was obtained for the areas located east of Marsalforn Bay and close to Da*h*¯let Qorrot Bay, including Ramla Bay, mainly owing to the presence of two sites protected as Special Areas of Conservation within the Natura 2000 network (G*h*¯ajn Barrani and Il-Qortin tal-Magun u l-Qortin il-Kbir sites). A high vulnerability (7.3%) resulted for the main roads, while 20.4% of the total surface is characterized by low vulnerability areas mainly corresponding to abandoned agricultural fields and bare rocks outcrops.

The type of analysis presented here can be easily replicated for other coastal regions because the data sets required for the proposed method are almost always freely available or relatively easy to compile. Furthermore, the wider application of this method to the entire coastal area of the Gozo Island as well as to the coastal zones of the Maltese Islands has the potential to contribute to an overall risk characterization for the entire Maltese archipelago and it provides a useful tool for the identification of the most exposed and vulnerable zones (hotspot areas) that require action for their protection as a matter of priority. Starting from readily available data, which can be processed and mapped in a GIS environment with relative ease, the method proposed and tested here can provide policy makers, as well as coastal management agencies, with a graphical representation of easily comprehensible and useful data to support decision-making processes at both operational (short-term) and strategic (medium-long-term) level, enabling them to devise adaptation and structural protection measures in a more specific and effective manner.

Moreover, the method is cost-effective and time-efficient on account of ease of data processing, and it can also be a precursor for more detailed qualitative and quantitative analysis of risk (from different hazard types) at different scales for the areas or sectors considered to be at greater risk.

Finally, beyond the site-specific results, the method represents an important contribution toward more comprehensive risk assessment, also in terms of its potential transferability and replicability to different hazard types, including the local effects of climate change from extreme weather/climate events and sea level rise that need to be taken into consideration in a timely and effective manner.

**Author Contributions:** Conceptualization: A.R., V.V., G.B., A.S.M., M.S.; methodology: A.R., V.V., G.B., A.S.M., M.S.; formal analysis: V.V. and A.R.; investigation: M.S., V.V., A.S.M.; data curation: A.R. and V.V.; writing original draft preparation: A.R. and V.V.; writing review and editing: A.R., V.V., G.B., A.S.M., M.S.; visualization: V.V. and A.R.; supervision: M.S.; project administration: M.S. and A.S.M.; funding acquisition: M.S. and A.M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study has been carried out in the frame of the Project "Coastal risk assessment and mapping" funded by the EUR-OPA Major Hazards Agreement of the Council of Europe (2020–2021). Grant Number: GA/2020/06 n◦ 654503. Scientific responsible: Anton Micallef (ICoD) and Mauro Soldati (Unimore).

**Acknowledgments:** The authors are grateful to Chris Gauci (Research and Planning Section, Marine and Storm Water Unit, Public Works Department, Floriana, Malta) for information provided about the social vulnerability aspects. The authors would also like to thank Chiara Bordin for helping with data processing in GIS environment.

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

#### **Appendix A**

The evaluation of the social vulnerability for the coastal districts located in the investigated area was based on the calculation of a synthetic index that takes into account a number of indicators that represent social protection schemes that provide, where possible, protection against a single risk or need. For a comprehensive description of the social protection gross expenditure, readers can refer to the social protection report for the period 2012–2016, published in 2019 and directly

#### downloadable at: https://nso.gov.mt/en/publicatons/Publications\_by\_Unit/Documents/A2\_Public\_ Finance/Social%20Protection%202016.pdf

The social indicators taken into account in this study are:


#### **References**


© 2020 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* **Barrier Islands Resilience to Extreme Events: Do Earthquake and Tsunami Play a Role?**

**Ella Meilianda 1,2,\*, Franck Lavigne 3,4, Biswajeet Pradhan 5,6,7, Patrick Wassmer 4, Darusman Darusman <sup>8</sup> and Marjolein Dohmen-Janssen <sup>9</sup>**


**Abstract:** Barrier islands are indicators of coastal resilience. Previous studies have proven that barrier islands are surprisingly resilient to extreme storm events. At present, little is known about barrier systems' resilience to seismic events triggering tsunamis, co-seismic subsidence, and liquefaction. The objective of this study is, therefore, to investigate the morphological resilience of the barrier islands in responding to those secondary effects of seismic activity of the Sumatra–Andaman subduction zone and the Great Sumatran Fault system. Spatial analysis in Geographical Information Systems (GIS) was utilized to detect shoreline changes from the multi-source datasets of centennial time scale, including old topographic maps and satellite images from 1898 until 2017. Additionally, the earthquake and tsunami records and established conceptual models of storm effects to barrier systems, are corroborated to support possible forcing factors analysis. Two selected coastal sections possess different geomorphic settings are investigated: (1) Lambadeuk, the coast overlying the Sumatran Fault system, (2) Kuala Gigieng, located in between two segments of the Sumatran Fault System. Seven consecutive pairs of comparable old topographic maps and satellite images reveal remarkable morphological changes in the form of breaching, landward migrating, sinking, and complete disappearing in different periods of observation. While semi-protected embayed Lambadeuk is not resilient to repeated co-seismic land subsidence, the wave-dominated Kuala Gigieng coast is not resilient to the combination of tsunami and liquefaction events. The megatsunami triggered by the 2004 earthquake led to irreversible changes in the barrier islands on both coasts.

**Keywords:** Sumatra; barrier island; earthquake; tsunami; land subsidence; liquefaction; morphology resilience; GIS

#### **1. Introduction**

Major earthquakes that have occurred in history have caused fatalities and losses beyond the direct earthquake shaking. Looking at historical losses from 1900 onward,

**Citation:** Meilianda, E.; Lavigne, F.; Pradhan, B.; Wassmer, P.; Darusman, D.; Dohmen-Janssen, M. Barrier Islands Resilience to Extreme Events: Do Earthquake and Tsunami Play a Role? *Water* **2021**, *13*, 178. https:// doi.org/10.3390/w13020178

Received: 30 October 2020 Accepted: 8 January 2021 Published: 13 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

around 1% of damaging earthquakes in history have caused about 93% of deaths worldwide. Within those events, secondary effects, such as a tsunami, liquefaction, fire, and the impact of nuclear power plants, have caused 40% of economic losses and deaths [1]. Coastal areas situated at the peripheral of tectonic subduction zones, such as the Sunda Trench, Peru-Chile Trench, and Japan Trench, are subject to multiple threats of co-seismic and post-seismic secondary effects. The triggered massive landslides, liquefaction, and tsunami have caused severe casualties and damage to the affected areas. The Central Sulawesi tectonic earthquake in September 2018, for instance, has triggered multiple secondary effects, including localized tsunami, landslides, and liquefaction [2]. The Great Eastern Japan Earthquake in April 2011 triggered the tsunami, liquefaction, fire, and the impact of a nuclear power plant [1]. The Great Sumatran Earthquake in December 2004 triggered both mega-tsunamis and land subsidence, all of which have caused significant fatalities and losses and required long-term recovery.

Eighty percent of the world's coasts are rocky [3], characterized by a sedimentarydeficit, and a complex morphology [4]. On the basis of satellite inventory, Stutz and Pilkey [5] reveal that 20,783 km of shoreline are occupied by 2149 barrier islands worldwide, which is 37% longer than formerly thought. Indonesia is an archipelagic country with an estimated 91,363.65 km of coastline [6], making it the second longest coastline in the world after Canada [4]. The Indonesian archipelago has 44 barrier islands, which equals 2% of the global figure. Their distribution is strongly related to sea-level history in addition to the influence of tectonic settings. Administratively, 80% of districts and municipalities in Indonesia are situated at the coastal areas, making the coastline the multi-purpose zone in the daily life and socio-economic development of the archipelagic nation.

The shoreline along the Sumatra island of Indonesia mostly consists of shore-parallel lagoonal ecosystems separated by chenier-type barrier islands, with crest height on average no more than 2 m [7]. In the temperate zone, such as in the United States, barrier islands are under tremendous pressure of rising sea level and increasing storminess. On the other hand, the barrier islands and beach ridges located in a tectonically active coastal region, such as the entire length of the west and north Sumatra of Indonesia, are equally under tremendous pressure of the potential tsunami, co-seismic and post-seismic land level changing, and liquefaction as the secondary effects of the tectonic fault and subduction earthquakes.

Studies about the influence of the secondary effects of tectonic activities on the barrier island and spit morphological development are currently relatively limited. This study aims to investigate the morphological resilience of the barrier islands in responding to tsunamis, co-seismic subsidence, and liquefaction. Spatial analysis was utilized to detect shoreline changes from the multi-source datasets of the multi-decadal time scale, including old topographic maps and satellite images from 1898 until 2017. As the next step, we identified the history of major or moderate earthquakes that occurred in each consecutive period. Subsequently, we corroborated those earthquake events with the reported impact of the earthquakes on the coastal areas. We also utilized the established conceptual models of storm effects to barrier systems to interpret possible controlling factors causing the remarkable morphological changes observed at two selected coastal sections. The results will be substantial in building a conceptual foundation for interplaying seismic-related forcing factors with the littoral transport regime that determines the morphological resilience of a coastal area.

#### *1.1. Impact of Earthquakes on Coastal Systems*

Investigation of a massive earthquake event has always been appealing for scientists. It is the opportunity to improve the understanding of various factors resulting in casualties and damage. Short-term responses of coastal morphology to tsunamis have been widely discussed in the previous studies through various approaches. Some works have combined models of tsunami wave propagation and inundation, and spatial analysis of the damage in the tsunami-affected area [8–10]. Others focused on the geomorphological [11–14], or ecological [10,15] impacts of the tsunami. Additionally, a few tens of meters of land subsi-

dence were also observed as the immediate response of coastal areas to seismic activities was observed early after the occurrence of the megathrust earthquake of December 2004 along the north and west coast of Sumatra Island [16–19]. An extended period investigation to the post-tsunami coastal recovery in Aceh has also been conducted, including further monitoring on changes of shoreline position and land-level changes [20–22], of land use [23], and the vulnerability and living condition of the coastal inhabitants [8,24].

Liquefaction has yet to be discussed thoroughly, concerning the impact of the megathrust 2004 Sumatra earthquake in Banda Aceh. Nevertheless, some indication of liquefaction indeed observed in Port Blair, the capital city of the Andaman and Nicobar Island of India, where various types of buildings, particularly of the typical reinforced concrete ones, have experienced settlement [25]. Such building failures by settlement were also observed and reported in Banda Aceh [26]. Recently, the soil mechanism of the liquefaction potential in Banda Aceh city after the earthquake was evaluated by [27] using the semi-empirical Idriss Method to quantify the liquefaction potential. The result suggests that most of the highly built subdistricts in Banda Aceh are prone to liquefaction. It is noteworthy that the liquefaction impact is not only triggered by the shaking from the megathrust Sumatra–Andaman earthquakes but also potentially by the activation of the Great Sumatran Fault System [28].

Following the megathrust Sumatran earthquake and tsunami of 26 December 2004, a series of geodetic monitoring campaigns were conducted between 2005 and 2015, to monitor the development of the land-level changes of the west coast of Aceh since the 2004 tsunami [21]. The results showed that after experiencing abrupt co-seismic subsidence, the beach experienced uplift with a rate of 27 mm/year since late 2005. This number is an order of magnitude higher than the rate of eustatic sea-level rise, which is around 4–12 mm/year at the Indonesian waters [29]. The results demonstrate a possible relationship between the development of the new frontier beach ridge, the immediate co-seismic land subsidence, and the probable post-seismic rebound (uplift) associated with the viscoelastic mantle relaxation a few years following the mega earthquake.

Overall, the previous studies suggest that the secondary effects of seismic activity, either associated with the off-shore megathrust subduction or mainland tectonic faults activity, play a crucial role in altering the morphological development of the tectonically active coastal area.

#### *1.2. Barrier Islands and Spits Morphological Resilience*

Several studies of large barrier systems, including barrier island or beach ridges and sand spits, have been well-studied, primarily when associated with extreme storm events [30–35]. They have shown to constitute a valuable archive of coastal evolution and their morphology, and internal structures contain information on both past relative sea levels and past storminess activity. Barrier islands and sand spits are considered exemplars of coastal resilience [30,31] and critically essential ecosystems to protect the low-lying area behind them [36]. Under the influence of storms, large barrier systems, such as barrier island and beach ridges, are inherently resilient landforms as long as they can internally recycle sediment to maintain overall landform integrity [37], the rate of sea-level rise is not excessive, and there is no sediment deficit [38]. Nott et al. [33] suggest that the building of beach ridges is being slowed down by a decrease in sediment supply to the coast due to less frequent river floods during periods of reduced storminess. While during a storm event barrier islands may experience breaching due to the intensive funneling of the overwash flows into specific throats [39], a calmer wave regime promotes littoral drift, allowing sand spits to grow extensively [30].

Arguably, a tsunami event may contribute to the geomorphological imprint in the evolution of a tectonically active coast in a way similar to a major storm event. Despite different sources and mechanisms, the effect of the destruction caused by a tsunami may have been analog to that of a storm surge, in particular along the coastline, where barrier islands are the first line of natural coastal protection. Studies of short-term shoreline changes of barrier islands at Banda Aceh coast of Sumatra Island, Indonesia after the 2004 tsunami [13], and those at the Dauphin Island, US after the Hurricane Katrina in 2005 [40] are two comparable studies which demonstrated similar variability of shoreline changes caused by extreme events.

Despite the great diversity of storm mechanisms that strike coastal areas, we consider that the impact of the surging waves to barrier islands is equivalent to that of tsunamis. Sallenger [35] identifies four types of storm regimes and how barrier islands respond to those different forcing magnitudes, which later were re-described in detail by [34]. The four regimes are briefly described in [32] which consist of the "Swash regime", "Collision regime", "Overwash regime", and "Inundation regime," which involve the morphological effects of migration, escarpment, breaching, and submergence of the barrier islands, respectively. The case of sinking effects due to co-seismic subsidence or liquefaction is relatively easy to observe by identifying the narrowing subaerial parts by comparing a pair of maps or satellite images of consecutive years. However, it is almost impossible to determine the cause of the sinking by solely relying on the spatial analysis from the multi-temporal maps and satellite images. Therefore, we also corroborate the results from previous studies and the other sources of information and accountable reports to support the analysis.

Successful coastal management that includes mitigation of possible negative impacts must be based on an understanding of these patterns of change as natural responses to high-intensity events [41]. It is also essential to understand the coastal geomorphic setting and the geological boundaries before attempting to model the large-scale behavior of these types of coastal systems [42]. Herein, historical data inevitably play a major role in identifying any remarkable, even more so, the irreversible changes in the long-term past.

#### **2. Study Area**

We investigate morphological changes of the seaward-most barrier islands and spits along the coastline of the north tip of Sumatra Island, which is situated between 05◦16 15" N and 05◦36 16" N, and between 95◦16 15" E and 95◦22 35" E (Figure 1a,b). The low-lying coastal area behind them is where the capital city of Aceh Province, Banda Aceh situated in the central part, and surrounded by the Aceh Besar district occupying ca. 125 km2 northern valley of the Barisan mountain range at the north tip of Sumatra Island, Indonesia. The Barisan mountain range is formed as the backbone of the leading-edge Sumatra Island, parallel to the Great Sumatran Fault running parallel with the Sumatra-Andaman subduction zone, or also known as Sunda Trench (Figure 1a). The shoreline stretches ca. 25 km connecting Ujong Pancu headland in the southwest and the Ujong Batee headland in the northeast (Figure 1c). The brackish back-barrier wetland ecosystem, aquacultures, and lowland coastal villages occupy the area of ca. 4 km width from the coastline with elevations are varying of −0.5 m to +2.0 m from the mean sea level. The entire coastal area was severely devastated by the megathrust earthquake of M 9.0, followed by the gigantic tsunami event on 26 December 2004.

The tsunami disaster caused severe casualties and damaged properties, as well as erosion at the coastline, altering the functionality of the ecosystems and the livelihoods of coastal communities [44]. In Banda Aceh, the shoreline retreated as far as ca. 200 m inland, and erosion rates of 30 m3/m on average (up to 80 m3/m locally) were calculated [11]. Six months after the tsunami, the shoreline experienced about 15% further retreats from the initial erosion by the tsunami [13].

Major infrastructure, e.g., port basin, coastal revetments, and coastal roads, collapsed or were destroyed, whereas the coastal environment was also profoundly altered [8]. A comparison of a coastal revetment of before and after the 2004's tsunami at Lambadeuk is depicted in Figure 2a,b, respectively. The remnants of the same revetment were stranded offshore after the disaster event (Figure 2b), suggesting the occurrence of local land subsidence. Several preliminary studies suggest that the local land subsidence at Lambadeuk and Ulee Lheue (Figure 1c) was less than 50 cm [16,18]. On the other hand, the barrier islands which used to protect the built areas behind the Kuala Gigieng coast were breached at their weakest sections by the tsunami waves. Figure 2c,d exemplify one of the breaching

sections of the barrier island at Kuala Gigieng on the fourth day and six months after the tsunami, respectively.

**Figure 1.** Banda Aceh coast, at the northern tip of Sumatra Island, Indonesia: (**a**) Tectonic settings of Sumatra, compiled from various sources by Hurukawa et al. [43], showing the configuration of the Sumatran fault, slip rates, the plot of Global Centroid Moment tensor (CMT) solutions of shallow earthquakes (Mw ≥ 6.0 and depth ≤ 60 km) between 1976 and 2012; (**b**) Segments of the Great Sumatran Fault System at the northern Sumatra consists of two bifurcated fault segments, i.e., Aceh Segment and Seulimum Segment. The rectangle shows the investigated coastal area in front of Banda Aceh, the capital of Aceh Province, which is situated between the two fault segments; (**c**) The two investigated coastal sections in this study are Lambadeuk at the southwest and Kuala Gigieng at the northeast since they are the most dynamic and there is less interference by hard-structure coastal protection.

**Figure 2.** Indication of land subsidence at Lambadeuk, southwest of Banda Aceh coast: (**a**) The revetment located at the upper shoreface to protect the coast from further erosion nine months before the tsunami hit the coast on 26 December 2004; (**b**) The remnant of a seawall post-tsunami. The yellow arrows in both pictures show the same revetment; (**c**) The breaching of barrier island at Kuala Gigieng four days after the tsunami; (**d**) The development of sand spit growth six months after the tsunami, to reconnect the breached barrier islands at Kuala Gigieng.

At the northern Sumatra, the fault system splits into two active segments, i.e., Aceh Segment and Seulimum Segment (Figures 1b and 3a). Natawidjaja and Triyoso [28] found that currently, those segments possess a seismic gap of 325 km and 70 km long in the last 100 years, respectively, and considered them as an alarming hazard potential in the future. The highly populated Banda Aceh city is situated between these two segments (Figure 3a). The coast is facing the Andaman Sea and is semi-embayed by the entraining forearc small islands at the north off-shore from the rough, energetic waves of the Indian Ocean to the west. Lambadeuk exemplifies a relatively broad lowland coastal system overlying a major segment of the active tectonic fault, i.e., Aceh Segment. Kuala Gigieng displays a marine-dominated coastal system. The mainland consists of old parallel coastal ridges and swales. In the last few centuries and beyond, both coasts were naturally protected by the Late Holocene barrier islands as the coastal system's seaward-most land boundary. The alongshore multitemporal bathymetric profiles in Figure 3b illustrate the averagely shallow bathymetry in front of Kuala Gigieng, and the tilting southwest towards Lambadeuk, where considerably deep trenches are observed in two consecutive years in 1893 and 1924. In contrast, the bathymetry around the same location appears to have been remarkably shallow in 2006.

**Figure 3.** Banda Aceh is situated in the wedge between two tectonic fault segments of the Great Sumatran Fault System, i.e., the Aceh Segment on the southwest and Seulimum Segment on the northeast: (**a**) The map shows the topographic and bathymetric contours of the northern Sumatra region, and the two fault segments; (**b**) The multitemporal alongshore bathymetric profiles in front of Banda Aceh coast show higher elevation at Kuala Gigieng and lower down towards Lambadeuk. The bathymetric data of the year 1893 and 1924 was digitized from the Dutch Colonial Nautical Chart obtained from KITLV, and the 2006 bathymetry was obtained from the bathymetric survey conducted by UP-PSDA of Syiah Kuala University.

The tides along the coast are categorized into a micro-tidal regime which is on average less than one meter high. The typical equatorial monsoonal climate brings about seasonal prevailing wind-induced wave heights and periods variations, i.e., southwesterly during April to September and northeasterly during October to March [22]. The Aceh River is the primary natural river crossing the low-lying coastal city of Banda Aceh (Figure 1c). The river course has been artificially bifurcated at 10 km upstream to a 300 m wide artificial Alue Naga Floodway Canal since the early 1990s, another major outflow dissecting the coastline at present. Both major outlets are the primary sources of sediment supply to the coastal system of the investigated area. Currently, both are regulated by training jetties to maintain navigational depth for fishing boats. Another smaller river that flows at the northeast flank of Banda Aceh plain is the Angon River (Figure 1c), which has been nonmigratory throughout the last century, debouching into the shore-parallel lagoon behind the Kuala Gigieng inlet.

The coastal area was relatively densely populated prior to the 2004 tsunami event, occupied by approximately 250,000 people. Around half of the total population has been reported dead or gone missing [45] after the 2004 mega tsunami. More than a decade since the tsunami event, the coastal city has been rehabilitated and reconstructed. Despite the devastation due to the earthquake and tsunami, people are likely to return to their original living and business locations close to the coast. By 2009, the post-tsunami rehabilitation and reconstruction program in Banda Aceh had established resettlement at the coastal area up to 91% [46], along with all the necessary infrastructure such as a ferry port, religious and administrative buildings, schools, housing areas, and a sanitary landfill [23]. Apart from having been recovered from the tsunami event, the coastal area has been subjected to frequent flooding during high spring tides [47,48]. The coastal area is also vulnerable to hydrometeorological hazards, exceptionally high risk to the future sea-level rise [23]. We recently investigated the scenarios of coastal inundation due to the slow-onset projected sea-level rise in the next couple of centuries. The results show that the increasing number of built areas closer to the coastline, despite past tsunami experience, are potentially subject to tremendous loss due to seawater inundation in the next couple of centuries [23]. Such conditions provoke socio-economic and environmental vulnerability for the entire coastal area [49].

#### **3. Materials and Methods**

Historical records of particular events leading to morphological alteration are indispensable resources for better understanding the chronological implication to the state of coastal morphology. Previous studies have taken a similar approach by using the history of a few large-scale tectonic events to investigate the meso-term coastal changes (e.g., [50–52]). Herein, we consider that a multi-decadal timescale is appropriate to capture any remarkable changes in the coastal morphology induced by those extreme events associated with seismic activity, and also suitable for coastal management planning [53].

#### *3.1. Spatial Data*

We digitized shorelines from various data sources, i.e., from Colonial topographic maps of the 19th century and more recent satellite images of different spatial resolutions, the details of which are listed in Table 1.


**Table 1.** Spatial data sources for shoreline change detection.

All data sources were geo-referenced to a master map (i.e., ortho-rectified aerial photo acquired in June 2005 from the NORAD survey) processed in ArcGIS 10.6.1 to have a common horizontal datum, projection, and coordinate system. Figure 4 shows the samples of the multisource and multitemporal data used in this study.

#### *3.2. Historical Records of Earthquake and Tsunami Events*

To support our spatial analysis using the historical maps and satellite images, we corroborate the records and results of studies of the earthquake events in the Indian Ocean/Andaman Sea region and along the northern part of the Great Sumatran Fault System. The underlying active fault system along the Sumatra Island is one of the most significant active faults in the world, with slip rates ranges from 10 to 27 mm/year [28]. Natawidjaja and Triyoso [28] identified that the Great Sumatran Fault has at least 19 segments along the fault zone. More than a dozen large earthquakes have occurred historically in the past two centuries associated with the fault zone. Similarly, Hurukawa et al. [43] concluded in their study that almost all the earthquakes occurred since 1892 at Sumatra Island were located at Sumatran Fault, with high seismic hazard mainly occurred over the entire northern part of Sumatra Island during the period of 1942–2003.

Here, we limit the area of earthquake influence by setting up a regional boundary within which we identify any major (M ≥ 7.0) and moderate (6.5 ≤ M < 7.0) earthquake events records, as depicted in Figure 5. The historical earthquake records are obtained primarily from the NGDC/WDS Global Historical Tsunami Database of the National Centers for Environmental Information (NCEI) database [54]. Additionally, records were also obtained from the Catalogue of Significant and Destructive Earthquakes 1821–2017 for the earthquakes associated with the Great Sumatran Fault [55]. From the database we found five remarkable earthquake events which occurred within our investigation period, which are enlisted in Table 2. The locations of the earthquake epicenters are depicted in Figure 5.

**Figure 5.** Regional earthquake occurrences at the northern Sumatra during the investigation period in the present study (1898–2017). The dashed-line rectangle shows the regional boundary of influence of the major (M ≥ 7.0) and moderate (6.5 ≤ M < 7.0) earthquake records. Records of earthquakes M ≥ 7.0 accompanied with the validity of tsunami occurrences are retrieved from the earthquake catalog by [54], and the tectonic earthquakes associated with the Great Sumatran Fault obtained from [55].

**Table 2.** Earthquakes leading tsunami events affecting Banda Aceh in the period of 1804 to 2004. Geographical boundaries displayed in Figure 4.


Note: \* Data source from [54]; 4 = definite tsunami; 3 = probable tsunami; 2 = questionable tsunami; 1 = very doubtful tsunami; 0 = event like a seiche. \*\* Data source from [55].

> The NCEI database is accompanied by the reports of the past earthquakes and their associated resulting hazards (e.g., tsunamis, earthquake damage, liquefaction, etc.) that were mostly archived by Soloviev and Go [56] and also by several other scientific documents and papers. Here, the tsunami event records were gathered from scientific and scholarly sources, regional and worldwide catalogs, tide gauge reports, individual event reports, and unpublished works. Descriptions of tsunami events with various occurrence validity levels associated with some earthquake records were used in this study to identify the number of possible tsunami events in the last couple of centuries [57,58]. The level of validity ranges from 1 (less probable occurrence) to 4 (most probable occurrence). It also specifies the tsunami events that occurred at specific locations from where any tsunami waves would have propagated towards their surrounding coastlines. A complete discussion of possible errors can be found in [57].

#### *3.3. Storminess*

There is merely a little overview of the climatic condition in the last century in the investigated area. Verstappen [7] reported that, in general, the climatic condition at the northern tip of Sumatra island was relatively dry in the last century. The equatorial position of Indonesia shelters it from tropical cyclones that often devastate the Philippines, Sri Lanka, Bangladesh, and the coastal zone of western Australia. The recent studies [59–61] of the historical cyclones over the Bay of Bengal and the Andaman Sea reveals that the pathways of the cyclones revolve around the latitudes of higher than 8◦ N. Accordingly, the coastal areas at the northern tip of Sumatra (5◦ N) are assumed to experience much weaker tropical storm events. Thus, remarkable breaching events of barrier islands, such as those identified in the present study, are unlikely as the results of tropical storms, which is particularly important to keep in mind when analyzing the possible forcing factors responsible for observable morphological changes in this study.

#### *3.4. Regular Wave Climate and Littoral Transport Rate*

We obtained the wave climate of the Banda Aceh coast by converting the eleven-year daily wind data records from 1995 to 2005 from the National Meteorology, Climatology and Geophysics Agency (BMKG), and subsequently translated into statistical wave heights and periods as functions of wind velocity, duration, and fetch. Surface and tidal currents are not well-recorded in this coastal region. The tides moderately range of 1.00 m from MHWL to MLWL. From the data, we found that the prevailing monsoonal wave directions were from the northwest and the northeast. Based on the resulting wave data, we then roughly estimate the littoral transport rate by using the simple CERC formula [62] at both Lambadeuk and Kuala Gigieng coasts. Estimates of net littoral transport along both coasts will be discussed in Section 4.

#### *3.5. Analysis of Morphological Changes*

Morphological changes of the seaward most barrier islands and spits in seven consecutive periods were analyzed in this study at the Lambadeuk and Kuala Gigieng coasts. Each period consists of shorelines of the seaward most barrier islands and spits of two consecutive years. Any changes found to be significant were measured as indicative figures to estimate the amount of displacement or growth of the morphological features.

#### **4. Results**

There is a substantial difference in geomorphic settings between Lambadeuk and Kuala Gigieng coasts. Lambadeuk, at the southwest flank of Banda Aceh city, is a transgressive coast built up during the Late Holocene. Typically, a transgressive Holocene coastal barrier underwent continuous erosion and roll-over so that the oldest washover stratigraphies are likely to have been erased [32]. On the other hand, the Kuala Gigieng coast at the northeastern flank is a regressive coast, which typically consists of coast-parallel ridges and swales environment. One can also observe consecutive ridge-swale morphology preservation at the highly energetic and marine-predominant coastal areas along the western coasts of Aceh, Sumatra, e.g., in [21,22].

Figure 6 displays the consecutive observation periods in this study, along with the identified earthquake events (or the absence of those) that occurred in each period. Each of the barrier islands and spits at the Lambadeuk and Kuala Gigieng have uniquely evolved and altered its shape, position, and vulnerability to seismic-related hazards impacts. The most significant changes are evident from a sequential comparison of the island and spit geometries and rates of areal changes. Figure 7a–j depicts the overall results of the delineation of both coasts' morphological features, covering the shorelines of barrier islands, sand spits, and the back-barrier intertidal areas. The sub-aerial parts of the morphological features are color-coded, and the sub-aqueous parts appear in white. Here, we display the shoreline delineations of 1898 and 2005 in every figure panels to provide references of the morphological states at the initial date and the date a few months after the 2004 tsunami

event. The latter is represented by the shorelines delineated from the high-resolution aerial photograph acquired in June 2005 (Table 1). In between those years, we then observe and analyze the successive pairs of barrier shoreline changes.

**Figure 6.** The observation of morphological changes of the northern tip coast of Sumatra Island, consists of seven periods with various intervals depending on the data sets available for comparison. The horizontal axis shows the consecutive periods of observation from 1898 until 2017, and the vertical axis is the length of each successive period of observation. Periods with remarkable earthquake events are marked by the bold black line indicating the year the event occurred, which shows how far away from the pair of observation years the event is. The horizontal red line in 2000–2005 indicates 26 December 2004, marking the earthquake that triggered the mega-tsunami.

> We quantify the rate of change of the barrier islands' morphology by simply substituting the polygon areas of a pair of barriers' perimeters from two consecutive observation years, and dividing the result by the length of the period in between. To avoid misleading quantification of the rate of change, we exclude the quantification of the areal changes of the back-barrier lagoon and intertidal morphology. The misled quantification may come from sediment deposition to the back-barrier area from the barrier breaching and from the water inundation extend over the intertidal areas at the time the satellite image was acquired. Nevertheless, for the observation purposes, we still display the delineation of the morphological features at the back-barrier during each observation years, as was delineated from the datasets.

#### *4.1. Rate of Change of Barrier Islands Morphology*

From the wave data analysis, we obtained the seasonal prevailing monsoonal waves annually. The southwest monsoon occurs between April and September and is characterized by relatively rough waves coming from the northwest at the Banda Aceh coast. About 53% of the waves approach the coast with a significant wave height of 1.0 m with a period of 3 s. During the northeast monsoon between October and March, the climate tends to be milder, with 30% of the waves approaching from the northeast with a significant wave height of less than 1 m and a period of 4.5 s. Based on the resulting wave data, we then roughly estimate the littoral transport rate by using the simple CERC (CERC stands for Coastal Engineering Research Center, U.S. Army Engineer Waterways Experiment Station, Vicksburg, MS) formula [62], which results in the estimate of net littoral transports at a rate of +0.30 hectares/year and +1.08 hectares/year at Lambadeuk and Kuala Gigieng, respectively. The prevailing direction of net littoral sediment transport at both coasts is directed southwest, with positive rates suggesting the coasts are accretional.

Tables 3 and 4 display the rates of land gain or loss at Lambadeuk and Kuala Gigieng as the results of the barrier islands' morphological responses during the consecutive periods of observations. The largest amount of land loss is observed for all the barrier systems in 2000–2005 due to the impact of the 2004 earthquake and tsunami event, i.e., −73.54 hectare and −128.73 hectare at Lambadeuk and Kuala Gigieng, respectively. These are equivalent to a 100% and 59% loss of land with rates of −14.71 hectares/year and −25.75 hectares/year, respectively. Both rates are an order of magnitude higher than the estimate of the littoral transport rate under regular wave climate.

**Figure 7.** *Cont*.

**Figure 7.** Observation of barrier island morphological changes from the consecutive pairs of shoreline delineation of Lambadeuk and Kuala Gigieng coastal section, at the northern tip of Sumatra Island: (**a**,**b**) Changes at Lambadeuk and Kuala Gigieng, respectively, in two periods 1898–1924 (26 years), and 1924–1944 (20 years). Major earthquakes occurred in both periods; i.e., on 4 January 1907 of *M* 7.6 which is most likely to have triggered a tsunami event (tsunami validity = 4), and on 26 June 1941, which triggered a far-field tsunami event; (**c**,**d**) Changes at Lambadeuk and Kuala Gigieng, respectively, in 1944–1967 (23 years), during which a major earthquake occurred on 2 April 1964, of which the epicenter was relatively close to Kuala Gigieng, and reportedly triggered tsunami and liquefaction; (**e**,**f**) Changes at Lambadeuk and Kuala Gigieng in 1967–1989 (22 years), during which a moderate earthquake occurred on 4 April 1983 of *M* 6.9, and the epicentre was just off-shore from the Lambadeuk coast. No reports on any tsunami event. However, Lambadeuk reportedly experienced tectonic land subsidence; (**g**,**h**) Between 1989 and 2000 (11 years), no significant earthquake occurred at both Lambadeuk and Kuala Gigieng; (**i**,**j**). Dramatic changes occurred in the period of 2000–2005. During this short period, the megatsunami triggered by one of the largest magnitude earthquakes to have occurred in the modern history occurred on 26 December 2004, of M~9.0.




**Table 4.** Average long-term historical rates of land area change for both coasts for selected periods at Kuala Gigieng coast. Rates are in hectares/year. Positive numbers indicate land gain and negative numbers indicate land loss.

In the century before the mega-tsunami, Lambadeuk experienced significant land loss during 1944–1967 by −41%; nevertheless, it was compensated by 21% and 31% of land gain in the subsequent two consecutive periods (Table 3). Overall, Lambadeuk on average gained 5% extra land during the last century, only to face a complete loss of barrier island due to the 2004 tsunami until the present. In contrast, Lambadeuk experienced land loss of merely −1% averagely during the last century prior to the 2004 tsunami. The loss of barrier island area by −59% during the mega-tsunami has been compensated by 49% in 2017 since the tsunami.

#### *4.2. Multitemporal Morphological Changes*

The following is our analysis of the morphological changes of barrier islands and spits observable in Figure 7, and estimates of the total land loss or gain as the results of those changes based on the quantitative analysis described in Section 4.1. Following this analysis, the interpretation of possible controlling factors for those changes will be discussed in Section 5.

#### 4.2.1. Period 1898–1924–1944

In these consecutive periods, the main barrier island at Lambadeuk is in a relatively stable position in 1898–1924 and 1924–1944 (Figure 7a). Despite a major earthquake occurring in 1907, the barrier islands at Lambadeuk and Kuala Gigieng show a relatively stable state of morphology, which is observable by comparing the shorelines of 1898 and 1924 (Figure 7a,b). The growth of barrier spits extended from the eastern end of the barrier island at Lambadeuk suggests no extreme forcing, which could have caused remarkable changes in the barrier's morphology, such as landward migrating or breaching. A relatively remarkable barrier breaching occurred, however, at the adjacent coastal section at Ulee Lheue (the eastern barrier island next to Lambadeuk in Figure 7a). The possible cause is inconclusive in this study, as there was no report of any tsunami or land subsidence occurred as a result of the major earthquake, at the investigated area; despite a wellknown tsunami event which was reported at Simeulue Island at the southwestern offshore of Sumatra.

In 1924–1944, a barrier spit was growing further east as far as ca. 500 m, and almost connected to the remnant of a major breach of the Ulee Lheue barrier island (Figure 7a). In 1944, the spit was modified, along with slightly narrowing and counterclockwise reorientation of the barrier island's eastern part. Meanwhile, at Kuala Gigieng, the barrier islands (appear in blue in Figure 7b) had undergone significant changes. Two new inlets were created while the entire sub-aerial part of the barrier island has shifted landward at a maximum distance of ca. 300 m at the central part, and the section at the northeast at the same time moved seaward at more or less the same distance. Furthermore, this new barrier island formation seems to be narrowing entirely, resulting in a total land loss of −49.06 hectares with an average rate of about −2.45 hectares/year (Table 4).

#### 4.2.2. Period 1944–1967

In 1944–1967, the barrier islands and spits of both investigated coasts were significantly changing compared to the previous period. In 1967, both coasts experienced a relatively far landward shift, i.e., around 300 m distance, while the narrowing subaerial parts of the barrier features are remarkable (Figure 7c,d). Additionally, the barrier spits were breached at multiple locations at Kuala Gigieng (Figure 7d). A tectonic earthquake associated with the Seulimum segment activation occurred in 1964 with a magnitude M 6.7 [28,43], of which the epicenter was merely 73 km northeastern off-shore of Kuala Giging coast (Figure 5). The remarkable features of narrowing and breached barrier islands and spits remain observable in 1967, i.e., in merely three years after the earthquake event.

#### 4.2.3. Period 1967–1989

In 1967–1989, most of the barrier islands and spits on both coasts became relatively stable, i.e., no barrier island migration was found. Nevertheless, the westernmost part of the barrier at Lambadeuk that appeared in 1944 disappeared in 1967 (Figure 7e), suggesting submergence of the barrier island. In contrast, at Kuala Gigieng, the barrier islands were in a stable position throughout the whole period (Figure 7f). It is noteworthy that in this period, there was a moderate earthquake of magnitude M 6.9 which occurred in 1983, of which the epicenter was located just off-shore from Lambadeuk (Figure 5 and Table 2), suggesting that the earthquake was associated with the Aceh Segment which was underlying the coast.

#### 4.2.4. Period 1989–2000

There were no major or moderate earthquakes recorded in 1989–2000 (Table 2 and Figure 5). At Lambadeuk, the barrier island maintained its position during the 11 years, showing a stable and mature barrier spits at the eastern end. The previously submerged and breached barrier island in 1989 at the western end has been reconnected since (Figure 7g). At Kuala Gigieng, elongated barrier spits were developed, interestingly, in the opposite direction from the development that occurred during the last period. Although the barrier islands and spits position show no significant migration, a narrowing land area is apparent. The further growth of the right-hand barrier spit to the southeast appeared in 2000 (appeared in pink in Figure 7h) extending further southeast as far as 700 m from the tip of barrier spit in 1989.

#### 4.2.5. Period 2000–2005–2017

The years between 2000 and 2005 were the remarkable period where the great earthquake triggering the mega-tsunami occurred and severely destroyed the northern tip coast of Sumatra Island on 26 December 2004. The change of morphological features along the coast was expectedly enormous. Figure 7i,j display the massive overwash of the entire investigated coast leading to several breaches, and the remnants of the barrier islands were shifted landwards. The stabilized barrier spit which appeared in 2000 had disappeared entirely (Figure 7i,j).

Twelve years after the 2004 tsunami, the shoreline in 2017 reveals a striking difference in both of the investigated coasts in responding to the tsunami waves. At Lambadeuk, after the disappearance of the entire barrier island (Figure 7i), the shoreline started to develop from the former lagoon's inner shore, which has since been exposed to the sea. Kuala Gigieng now has two relatively permanently open inlets because of the barrier island and barrier spit breaching during the 2004 tsunami (Figure 7j).

#### **5. Discussion**

Our investigation on the morphological changes of barrier islands and spits at two coastal settings along the northern tip of Sumatra Island revealed several new findings on the possible controlling factors associated with the seismic activity that change barrier island morphology. The intervals between periods of observation of morphological

changes are of less than 30 years, which fits the return period of major earthquakes in the region of Sumatra. The impact on the coastal morphology can be in the form of barriers breaching, sinking, landward migrating, and landward bending of spits development. The possible interplaying forcing factors to those morphological changes are discussed in the following subsections.

#### *5.1. Tsunami Overwash*

In 1944, the barrier island at Lambadeuk remained stable compared to its morphological state in 1924, with the additional growth of barrier spit (Figure 7a). In contrast, the barrier islands at Kuala Gigieng experienced a landward migration, breaching, and the growth of barrier spit bending landwards (Figure 7b). Such morphological pattern is typical for barrier islands washed over by storm events or hurricanes, e.g., Hurricane Sandy and Hurricane Katrina [40]. Energy dissipation is then mostly achieved through the overwash bore running over the barrier, flattening the barrier crest profile, and depositing off-shore originating material over the back-barrier zone [32]. Provided that the height of the barrier island was typically less than 2 m high [7], the forcing factor which possibly controls such morphological impact is a wave coming from the sea. Waves equal to or slightly higher than the height of a barrier island, with a fairly long period, may have overtopped some of the lowest points of the barrier island. Analogously to a storm event, among the four regime types of impacts caused by storm waves on a barrier island [32], the observable morphological changes at Kuala Gigieng in 1944 may fall into the "Overwash regime".

In 1924–1944, there were no earthquake events associated with either the Sunda Trench and the Great Sumatran Fault system that occurred near the northern tip of Sumatra Island in this period. However, a major earthquake of M 7.6 occurred in 1941 at the western coast of Car Nicobar Island, with its epicenter at 12.1◦ N and 92.5◦ E, or around 834 km northwest from Banda Aceh (see location in Figure 5). NCEI [54] recorded a high score validity tsunami event associated with this earthquake (Table 2). From the historical reports [63,64], the earthquake generated a tsunami throughout the Andaman Sea and the Indian Ocean. Despite the lack of reliable records, for instance from the tidal gauge, which was in operation at the time, some local newspapers in India reported that a tsunami was witnessed along the eastern coast of India. At the Nicobar and Andaman Islands, the wave heights were reported as high as 0.75–1.75 m. It was estimated that 3000 to 5000 people were killed in Sri Lanka and on the east coast of India [65]. Local newspapers believed to have mistaken the reported deaths and damage to a storm surge.

In support of these reports, a numerical model developed by [66] was to simulate the Andaman–Sumatra tsunami propagation of tsunami triggered by the earthquake that occurred in 1941. The result reveals an agreement with the documented observation reports that the earthquake felt over a wide area covering the eastern and southern Andaman Sea region (i.e., the northern coast of Sumatra). The model results show that the tsunami has reached the coast of Nagapittinam on the west coast of India after 165 min, with run-up heights in the range of 0.95–1.25 m, and the north tip of Sumatra at 120 with tsunami wave height was less than 1.00 m [66].

It is noteworthy to mention that following the earthquake, two events with magnitude 6.0 struck within 24 h after the main shock of 27 June 1941, and there were 14 aftershocks of magnitude up to 6.0 until January 1942 [67]. At Lambadeuk coast, McKinnon [68] reported that a local informant indicated that the location of the village of Lambaro, located at the western end of the barrier island of Lambadeuk (Figure 7a), had been moved three times within living memory, the last time being early in the Japanese occupation (1941–1942) when the inhabitants were forcibly removed from the shoreline and resettled at the foot of the surrounding hills. Thus, the far-field tsunami event in 1941 was most likely responsible for changes of the barrier islands and spits at Kuala Gigieng in 1924–1944. The nearshore bathymetry may have been controlling the tsunami arrival at the shoreline of Lambadeuk and Kuala Gigieng. The deepening bathymetry (Figure 3b) may help to dampen the

incoming tsunami wave energy at Lambadeuk, while the shallow bathymetry at Kuala Gigieng may have amplified the tsunami wave force.

#### *5.2. Combined Liquefaction and Tsunami Overwash*

The observation period of 1944–1967 exemplifies the period where the impact of both liquefaction and tsunami occurred, and the effects are observable on both coasts. The barrier islands and spits delineated from the satellite image in 1967 in Figure 7c,d show the state of the coast three years after a major earthquake event on 2 April 1964 of magnitude *M* 7.0, which epicenter located at 5.6◦ N and 95.4◦ E (see Table 2), or ca. 12 km northeastern off-shore of Banda Aceh City. The barrier islands and spits at both coasts experienced hundreds of meters of landward migration (and slightly reoriented clockwise at Lambadeuk), multiple-barrier breaching, especially at Kuala Gigieng coast, as well as a narrowing of subaerial parts of barriers, most probably as a result of land subsidence induced by liquefaction (Figure 7b,c). Surprisingly, the multiple inlets as the results of breaching events remained open even after three years of development, suggesting that littoral transport and sediment supply under the regular wave climate have not made up for the loss of land caused by the tsunami and liquefaction.

Soloviev [56] reported that an earthquake that occurred on 2 April 1964 had caused considerable damage to adobe buildings, the ground cracked open, the ground subsided, mud and sand gryphon appeared. Such impacts of ground shaking are comparable to our observation on the 2016 Pidie Jaya tectonic fault earthquake of *M* 6.5, in Pidie Jaya district on the northern coast of Aceh Province [69]. In Pidie Jaya earthquake, black sandy mud emerged through small cracks opening in the internal floors of houses. Both coasts are comparable for their similar type of wave-dominated sandy coastal area with typical soil type of alluvium containing gravels, sands, and muds [70]. Such soil structure provides some degree of tremor amplification [69], which creates saturation on sandy layers and increases pore water pressure, leading to liquefaction [71].

NCEI [54] recorded a high-validity event of the tsunami (validity 3), mostly based on a report by [56], in which a combination of a tidal (tsunami) wave and the locals observed land subsidence of half a meter at Ulee Lheue, east of Lambadeuk barrier island (Figure 7c). The multiple breachings that appeared at Kuala Gigieng in 1967 (Figure 7d) may have been further submerged by the already breached barrier islands in the earlier period (1924–1944), suggesting an effect of liquefaction combined with a tsunami overwash.

#### *5.3. Co-Seismic Tectonic Subsidence*

The observation period of 1967–1989 reveals a unique contrasting development of barrier system morphology between the two observed coasts. Contrary to the long-distance growth of sand spits to their distal length observed at Kuala Gigieng (Figure 7f), the western side of barrier island at Lambadeuk (Figure 7e) appeared to have been heavily subsided, that the slightly more than 1 km barrier island connected to the Ujong Pancu headland which appeared in 1967 and had disappeared entirely in the satellite image of 1989. This left a piece of subaerial part detached from the main barrier island. Such contrasting changes most likely have something to do with tectonic activities within the period.

A moderate tectonic earthquake occurred on 4 April 1983 with magnitude *M* 6.9 [72], and the epicenter was located at 70 km northwest off-shore of Banda Aceh (Figure 5, Table 2). The earthquake caused significant damage to buildings in Banda Aceh city, with the Mercalli scale recorded as category VI [55]. There was no report from the locals of any tsunami occurrence or liquefactions such as those which occurred in 1964. It is also noteworthy that the submerged area is close to the underlying Aceh Segment, suggesting tectonic subsidence responsible for the submergence. McKinnon [73] investigated the archaeological artifacts at this location, which suggests localized tectonic subsidence in the vicinity of the Sumatran Fault System, which appeared to be dramatically evident around Lambadeuk. Strong evidence of submergence came from the rectangular structure discernible during the low tides, which happened to be a former mosque foundation that

remained visible underneath the water in an aerial photograph from 1978. An interview with a local informant suggested that in the last 80 years towards the first Japanese occupation in 1942, a sunk off the coast at least two to three meters had been occurring [73]. The shoreline had retreated about 150–200 m at a village called Lambaro, which was originally located at the Lambadeuk barrier island. This village, which had already been moved three times, was eventually rebuilt about 12 km inland at present. Bearing in mind that this area has evidently been submerged in two consecutive periods, i.e., 1944–1967 and 1967–1989, this suggests that Lambadeuk is not resilient to multiple tectonic subsidence.

#### *5.4. Dormant Period of Major and Moderate Earthquake*

Among the entire periods of our investigation, there was also a period where major or moderate earthquakes were absent, which was during 1989–2000 (Figure 7g,h). Based on the work of Nott et al. [33], one might have expected that less intense storm periods would have induced a slowing down of the development of beach ridges (in this case, the barrier islands). Instead, we note that the reconnection of the previously breached barrier islands occurred during this period (Figure 7g,h). Moreover, the further growth of the barrier spit at Kuala Gigieng towards the opposite direction from the previous period (Figure 7h) suggests that the prolonged net littoral drift prevailing southwestern direction is most likely responsible for the observed changing course spit growth.

Despite the stable development of barrier islands morphology during the non-event periods by the littoral process, the development of barrier islands depends on the continuous sediment supply from the rivers. A consistent prevailing longshore drift can promote such stable elongated barrier spit growth in a considerably long period, without any remarkable disruption such as by tectonic events [74]. A continuous sediment supply from the major outlets may have facilitated the growth of long barrier spits to their distal lengths, whereas the seasonal change of littoral drift between the alternating rough west monsoon and calm northeast monsoon may promote the spits' balanced growth.

#### *5.5. Mega Tsunami Overwash*

The great earthquake of 26 December 2004 has generated not only a tsunami but also land subsidence along the northern tip coast of Sumatra, particularly along its western coast. There were only a few meters of land subsidence that occurred in Banda Aceh [16,18]. The locals also observed liquefaction effects of clayed soils in many areas at Banda Aceh city [8]. Despite these, the far-field gigantic tsunami waves of more than 10 meters in height arrived at the northern tip of Sumatra and were most definitely the controlling factors of the remarkable changes of the barrier islands at Lambadeuk and Kuala Gigieng. This includes the extensive disappearance of the barrier islands (Figure 7i,j), which had long been preserved as the natural coastal protection to the back-barrier ecosystem since the Late Holocene.

Even after 14 years of post-tsunami develpment, the coastal area remains exposed to the ocean with the absence of a new barrier island. Clearly, the 2004 tsunami event has been the primary forcing factor responsible for the disappearance of most parts of the barrier islands at Lambadeuk and Kuala Gigieng (Figure 7g,h). For the case of a tsunami, out of the four types of storm regimes and prediction of beach changes proposed by [35], the presumable tsunami overwash that occurred in the 2004 event most probably falls into regime 4: "Inundation regime". Here, the barrier islands were overtopped by the tsunami, of which the height is considerably higher than the crest of the barrier island, overtopping the barrier island and flattening the barrier topography.

#### *5.6. Fluvial and Lagoon Systems*

The three major outlets supplying sediments to the coastal area are the Aceh River, Alue Naga Floodway Canal, and Angon River. All of them have been generally stable in their position throughout the entire investigated period, even after the 2004 tsunami. Small river streams flowing into the lagoon system have not been migrating from their original

position during the period of investigation, particularly those the locations of which were not directly associated with the tectonic fault (e.g., Sumatran Fault) and the length of the river streams are relatively short, i.e., less than 2 km. Combined with wave actions, in particular at the wave-dominated Kuala Gigieng coastal system, the sediment supply may have been actively contributed to the relatively quick reconnection of the barrier island breaching and the overall recovery of the associated barrier islands.

Verstappen [7] characterized the Sumatran river system, whereby the river discharge was affected by the monsoon. The sand fraction of the material carried off by the rivers and reaching the sea usually is considerably small as a result of the intense chemical weathering inherent in the prevailing humid tropical environment. Moreover, in the locations where neo-tectonism is dominant, the prevailing climatic conditions are less important compared to those at the more seismically stable locations.

#### *5.7. The Morphological Resilience of Tectonically-Active Coasts*

The present study reveals that major and moderate earthquakes with a return period of 20 to 30 years at the north tip of Sumatra Island of Indonesia have the potential to trigger tsunamis and co-seismic land subsidence and liquefaction. Figure 8 summarizes the varying land areas of barrier islands from one period to another subject to the presence of major or moderate earthquakes or the absence of them at Lambadeuk and Kuala Gigieng. In the previous discussions, we have been able to identify several different cases of morphological changes of barrier islands, which were driven by various controlling factors associated with major and moderate earthquake events in the last century.

**Figure 8.** Historical land area trends for the Lambadeuk and Kuala Gigieng barrier islands to the timing of major earthquakes that impacted the islands.

The mega-tsunami of 2004 remarks the fundamental changes in terms of coastal morphological characteristics, i.e., from a low-lying wetland ecosystem naturally protected by barrier islands and spits into a wave-exposed intertidal coastal area. The sediment transport under the present day's wave regime is unlikely capable of restoring the removal of those barrier islands that was built during the Late Holocene, most probably under an extremely different wave regime. Thus, we may conclude that both Lambadeuk and Kuala Gigieng are not resilient to a mega-tsunami.

For the case of barrier islands directly associated with the tectonic fault system, a combination of co-seismic subsidence, liquefaction, and tsunami can dramatically change the barrier island morphology, such as Lambadeuk. As for barrier islands located away from the tectonic fault system, such as Kuala Gigieng, the tsunami and liquefaction are responsible for modifying the barrier morphology, either independently or simultaneously. The breaching, landward migration, and the sinking barrier islands which occurred in the consecutive periods have been proven irreversible, despite a temporary growth of extended barrier spits at a sediment-rich environment. The unrecovered breaching, the submergence, as well as the disappearance of barriers observed even before the event of mega-tsunami of 2004, suggest that the barrier islands are not resilient enough against the recurrence of seismic events.

Evidently, within a century, the recurrence of several major earthquakes, either associated with subduction or tectonic fault zones, may trigger secondary effects, which potentially reduce the functionality of barrier islands and spits as natural coastal protection. The recurrence of small tsunamis, tectonic land subsidence, and liquefaction falls within an engineering timescale, which, therefore, should ultimately be taken into account in managing such a tectonically active coastal settings.

#### **6. Conclusions**

The present study investigates the resilience of barrier islands and spits because they protect low-lying coastal areas with high economic and environmental value. This is particularly true in archipelagic countries like Indonesia. Spatial analysis of barrier islands and spits is delineated in GIS, utilizing a multitemporal and multisource data series from the old Colonial topographic maps to the most recent high-resolution satellite images of 1898 to 2017, which encompasses a multidecade timescale observation. The earthquake and tsunami records and established conceptual models of storm effects to barrier systems are corroborated to support possible forcing factor interpretations. We found that tsunami, co-seismic subsidence, and liquefaction are the secondary effects of moderate or major earthquake occurrences that mostly control the state of the modern barrier islands and spits morphology in the investigated coastal area. Those controlling factors intermittently disrupt and interplay with sediment transport otherwise induced by regular wave climate, and eventually alter the coastal morphology development trend in the long term. The records of earthquake occurrences in front of the Banda Aceh coast suggest a return period of 20 to 30 years of major and moderate earthquakes, most of which triggered tsunami events (e.g., earthquakes in 1907, 1941, and 1964), and have far weaker wave energy than the one that occurred in December 2004. Evidently, liquefaction and co-seismic subsidence are the other controlling factors contributing to the remarkable morphological changes, which in some cases may be coupled with tsunami events. The results demonstrate that the semiprotected embayed Lambadeuk coast has been progressively lost its barrier islands due to repeated co-seismic subsidence. The wave-dominated Kuala Gigieng coast is not resilient to the combination of tsunami and liquefaction events. The mega-tsunami triggered by the 2004 earthquake has led to irreversible changes in both coasts' barrier systems. The barrier system is supposed to be a natural protection measure for the ecosystem and the economic value of the protected back-barrier domain. The exposed coastal mainland due to the barrier system disappearance may increase the risk of disaster by tsunamis, co-seismic subsidence, and liquefaction in the future. A further comprehensive evaluation of the resilience of the tectonically active coast, therefore, needs to be done. For such a dynamic tectonically active coastal area, investment in coastal protection to protect the invaluable ecosystem, and economic activities on the back-barrier domain will, therefore, be challenging. The results of the present study imply that in managing the coastal area where seismicity is highly active, despite small magnitudes, tsunamis, liquefaction, and land subsidence should be considered well.

**Author Contributions:** Conceptualization, E.M., F.L. and B.P.; methodology, E.M., F.L., M.D.-J.; software, E.M., F.L.; validation, E.M., F.L., P.W.; formal analysis, E.M.; investigation, E.M., P.W.; resources, E.M., D.D.; data curation, E.M.; writing—original draft preparation, E.M.; writing—review and editing, B.P., F.L., M.D.-J.; visualization, E.M.; supervision, M.D.-J., F.L.; project administration, D.D.; funding acquisition, E.M., D.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Indonesian Ministry of Research, Technology and Higher Education (Ristekdikti) through Scheme for Academic Mobility and Exchange (SAME) Program and World Class Professor Program (WCP) Scheme A.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

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

#### **References**


#### *Article*
