Incorporating dynamics of land use and land cover changes into tsunami numerical modelling for future tsunamis in Banda Aceh

This study aims to investigate the tsunami exposure to current land use land cover (LULC) with the LULC predictions for the next 50 years, including the impact of sea-level rise (SLR) in Banda Aceh. This study began with a spatio-temporal dynamic analysis using LULC changes projection. Firstly, Markov Chain was used to simulate the changing trends in land use from 2011 to 2015. The changing trends were used to simulate land use predictions for 2015-2065 using Cellular Automata (CA). There were six main types of LULC classifications, i.e., ponds, built-up areas, mangroves, bare land, urban forests, and water areas. Furthermore, the map resulting from the land use simulation was used as the Manning Coefficients for tsunami simulation using the COMCOT model. The initial tsunami waves were generated based on the 2004 Indian Ocean Tsunami caused by 9.15 Mw earthquake. As a dynamic tsunami hazard approach, a sealevel rise of 0.35 m was considered for the next 50 years. The study results show the built-up area will be affected by the tsunami, about two-thirds of the city's built area. Growth tends to occur in the southern part of the city outside the tsunami hazard zone. But there will also be the growth of built-up areas near the coast. Although much development is observed in the southern part of the city, the coastal area also needs attention because the sea level rise factor can exacerbate the impacts of tsunamis in the future.


Introduction
Population in coastal areas grows rapidly from time to time. The primary effect of high population growth is a change in land use and land cover (LULC). Population growth in coastal areas is followed by potential associated coastal hazards. The threat of hazards facing coastal areas becomes more concerning due to the impacts of global warming inducing sea level rise (SLR) [1]. Several studies have proven that SLR can cause flooding [2] and exacerbate the impact of tsunami waves [3], [4].
The 2004 Indian Ocean tsunami caused massive destructions on coastal cities in several countries in the Indian Ocean basin. The worst affected area was Banda Aceh in Indonesia, with around 165,000 people killed. At the time, the tsunami inundated 60% of the city of Banda Aceh, destroying everything in its path due to the extreme wave velocity. Coastal sub-districts directly adjacent to the beach experienced complete damage by the tsunami [5]. During the rehabilitation and reconstruction period between 2005-2009, Banda Aceh underwent massive reconstruction with 90.8% growth in the built-up area [6]. As the capital city of Aceh Province, population growth is mainly an urbanization factor. Urbanization commonly occurs in the main city like Banda Aceh because it is the center of government, economy, and education.
A tsunami is a significant disaster event with low frequent but high impacts events. Imamura et al. [7] * Corresponding author: syamsidik@tdmrc.org explain that future disaster impacts are increasingly complex. Therefore, multi-disciplinary researches are needed to provide comprehensive and sustainable disaster risk reductions [7]. The phenomenon of the tsunami threat is not a static event. There are other coastal hazards that will contribute to more significant impacts on the community. On the other hand, the element at risk is also dynamic along with the changes in the population. Preparing disaster mitigation is only based on static historical data; without taking into account the dynamic factors. However, such an approach provides a limited understanding of longerterm tsunami risk. For this reason, it is necessary to incorporate dynamic hazard components and the dynamics of the elements at risk.
In a previous study, Sleeter et al. [8] observed the human exposure in 2061 by the tsunami with predicting land use in the US-North Pacific without considering the SLR. Another study examined the impact of flooding due to SLR with projected land use in 2030 and 2080 in Bay County, North Florida [9]. Up to the present, only limited studies could be identified evaluating impacts of the future tsunami based on dynamic LULC changes.
The objective of this research is to investigate the impact of land-use exposure in 2065 in tsunami-affected areas in Banda Aceh. The 2065 projection was chosen to highlight the effect of sea-level rise in the next 50year. The year 2065 has also entered the 3rd period of Banda Aceh's spatial planning regulation after the tsunami recovery. The affected area was analyzed based on the city administrative boundaries. The results of this research are expected to fill the gap from long-term tsunami mitigation efforts, especially for disaster risk reduction-based spatial planning in Banda Aceh. Fig. 1. Banda Aceh's study area (shown by the blue box). This map is visualized from World Imagery.

Study area
As illustrated in Fig. 1, Banda Aceh is located at the northern extremity of the Sumatra mainland. Banda Aceh covers 60 km 2 and has a 15-km shoreline. Banda Aceh city is flanked on the south and east by the Aceh Besar district. Banda Aceh is exposed to coastal hazards due to its flat inland geography and average height of 0.8 meters above sea level. Based on the study of 100years SLR projection [10], Banda Aceh will experience coastal flooding and retreat of the shoreline up to 98 m. The city also is situated closely to the Sumatra-Andaman subduction zone, which can produce a large tsunami-like event similar to the 2004 Indian Ocean tsunami. If a tsunami occurs under the sea level rise situation, the tsunami will spread further and exacerbate the impact of human exposure more significantly than the 2004 tsunami [4].
Demographically, the study area includes urban and rural settlements. Meuraxa, Kuta Raja, Kuta Alam, and Syiah Kuala are the four seaside sub-districts of the city. Two of them -Meuraxa and Kuta Raja -were damaged by the 2004 tsunami. As the capital city of Aceh province, Banda Aceh is the center of government and the economy with the highest population density compared to other districts in Aceh. A study revealed that the community preparedness of Banda Aceh toward coastal hazards is a low category. The low category also includes an overall assessment of the city's resilience factors from coastal hazards [11].

Method
This research has three main steps: LULC change modeling, tsunami modeling, and land use exposure assessment. The first step is land use projection in 2065 using Markov-Chain and Cellular Automata [12] based on historical data, LULC 2011 and 2015. The second step is tsunami simulation using Cornell Multi Grid Coupled Tsunami Model (COMCOT) with tsunami wave generation from the 2004 earthquake with 9.15 Mw. The tsunami simulation was run on two scenarios, current condition without sea-level rise with LULC 2015 and future condition with 0.35 m of SLR 50 years with LULC 2065. The LULC maps are used for Manning's roughness coefficient on tsunami simulation. By overlaying the tsunami inundation map and the LULC map, we estimate the land use exposure in five land use classes.

Land use land cover classification
The study area's land use map was created using data from Bappeda (Banda Aceh Development and Planning Agency) for the years 2011 and 2015. These maps were interpreted from visual image analysis produced by Worldview-3 image. Water bodies, built-up areas, bare land, fish ponds, city forests, and mangrove forests were grouped into six key land use types for this study. This land cover classification was performed using ArcGIS software, adjusting to the administrative boundaries of Banda Aceh city. The land category and the descriptions of each land use can be seen in Table 1.  Fig. 2. To validate the land use simulation, an actual 2019 map was used. This map was digitized from the SPOT 5 satellite image. The current scenario (LULC 2015) and the future (LULC 2065) were used as manning roughness coefficient maps for tsunami simulations. The classification of Manning coefficient are waterbody (0.013), fishpond (0.017), bare land (0.02), mangrove (0.035), built-up area (0.04), and city forest (0.035). Land use is assigned to each grid cell in the research region, and Manning's coefficient is assigned to each land use.

Land use modelling
To forecast future land use in 2065, a Markov Chain and Cellular Automata model was used. This model is already part of IDRISI Selva software that is capable of land-use simulation [14]. For this study, the model was developed to produce a spatial projection of change based on historical data (2011-2015). The information was collected from Bappeda Banda Aceh The selection of data for that year was because the tsunami recovery process was considered to have returned to normal condition so that the trend of changing land use can be used to predict land use in the next 50 years. This data was used to determine the transition matrix of the two maps. The use of two actual maps in projection land use has been done by Mansour et al. [15]. Their projection map shows the good agreement of validation. The overall land use projection flowchart can be seen in Fig.  3.

Markov Chain
The Markov chain model analyzes change by calculating the probability transition area over a period of time. A probability matrix was applied to predict land use [16]. The Markov chain is a stochastic model that uses the trend of LULC changes in the preceding year to estimate the probability of land-use change. The model in this research was based on transition probabilities. From 2011 to 2015, the LULC's transition probability matrix, which was used to forecast land use, changes. [17].
Markov Chain uses the following mathematical equations [18]: is a transition matrix from land-use i to land-use j, where ( +1) and ( ) are land-use at time t+1 and t, respectively. And i and j are land use at first and second times.

Celular Automata-Markov (CA-Markov)
Simulating dynamic spatial phenomena and predicting future land-use change has been done using the CA-Markov model. The combined CA-Markov chain model also takes advantage of the Markov chain expected amount of land-use change and the CA models dynamic explicit spatial simulation. As a result, the CA-Markov model can be used to model spatial land-use change [18].
The equation used in the Cellular Automata module is as follows Where S is the limit of limited and discrete cellular states, N denotes the cellular field, t and t +1 denote the time difference, and f denotes the transformation rule. In

Tsunami numerical modelling
Tsunami numerical modeling was run using the Cornell Multi-grid Coupled Tsunami Model (COMCOT). The tsunami simulation includes tsunami generation from the earthquake source to the tsunami propagation and inundation. Some studies have been successfully demonstrated tsunami propagation along the Sumatrasubduction zone using COMCOT [19] [20]. In this study, the earthquake source is data from the 2004 earthquake with an earthquake magnitude of Mw 9.15 with a thrust fault earthquake mechanism. The source of the earthquake as the center of seabed deformation is reproduced by a fault model that provides information on earthquake parameters (Table. 2). The earthquake parameters, including magnitude, depth, fault length, fault width, epicenter coordinates, and angle of movement, have been developed by Koshimura et al. [21]. These have been validated for tsunami heights on land. With a total fault length of 1,155 km and a depth of 10 km, this fault model stretches from the Andamans to the northern half of the island of Simeulu.   Land use maps of the current and future conditions assigned the Manning coefficient are used as roughness variables in the tsunami numerical simulation using COMCOT. The Manning coefficient does not apply in Layers 1 to 5, and the various Manning coefficient is only used on the innermost layer (Layer 6). Tsunami waves were simulated using an earthquake magnitude Mw 9.15 using six simulation layers with a grid size, 8.5 m x 8.5 m on the smallest layer. The grid settings as shown in Table 3, and the domain computations are shown in Fig. 4.
COMCOT is a two-dimensional horizontal model that uses shallow water equations (SWE) to simulate depth-averaged velocities [22,23]. Leapfrog methods are used to solve the equations. SWEs were switched between linear and nonlinear in the model. The linear SWEs in a spherical coordinate can be seen in Equations (4) to (6).
Equations (7) to (13) The bottom friction terms using Manning's formula are described as follows: Where  is the water surface elevation; (P, Q) denotes the volume fluxes in X (West-East) and Y (South-North) directions, respectively; ( , ) denote the latitude and longitude of the Earth; The Earth's radius is R; the gravitational acceleration is g, and h is the water depth. f is the Coriolis parameter, and Ω is the Earth's angular velocity.

Land use land cover exposure
Land use land cover (LULC) current condition and future prediction data are overlaid with a simulated tsunami inundation map with an earthquake magnitude of Mw 9.15. The map is specified in the coordinates WGS1984-UTM zone-46N. The LULC in the tsunamiaffected area was calculated for each class of land cover. The total affected area and the safe area are also calculated spatially. This process is done with ArcGIS. Furthermore, the results of the current and the future conditions are compared.  The 2065 LULC shows that the built-up area of Banda Aceh City will increase significantly. Compared to the 2015 LULC, the bare land area has changed a lot to build-up areas in the central and southern parts of the city. The simulation results demonstrate a quantitative change in the area (Fig. 6) in each land cover type. The largest land cover is the build-up area, with an area of more than 32 km 2 . Tsunami wave propagation was simulated by a numerical model using COMCOT. The initial water level of the 2004 tsunami is shown in Fig. 7. This change in water level occurs due to differences in deformation along the fault zone. Blue color indicates a decrease in water level elevation, while an increase in water level elevation is indicated by red. The decrease in elevation in the western part of the fault close to the west coast of Aceh caused the seawater to retreat shortly after the earthquake. Then, 30 minutes after the earthquake, the tsunami waves came to the coast of Banda Aceh and inundated the 3 km coastal area of the city. In the 2065 future scenarios, the tsunami totally flooded five subdistricts: Meuraxa, Kutaraja, Kuta Alam, Syiah Kuala, and Jaya Baru, as shown in Fig. 8.

Tsunami exposure on LULC
Simulations were run for both inundation scenarios, including current conditions in 2015 without SLR and future conditions in 2065 with the 50-year SLR incorporated in, as shown in Fig 9. Under current conditions, the inundated area was 41.68 km 2 or 69.5% of the total city area. In the future scenario in 2065, the inundated area increased to 44.70 km 2 or 74.5%. Only 15.30 km 2 of a safe area was left in the southern part of the city. The overlaying of these two scenarios indicates that all of the LULC classes were affected by the tsunami. These conditions occur because the various LULC in Banda Aceh are more concentrated in coastal areas than the southern part of the city. Figure 10 depicts the land-use area that was flooded by the tsunami. The fish ponds area, in both scenarios, were all affected by the tsunami. In the 2004 tsunami event, fish ponds took a significantly long time to recover to be used again by fishermen. It means that fishermen's livelihoods were disrupted for a long time. The land cover most affected by the tsunami in the scenario of 2065 was the built-up area with more than 20 km 2 or ½ of the tsunami-affected site. The built-up area here includes residential areas, public buildings, roads, offices, and economic regions. The larger the residential area affected, the larger the number of the affected population. Other vital assets located in the affected built-up area will result in economic losses of the city. This study examined on how the tsunami might affect Banda Aceh's land use in the future. The dominant builtup area grew in the southern part of the city or the safe zone area. Not only in the safe zone area, but regional growth also occurred in areas affected by the tsunami. The total built-up area affected by the tsunami was 2/3 of the total built-up area of the city. This reveals the high vulnerability of the community to tsunamis in the future. Other types of land cover that were also affected, such as bare land, ponds, and the affected water area, will also increase in future tsunami scenarios. Bare land will experience a negative growth trend because most of it will be converted into built-up areas. Furthermore, the SLR will exacerbate the impact of the tsunami on the city. Thus, the built area, previously a safe zone, will become an affected area in 2065.
The results of this analysis provide important information for policymakers, especially urban spatial planning, which will determine the city's resilience to tsunamis. For example, the tsunami safe zone boundary map used for city spatial reference needs to be updated. The results of this study are also useful as initial information in assessing tsunami risk from the aspect of economic growth in the future. Increasing the built-up area also means increasing population. Therefore, it will also contribute to a larger human exposure. This result gives insight into the needs for sustainable and longterm tsunami risk reduction for Banda Aceh.