Numerical simulations of the 2004 Indian Ocean tsunami deposits' thicknesses and emplacements

After more than a decade of recurring tsunamis, identification of tsunami deposits, a part of hazard characterization, still remains a challenging task that is not fully understood. The lack of sufficient monitoring equipment and rare tsunami frequency are among the primary obstacles that limit our fundamental understanding of sediment transport mechanisms during a tsunami. The use of numerical simulations to study tsunami-induced sediment transport was rare in Indonesia until the 2004 Indian Ocean tsunami. This study aims to couple two hydrodynamic numerical models in order to reproduce tsunami-induced sediment deposits, i.e., their locations and thicknesses. Numerical simulations were performed using the Cornell Multi-grid Coupled Tsunami (COMCOT) model and Delft3D. This study reconstructed tsunami wave propagation from its source using COMCOT, which was later combined with Delft3D to map the location of the tsunami deposits and calculate their thicknesses. Two-dimensional horizontal (2-DH) models were used as part of both simulation packages. Four sediment transport formulae were used in the simulations, namely van Rijn 1993, Engelund–Hansen 1967, Meyer-Peter–Mueller (MPM) 1948, and Soulsby 1997. Lhoong, in the Aceh Besar District, located approximately 60 km southwest of Banda Aceh, was selected as the study area. Field data collected in 2015 and 2016 validated the forward modeling techniques adopted in this study. However, agreements between numerical simulations and field observations were more robust using data collected in 2005, i.e., just months after the tsunami (Jaffe et al., 2006). We conducted pit (trench) tests at select locations to obtain tsunami deposit thickness and grain size distributions. The resulting numerical simulations are useful when estimating the locations and the thicknesses of the tsunami deposits. The agreement between the field data and the numerical simulations is reasonable despite a trend that overestimates the field observations.

1 Introduction

In recent decades, sediment transport dynamics due to extremely long waves, such as tsunamis, has increasingly caught the interest of researchers, especially from countries in the direct path of these phenomena. Studies of tsunami impact on coastal morphology, involving a comparison with actual field data, are rare for previous tsunami events. This limits our ability to understand the physical processes occurring during erosion and sediment deposition. The 2004 Indian Ocean tsunami, triggered by a magnitude M_w=9.1 earthquake, dramatically changed coastal areas around Aceh of Indonesia (Borrero et al., 2006). Erosion processes and its effects were quite profound and severe along most of Aceh's west coast (Syamsidik et al., 2015a). Severe erosion associated with the 2004 tsunami actually created a new island (Al'ala et al., 2015). Tsunami waves transported large volumes of sediment of several hundred meters at certain locations, resulting in a large-scale coastline recession. Large shear stresses generated by the long waves overpowered the critical shear stress of bedload materials throughout the coastal zone. However, full investigation into the physical processes of longwave shear stresses during the 2004 tsunami did not occur. Owing to the lack of sufficient coastal monitoring equipment in tsunami-prone areas and the rare occurrence of tsunamis, insights into sediment transport mechanisms during tsunami events remain poorly understood. Overland flow velocities extracted from spontaneous eyewitness videos recorded during the 2004 Indian Ocean tsunami are limited to a few inland locations in Banda Aceh (Fritz et al., 2006). The coupling of tsunami hydrodynamics and sediment transport warrants future studies in order to improve our fundamental understanding of tsunami processes. Before the 2004 tsunami, paleo-tsunami studies were quite rare in this region (Andrade et al., 2014).

Two major factors, i.e., episodic land subsidence due to tectonic movement and massive sediment transport during tsunami propagation, play various roles in modifying coastal morphology (Moore et al., 2006). Several studies conducted in this region include the investigation of tsunami deposits at Khao Lak in Thailand (Dawson and Stewart, 2007; Jankaew et al., 2008) and Meulaboh in the Aceh Province (Monecke et al., 2008). Another study reported the presence of tsunami sediment deposits in a coastal cave, which is less than 5 km from our study area in Lhoong of Aceh Besar (Rubin et al., 2017). Both studies revealed that information from tsunami deposits is capable of improving our understanding of the physical processes involved during the transport of sediment. Paris et al. (2007) performed their study around Lhok Nga of Aceh Besar, which is about 20 km to the east of our study areas. Their study exhibited the influence of local topography on the sediment thicknesses found in the area. The thickest deposits were found in a low topography situation, and steep slopes gave varied results in spatial distribution of the tsunami deposit.

Previous studies on tsunami deposits have inspired tsunami mitigation efforts in the USA and in Indonesia (Dunbar et al., 2008; Rubin et al., 2017). Data on tsunami deposits would provide more scientific evidence for past tsunamis for a significant period of record, such as in the case of the coastal cave in Aceh that preserved a very long record of past tsunamis in the region (Rubin et al., 2017). A set of tsunami deposit data could help us to estimate the recurrence period of the tsunami from the subduction zone (Fujiwara, 2007; Minoura et al., 2001). Finding more tsunami deposits would strengthen the validation process and help find some ways to estimate the sources of tsunamis through inverse mechanisms (Buckley et al., 2012; Moore et al., 2011).

The recurrence rate of tsunamis throughout this region, information critical for hazard assessment, is inferable from events recorded in the tsunami deposits. Sediment transport due to tsunami waves has certain characteristic properties and hydrodynamic regimes during both the tsunami runup and backwash. During transport, material that is mobile depends on grain size and density, which vary greatly because of the erratic nature of tsunami sediment transport. Very few studies investigated the energy involved during sediment transport produced by the 2004 tsunami in Aceh. Furthermore, only a few of these published studies confirmed the thickness of the sediment deposits after the tsunami via geological surveys. This is the first study in northern Sumatra which applies numerical simulations in order to estimate the locations of tsunami deposits prior to the field surveys. Earlier investigations concentrated on forward tsunami modeling, such as in Kuala Meurisi (Apotsos et al., 2011a), Lhok-Nga (Gusman et al., 2012; Li et al., 2012), and Ulee Lheue Bay (Syamsidik et al., 2017). This research aims to couple two hydrodynamic numerical models in order to reproduce the spatiality of the 2004 tsunami sediment deposits, i.e., their locations and thicknesses. We analyze the possible sources of large shear stresses produced by the tsunami which transported sediment grains of various sizes. Several studies observed that backwash processes, associated with tsunami waves, deposited large volumes of sediments offshore (Jiang et al., 2015). Sediment properties, such as grain size, correlate poorly with the thicknesses of tsunami sediment deposits (Cheng and Weiss, 2013). Another poorly characterized physical process related to the tsunami is the location and extent of inland sediment deposition, which differs from the extent of tsunami inundation (Goto et al., 2012).

We combine both numerical simulations and field surveys to investigate sediment transport induced by the tsunami at Lhoong of Aceh Besar. We also conducted topographic surveys and sediment deposit pit tests in the vicinity. Studies rarely focus on changes along the shoreline due to the impacts of tsunami waves and subsequent erosion and deposition. Morphologic changes to the shoreline and sediment thickness are attributable to the amount of energy involved in tsunami propagation and inundation. This study estimates the energy required to transport sediments via a tsunami wave. The trajectory of the 2004 Indian Ocean tsunami was numerically simulated using the Cornell Multi-grid Coupled Tsunami (COMCOT) model, from its rupture area to the shoreline, and we applied Delft3D-FLOW to simulate sediment transport processes by inland flooding in the study area. The combination of these two models provides an ideal re-creation of the 2004 tsunami scenario, since COMCOT does not include sediment transport, and the original version of the Delft3D model does not have the option of multi-fault scenarios to generate tsunami waves (Syamsidik and Istiyanto, 2013). The numerical model is described briefly in Sect. 3, and the model results are presented in Sect. 4. Field observations are also described in Sect. 4. The observations and model results are discussed in Sect. 5, and conclusions are presented in Sect. 6.

2 Study area

To carry out the objectives of this study, we selected an area where human intervention has remained minimal during the decade following the 2004 tsunami in order to obtain well-preserved tsunami deposits. The study area is located at Lhoong in the Aceh Besar District, Indonesia (Fig. 1). The 2004 Indian Ocean tsunami destroyed a wide swath of coastal area, including the Aceh Besar District. The study area is situated on the west coast of Sumatra, approximately 40 km south of Banda Aceh, the city most devastated by the tsunami. Several prior studies have been conducted in the general study area (Jaffe et al., 2006; Rubin et al., 2017). Our study area is situated in a tropical area where the rainy season occurs about 4–5 months in a year, with a high precipitation rate. After more than 10 years, the 2004 Indian Ocean tsunami deposits in this area have experienced some natural processes despite the selection of the study area as the most sediment-preserved area in the Aceh Besar District. Szczucinski (2012) also argues that tsunami inundation less than 3 m would unlikely preserve the sediment deposit years after the tsunami event. Coastal mountains, representing extension of the Bukit Barisan Mountains, surround the study area. We selected four specific areas for detailed sediment transport investigations. The selected areas are shown in red boxes in Fig. 1, namely, Birek, Pasie Janeng, Jantang, and Saney. Among these locations, only Jantang has a wide coastal plain, and the other three sites have more rugged terrain. Birek and Pasie Janeng have V-shaped coastal areas characterized by steep hill slopes 100 to 300 m from their coastlines. Saney headland features a narrow plain sandwiched between coastline both to the north and south. To the north of Saney, there is an area which could be classified as plain coast, about 3 km long. However, this area has undergone heavy anthropogenic intervention after the tsunami as it was formerly a paddy field and settlement area. Therefore, we determine the four study sites as the most tsunami-deposit-conserved area.

Figure 1 Field survey sites at Lhoong, Aceh Besar, Sumatra, approximately 40 km south of Banda Aceh, Indonesia.

The 2004 tsunami heavily impacted the Lhoong region, leaving a number of areas conserved after the tsunami because of their relative isolation from major human activities. After more than a decade, it is difficult to find a suitable location to investigate tsunami deposits because many areas have been rebuilt, modifying the deposits during reconstruction processes. Birek, Pasie Janeng, and Saney are the sites that conserve the best sediment deposits because of the absence of anthropogenic influences that might disturb the tsunami deposits. Furthermore, researchers conducted a field survey of Jantang 3 months after the tsunami (Peters and Jaffe, 2010; Jaffe et al., 2006). Aerial photographs of the study area are shown in Fig. 2. Birek is surrounded by a range of hills covered by dense forests that regrew rapidly after the 2004 tsunami. The characteristics of Pasie Janeng are similar to Birek. Jantang (Fig. 2c and d), however, has a wider flat area that allowed the tsunami to propagate farther inland compared with Birek and Pasie Janeng. Since the 2004 tsunami, however, Jantang has recovered, new settlements have been established, and it is actively used for agricultural purposes.

Figure 2 Aerial drone photos of the study areas at Birek (a, b), Pasie Janeng (c, d), Jantang (e, f), and Saney (g, h) captured in October 2017.

3 Methods

3.1 Numerical simulations

During the simulation process, we combined two numerical models. COMCOT simulated the formation and propagation of tsunami waves. A series of laboratory experiments and prior tsunami events have validated the use of the COMCOT model (Liu et al., 1995). COMCOT was successfully implemented to investigate historical tsunami events, such as the 1992 Flores Island tsunami (Liu et al., 1995) and the 2004 Indian Ocean tsunami (Wang and Liu, 2006, 2007; Syamsidik et al., 2015b). Using the incident tsunami wave heights generated by COMCOT, Delft3D-FLOW models sediment transport from the ocean to the inundation zone, as well as the reverse flow. Both models coupled the hydrodynamic and morphodynamic models. The hydrodynamic model has been validated comparatively with numerical and laboratory data following several tsunami benchmark standards (Apotsos et al., 2011a). For the hydrodynamic model, Delft3D uses a finite difference scheme on a three-dimensional grid with nonlinear shallow water equations (NSWEs). The linear shallow water equations used in COMCOT in a spherical coordinate system are as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{R\cos \mathit{\phi }}}\left\{{\displaystyle \frac{\partial P}{\partial \mathit{\psi }}}+{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}\left(\cos \mathit{\phi }Q\right)\right\}=-{\displaystyle \frac{\partial h}{\partial t}},\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial P}{\partial t}}+{\displaystyle \frac{gh}{R\cos \mathit{\phi }}}{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{\psi }}}-fQ=\mathrm{0},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{gh}{R}}{\displaystyle \frac{\partial h}{\partial \mathit{\phi }}}+fP=\mathrm{0}.\end{array}$$

The equations governing nonlinear shallow water behavior in COMCOT for a spherical coordinate system can be expressed as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{R\cos \mathrm{\varnothing }}}\left\{{\displaystyle \frac{\partial P}{\partial \mathit{\psi }}}+{\displaystyle \frac{\partial }{\partial \mathrm{\varnothing }}}\left(\cos \mathrm{\varnothing }Q\right)\right\}=-{\displaystyle \frac{\partial h}{\partial t}},\end{array}$$

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle \frac{\partial P}{\partial t}}& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{R\cos \mathrm{\varnothing }}}{\displaystyle \frac{\partial }{\partial \mathit{\psi }}}\left\{{\displaystyle \frac{P^{\mathrm{2}}}{H}}\right\}+{\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial \mathrm{\varnothing }}}\left\{{\displaystyle \frac{PQ}{H}}\right\}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{gH}{R\cos \mathrm{\varnothing }}}{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{\psi }}}-fQ+F_x=\mathrm{0},\end{array}\end{array}$$

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\partial Q}{\partial t}}& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{R\cos \mathrm{\varnothing }}}{\displaystyle \frac{\partial }{\partial \mathit{\psi }}}\left\{{\displaystyle \frac{PQ}{H}}\right\}+{\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial \mathrm{\varnothing }}}\left\{{\displaystyle \frac{Q^{\mathrm{2}}}{H}}\right\}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{gH}{R}}{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathrm{\varnothing }}}+fP+F_y=\mathrm{0},\end{array}\end{array}$$

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }f=\mathrm{2}\mathrm{\Omega }\sin \mathit{\phi },\text{(8)}& {\displaystyle }& {\displaystyle }{F}_x={\displaystyle \frac{g{n}^{\mathrm{2}}}{H^{\mathrm{7}/\mathrm{3}}}}P{\left(P^{\mathrm{2}}+Q^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{F}_y={\displaystyle \frac{g{n}^{\mathrm{2}}}{H^{\mathrm{7}/\mathrm{3}}}}Q{\left(P^{\mathrm{2}}+Q^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\text{(10)}& {\displaystyle }& {\displaystyle }H=\mathit{\eta }+h.\end{array}$$

Meanwhile, for a Cartesian coordinate system, COMCOT applies the following equations for linear shallow water equations.

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+\left\{{\displaystyle \frac{\partial P}{\partial x}}+{\displaystyle \frac{\partial Q}{\partial y}}\right\}=-{\displaystyle \frac{\partial h}{\partial t}},\text{(12)}& {\displaystyle }& {\displaystyle \frac{\partial P}{\partial t}}+gh{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}-fQ=\mathrm{0},\text{(13)}& {\displaystyle }& {\displaystyle \frac{\partial Q}{\partial t}}+gh{\displaystyle \frac{\partial h}{\partial y}}+fP=\mathrm{0}.\end{array}$$

For linear shallow water equations, the following formulae were applied in COMCOT.

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+\left\{{\displaystyle \frac{\partial P}{\partial x}}+{\displaystyle \frac{\partial Q}{\partial y}}\right\}=-{\displaystyle \frac{\partial h}{\partial t}},\text{(15)}& {\displaystyle }& {\displaystyle \frac{\partial P}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left\{{\displaystyle \frac{P^{\mathrm{2}}}{H}}\right\}+{\displaystyle \frac{\partial }{\partial y}}\left\{{\displaystyle \frac{PQ}{H}}\right\}+gH{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}+F_x=\mathrm{0},\text{(16)}& {\displaystyle }& {\displaystyle \frac{\partial Q}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left\{{\displaystyle \frac{PQ}{H}}\right\}+{\displaystyle \frac{\partial }{\partial y}}\left\{{\displaystyle \frac{Q^{\mathrm{2}}}{H}}\right\}+gH{\displaystyle \frac{\partial \mathit{\eta }}{\partial y}}+F_y=\mathrm{0},\end{array}$$

where η is the water surface elevation, (P, Q) gives the volume flux in the X (west–east) direction and in the Y (north–south) direction, respectively, (φ, ψ) is the latitude and longitude, R is the Earth's radius, g is gravitational acceleration, and h is water depth. The $-\partial h/\partial t$ expression reflects the effect of transient seafloor motion, f represents the Coriolis force coefficient due to the Earth's rotation, Ω is the Earth's rotation rate, H is the total water depth, F_x and F_y represent bottom friction in the ψ and φ direction, respectively, and n is Manning's roughness coefficient. Since Delft3D-FLOW is unable to generate tsunamis formed from a multi-fault scenario, we therefore used COMCOT to propagate tsunami waves from layers 1 to 3. In the nested layer 4, we ran Delft3D-FLOW with boundary conditions based on water elevations produced in layer 3. We were able to vary Manning's roughness coefficient based on actual conditions, in reference to natural channels and floodplains (Arcement and Schneider, 1989). An open boundary condition for Delft3D-FLOW was placed about 14.5 km from shoreline. The hydrodynamic boundary condition was obtained from the COMCOT numerical simulation at layer 3. Here, the boundary condition was set to follow the water level. The water level boundary condition is a modified Riemann boundary condition to which Stelling added the time derivative of the water level and velocities (Stelling, 1984). This was meant to reduce the reflection process caused by the eigenfrequency of the simulation.

The four simulation domain layers applied to COMCOT simulation all used spherical coordinates. Figure 3 shows the layers of the simulations used in COMCOT. Layer 4 was later converted to the modeling domain for use in the Delft3D simulations. Coordinates, grid size, and types of shallow water equation (SWE) used in each individual layer are shown in Table 1. Delft3D-FLOW used the finest COMCOT layer as a domain, with coordinate transformation from spherical to Cartesian. We combined GEBCO and other nautical charts to produce improved bathymetry data. For the inland regions, we combined Shuttle Radar Topography Mission (SRTM) 30 with some cross-topography measurements to improve the topographical data.

Figure 3 Nested simulation layers incorporated into COMCOT and the boundary condition location for the Delft3D simulation (a solid blue line) in layer 4.

Table 1 Specification of layers applied in the COMCOT simulation.

To simulate the rupture area caused by the M_w=9.1 earthquake on 26 December 2004, we used a model proposed by Piatanesi and Lorito (2007) and Romano (2009). The rupture that generated the initial waves around the source area is shown in Fig. 3. We divided the rupture area into eight segments as described in Piatanesi and Lorito (2007). This model simulated a multi-fault scenario and was validated by Syamsidik et al. (2015b). Using NSWEs for layers 1 to 3 provides insignificant differences to the waves generated around the offshore area. On the other hand, applying nonlinear SWEs for layers 1 and 3 will reduce computational time without reducing the accuracy of the simulation results.

We initiated sediment transport using a grain diameter switch for each of the three sites. By producing shear stress, we calculated a uniform Manning's roughness coefficient of 0.02 for the seafloor and a coefficient of 0.04 for overland areas combined with a water level fluctuation boundary condition from COMCOT simulations. We used Delft3D to simulate sediment transport processes caused by tsunami waves. We modeled this at the innermost layer of the simulation domain (layer 4 in the COMCOT model). Delft3D uses NSWEs combined with finite difference methods. The NSWEs in the model apply the equations of conservation mass, momentum, and energy flux and have the capability to simulate rapidly varying flow conditions (Stelling and Duijmeijer, 2003). In its basic module, Delft3D computes sediment transport using bedload and suspended load characteristics as described in van Rijn (1993, 2007). Changes in topography caused by erosion and sedimentation processes are updated at each time step of the simulation (Lesser et al., 2004). A more detailed explanation of the sediment transport formulas used in Delft3D can be found in Apotsos et al. (2011a). One study proved that van Rijn (1993) formulas provided the best results compared to tsunami-induced sediment transport experiments in a wave flume (Li and Huang, 2013). Notwithstanding the reviews, we incorporated more sediment transport formulas in our numerical simulation, such as Engelund–Hansen 1967, Meyer-Peter–Mueller (MPM) 1948, and Soulsby 1997. Recently, researchers have used Delft3D to simulate and test extreme events that triggered morphological changes, such as in the case of the Indian Ocean tsunami (Syamsidik et al., 2017; Gelfenbaum et al., 2007; Apotsos et al., 2011a, b, c). Previously, Delft3D was used and validated based on conditions less extreme than a tsunami, such as wind–wave-driven sediment transport (Broekema et al., 2016; Lesser et al., 2004; van Rijn et al., 2011).

3.2 Field measurements

The results of the numerical simulations present estimations of the locations of the tsunami deposits, as well as their respective thicknesses. To assess the validity of our results, we conducted a series of surveys at three locations, i.e., Birek, Pasie Janeng, and Saney. We conducted these surveys in October 2015, January 2016, and April 2016 for Birek, Saney, and Pasie Janeng, respectively. Surveys were conducted during different periods of the year because of the intensive labor and time needed for each survey. At each location, we performed manually excavated pit tests.

To define the tsunami deposit, we followed the USGS procedure as described in Peters and Jaffe (2010). The same procedure was also applied during a tsunami deposit survey in Jantang conducted by Peters and Jaffe (2010). The identification of tsunami deposits relies on 10 criteria. Tsunami deposits are easily identifiable by their sharpness and erosional basal contact with original soil below or under the tsunami deposit. In this study, tsunami sediments were easily identified because of the characteristics of adjacent sediments, which have high organic contents from surface runoff processes. Tsunami sediments form thin layers, ranging between 1 and 30 cm thick in most cases (Srinivasalu et al., 2009). Beaches near the study area are composed of sandy-quartz materials. Patches of coral reefs exist around Birek, Saney, and Pasie Janeng. Therefore, it was difficult to distinguish between tsunami deposits and sediments produced by littoral transport, also highlighted in Jaffe et al. (2003). The presence of sea shells was identified through a microscopic observation of the sediment material. This could distinguish the tsunami-induced deposit from original topsoil material or other surface runoff process. Other criteria of the tsunami deposit were advised by Jaffe et al. (2003) and Peters and Jaffe (2010) who put the methods of the tsunami-induced sediment transport investigation into one practical guideline. We followed the steps of the guideline exactly and clarify it using some microscopic observations. We spent a significant period (from 0.5 to 2.0 h) for each sample to carefully identify the tsunami deposit and distinguish it from other sources of sediment. In total, we collected 14 samples of the tsunami deposits and 22 locations of sampling performed by Jaffe et al. (2006). Rip-up clasts were also found in the layers of the tsunami deposit layers, indicative of the energy that transported the material, such as in the case of a tsunami. Therefore, the existence of rip-up clasts further confirms the presence of tsunami deposits (Morton et al., 2007). Costa et al. (2012, 2015) proposed that the shape of the zircons in the sediment deposit could be used to interpret the number of waves and tsunami runup or backwash processes. Costa et al. (2015) stated that euhedral zircon could be associated with backwash processes. Meanwhile, rounded zircons could be attributed to deposition that occurred during tsunami runup processes. Detailed observation is absent in our present study.

4 Tsunami runup and tsunami deposits

We analyzed the impacts that the 2004 tsunami had on coastal morphology using field surveys and satellite images. Tsunami-inundated areas were observed using satellite images to measure changes in coastal morphology in the vicinity of Lhoong, Aceh Besar. We identified tsunami deposits during a field survey in October 2015. Surveys started with coastal topography measurements and then continued by location of potential deposits based on topography and surrounding conditions.

4.1 Tsunami inundation area

We analyzed the inundation area caused by the tsunami using satellite images captured in February 2005 using DigitalGlobe. After the 2004 tsunami propagated around the rupture source, waves moved inland, reaching the foothills approximately 200 m from the coastline at Birek and approximately 1 km inland at Jantang. The extent of inundation is shown in Fig. 4 by comparing the boundaries from numerical simulations and satellite images. The inundation limits generated by numerical simulation and those shown in the digitized images are generally consistent. On the basis of the inundation zones, we selected three locations to perform pit tests to sample the tsunami sediment deposits. For Jantang, we refer to descriptions in Jaffe et al. (2006).

Figure 4 The 2004 tsunami inundation limit at Lhoong based on the numerical simulation (dashed black lines) and satellite image digitization (dashed red lines).

4.2 Coastal topography

Surrounded by hills, the topography conditions in the Birek coastal area are classified as V-shaped, indicating that the elevation of the center region is lower than its surroundings. At Birek, a small creek channels surface runoff to the sea. Figure 5 shows the topography of Birek and cross profiles of the two transects, i.e., Cross 1 and Cross 2. The locations of the pit tests for the tsunami deposit sample collected are marked as B01–B03 in Fig. 5.

Figure 5 The topographic profiles at Birek and the locations of the pit tests (B01–B03), Pasie Janeng (P01–P04), and Saney (S01–S05).

Pasie Janeng has a limited plain area, stretching approximately 100 m from the shoreline. Steep slopes characterize the inundation zone further inland at Pasie Janeng (Fig. 5b). Given this topography, we estimate that surface runoff from rainfall, flowing from the hillside to the beach, significantly affects the current thickness of tsunami deposits. We assume that hydrologic runoff processes are the main mechanism that erodes the tsunami deposits. At Saney headland (Fig. 5c), topography is flat, but the study area is sandwiched by the coastline. We conducted pit tests to investigate tsunami deposits at the Point S01 to Point S05 locations.

4.3 Sediment characteristics

Tsunami deposits were identified by the sharp, clear presence of a sandy layer between new and old topsoil. The presence of seashells also strengthened the indication of tsunami deposits (Jaffe and Peters, 2010). Sediment deposition from the tsunami waves almost reaches the foothills. The energy produced by the earthquake was capable of generating tsunami waves that deposited sediments up to 6.3 m in elevation above sea level. The thicknesses of the tsunami deposits correlate with distance from the shoreline, becoming thinner the farther inland. At sampling location B01, we did not find tsunami deposits. At location B02, tsunami deposit thickness was approximately 29 cm. The location of Point B02 is approximately 320 m from the coastline. At sampling location B03, we identified a thinner tsunami deposit, approximately 7.5 cm thick. This location is approximately 400 m from the coastline (Fig. 5). Pit dimensions at these two locations were approximately 1.8 m in length, 0.9 m in width, and 1.50 m in depth.

The profiles of the pit tests at Birek, Pasie Janeng, and Saney are shown in Fig. 6. All of the profiles have similar trends in thicknesses of tsunami deposits.

Figure 6 Profiles of tsunami deposit samples B02 and B03 at Birek, PJ01 and PJ03 at Pasie Janeng, and S01 and S02 at Saney.

The majority of the sediment deposits at several locations in Lhoong contain sand, i.e., >70 % sand (Birek = 82.15 %, Pasie Janeng = 78.85 %, and Saney = 71.88 %). Waves carried the majority of the sand from the ocean bed, which mixed with deposits from each locality. Each sampling location contains less than 25 % coarse grains (Birek = 6.3 %, Pasie Janeng = 17.29 %, and Saney = 20.46 %), which is produced by local erosion. The tsunami flow depth drove the sediment transport process, which resulted in the characteristic sandy tsunami deposits. Additionally, we found rip-up clasts in the tsunami deposits at Birek and Pasie Janeng.

4.4 Model hindcast

To simulate both hydrodynamic and morphodynamic processes during the 2004 Indian Ocean tsunami, COMCOT and Delft3D were run simultaneously. Transient waves produced by COMCOT successfully drove the sediment transport simulation performed by Delft3D. Tsunami wave heights, shear stresses, and the sediment transport process produced by Delft3D-FLOW are useful when trying to explain the sediment transport mechanism during the runup of tsunami waves. Here, we present the results of the sediment transport process simulations in cumulative sedimentation and erosion thicknesses. Extreme changes in coastal morphology from plains to steep slopes significantly contributed to the wide distribution of tsunami deposits throughout the area affected by the tsunami.

4.4.1 Tsunami wave height and sedimentation

Tsunami occurrence induced tsunami waves, inundating areas far inland. These results came first from COMCOT, which produced the tsunami wave height that was used as a parameter for Delft3D simulations. The modeled tsunami consisted of two major waves: (1) the first wave, which peaked at a height of 5.23 m, and (2) the second wave, which peaked at a height of approximately 5.61 m (Fig. 7).

Figure 7 Delft3D-FLOW boundary flow input from COMCOT results.

On the basis of the COMCOT results, Delft3D-FLOW, at a nearshore observation point, produced two major tsunami waves: (1) the first wave peaked at 12.2 m, and (2) the second wave followed 14 min later peaking at 8.9 m, 3.3 m lower than the first wave. The fluctuation in runup and backwash generated significant sediment transport in coastal waters and overland flow.

Sedimentation processes occurred predominantly during backwash. Specifically, such a case only occurs on relatively mild slopes. On the other hand, steep slopes remove sediments during the backwash. Reduced energy after the second and subsequent waves carried more sediment deposits but had significantly less energy during the backwash process, therefore leaving more sediments behind (08:45 LT – local time). As shown in Fig. 8, backwash produced a sediment deposit that was 0.38 m thick during the second wave.

Figure 8 Relation between wave height and cumulative sediment on the observation points (shown in red dots).

4.4.2 Tsunami velocity fields

The sediment transport process is highly influenced by the shear process induced by tsunami currents along the bottom of the bed. The velocity fields during the tsunami runup process help in explaining the energy involved during sediment transport processes. Figure 9 shows the velocity fields at the Birek and Pasie Janeng locations.

Figure 9 Tsunami wave velocity fields at Birek during wave advance (a) and during wave retreat (b) and at Pasie Janeng during wave advance (c) and during wave retreat (d).

As shown in Fig. 9a and b, during the wave's advance, the velocity at Birek location is lower than that when the wave is retreating. We also found the same condition at Pasie Janeng (Fig. 9c and d). The maximum value of shear stress was 83 N m⁻², recorded 24 min after the first wave arrived, eroding 0.11 m of the bed layer. High shear stress values significantly erode land, still depositing sediments 42 min after the first waves arrive.

4.4.3 Tsunami deposit thickness

On the basis of the simulations using both Delft3D and field measurements, we compare the differences in tsunami deposit thicknesses between the four sample locations. The runup process mostly produces erosion, whereas most backwash does the opposite. The accumulation sedimentation map gave relatively reliable information on tsunami deposit locations, similar to the simulation results found for Lhoong. Figure 10 shows maps of spatial distribution of cumulative sedimentation and erosion after the tsunami waves based on the numerical simulation results. Negative values are for erosion, and positive values are for sedimentation thickness in meters. The following points compare the results of the numerical simulations and the field data from the four sample locations.

Figure 10 Tsunami deposit thickness from numerical models using Delft3D-FLOW at Birek, Pasie Janeng, Jantang, and Saney.

a. Birek

Hydrodynamic process analysis revealed that there is some wave amplification, along with a decrease in the current velocity. Elevated topography halted the tsunami runup at the study area. The elevated topography deflected tsunami waves and reduced their velocity. The reduced tsunami flow velocity also signifies that the tsunami waves are losing energy. This explains the presence of a thinner tsunami deposit, which we found around the foothills compared with areas farther from the hills. The deflected waves created a larger wave around the hill, but at the same time, there was a significant decrease in the shear stresses. The lower the value of the shear stress that is generated by the waves, the more likely it is for sediment deposition to take place. Figure 11a compares the results of the numerical simulation and field data at Birek.

Figure 11 Comparison of tsunami deposit thickness results from the numerical models using four different sediment transport formulae and the field survey at Birek, Pasi Janeng, Jantang, and Saney.

The coastline is marked by the origin along the x axis. Data from the three pit tests show that the numerical results overestimate the field data. It is important to note that this location has a V-shaped topography, surrounded by hills. The deposit thickness decreases with increasing distance from the coastline.

b. Pasie Janeng

Compared with Birek, we found slightly different results at Pasie Janeng (see Fig. 11b). Areas with mild slopes are shorter than those at Birek, limiting the area capable of trapping sediments. The steep topography of the hills surrounding Pasie Janeng funnels rainfall runoff into a small channel that drains directly to the sea. The path that runoff took eroded a small amount of the tsunami deposit that could have been deposited there for almost 12 years. The morphology of the site resembles an embayment, which potentially receives new sediment from coastal processes, making it difficult to find tsunami deposits nearshore.

c. Jantang

Jaffe et al. (2006) examined the thickness of tsunami deposits at Jantang during a survey conducted in 2005. The wide area, with a mild slope, was advantageous for capturing sediments during the tsunami event. The relatively short time period between the tsunami and the survey was beneficial given that there was little activity in the area (i.e., not until 2009) just after the 2004 tsunami, concerning rehabilitation and reconstruction. Transect data based on field surveys were compared to the numerical simulations. Comparisons between the field data and the numerical simulation results are shown in Fig. 11c.

The results from the numerical simulations and the field survey are in good agreement. Unlike the data from the other three locations, the results of the numerical simulations at Jantang underestimate the field data. The absence of variations in the roughness coefficient slightly affected sediment transport results. The lack of additional land cover parameters from corals or roads also affected several of the model-calculated deposition and erosion zones. Besides these limitations, the simulation results are in good agreement with the field data.

d. Saney

Saney is protected by a small hill facing the ocean situated at the edge of the land, which made it more vulnerable to tsunami wave energy. A settlement complex was completely swept away during the 2004 tsunami. The tsunami wave that came from the west struck the edge of the land and was amplified after encountering the mild slope area. The narrow headland of the site provided a small place for the tsunami deposit to settle down. Years after the tsunami, other factors may have eroded the deposit, leaving thinner tsunami deposit. The results of the numerical simulations overestimated values compared with the field data (see Fig. 11d). Regardless of the process, we could still find the tsunami deposit 13 years after the event with the support of numerical modeling.

5 Discussion

Estimations of the deposit locations and the thickness of tsunami deposits are key factors for site selection and target fieldwork efforts. Even 13 years after the Indian Ocean tsunami, numerical simulations still provide useful estimates in order to locate tsunami deposits. Finding suitable tsunami deposit locations may become increasingly difficult as massive reconstruction efforts significantly alter the tsunami deposits in certain areas. The tsunami deposits may ultimately disappear or be difficult to identify. It is important to discuss land use practices with the local community in order to target well-preserved field sites. Using a linear regression method, we found that the van Rijn (1993) formulas gave the best approximation of the field data. Similar results were found at the other three study sites, with a smaller correlation (r²). Our findings confirmed suggestions from laboratory experiments done by Li and Huang (2013).

We encountered the following limitations during the course of this study. First, field measurements could not provide the exact number of tsunami waves. On the basis of our models, two main waves occurred during the tsunami in this region. This is corroborated by several eyewitness accounts from the study area given the absence of instrumental wave recordings. The small number of waves determined the size of the material contained in the deposits, such as the rip-up clasts in the tsunami deposit layer. We found rip-up clasts at both Pasie Janeng and Birek. This confirms that there were only a few waves, as these differ from storm waves (Morton et al., 2007). Secondly, there are uncertainties about the model itself. For sediment transport, we used the two-dimensional horizontal (2-DH) Delft3D that did not account for turbulence kinetic energy, which may simplify the results. Furthermore, our simulations used depth-averaged velocity produced by the tsunami waves. This may, in fact, limit the maximum suspended sediment concentration because of the exclusion of suspended-sediment-induced density stratification (Jaffe et al., 2016). Other studies have attempted to run Delft3D with the 3-D model that allows interlayer settling and erosion processes during tsunami wave propagation (Apotsos et al., 2011b; Gelfenbaum et al., 2007). Our study area has a complex topography and morphology. All of the sites, except for Jantang, have headlands and patches of coral reefs. The complexity of the topography could produce a significant difference between numerical model results and field data, as observed in the case of the Sanriku coast in Japan (Goto et al., 2017). The roughness coefficient, which was set to a value similar to all grids of the sea area, could contribute to the different thicknesses found between the numerical simulations and the field data.

Regardless of the limitations and uncertainties, this study has successfully demonstrated the effectiveness of simulations in estimating the locations of tsunami deposits. Although natural processes such as rainfall, runoff, and aeolian sediment transport have taken place in the study area for more than a decade, location estimation produced by the numerical simulations is correct. The methods presented here are applicable to areas without any human activity or other extreme events, such as the effects of storms, that could not be considered in the simulations. Discontinuous and patchy tsunami deposit records may be the result of anthropogenic disturbances (Chague et al., 2018).

Further advances in the field of tsunami-induced sediment transport should incorporate collaborations between multidisciplinary researchers, as this may increase the comprehensiveness of the study. The limitations and uncertainties highlighted in this study remain major challenges in paleo-tsunami studies (Sugawara et al., 2014). A solid understanding of paleo-tsunamis, by exhibiting tsunami deposits to coastal communities, serves to allow us to better communicate the risk of tsunamis to them. We also acknowledge that the use of high-resolution topography and bathymetry data could significantly contribute to the reliability of the numerical simulation results.

6 Conclusions

This study has successfully demonstrated the coupling method of COMCOT and Delft3D-FLOW to provide a better understanding of tsunami-wave-induced sediment transport and to find the locations of tsunami deposits in a specific area more than one decade after the tsunami. Large shear stresses generated by tsunami waves transported sediment farther inland. Field survey topography analysis also contributed to increased accuracy in defining the deposit location inside the tsunami inundation zone. Combining both methods contributed to results that are in good agreement between numerical simulations and field data. Tsunami deposit surveys were conducted at three locations between 2015 and 2016, showing that actual tsunami deposit thicknesses are thinner than what the numerical results show. Surface runoff, soil consolidation, and aeolian sediment transport could contribute to these differences. On the other hand, tsunami deposit data, obtained just after the tsunami, provided a good fit to the numerical results. This study demonstrates the importance of rapid data collection in the immediate aftermath of a tsunami. Delft3D results displayed the location of sediment deposition, and field surveys confirmed the presence of tsunami deposits physically. Ultimately, the characterization of modern tsunami deposits will allow us to interpret earlier extreme events in sedimentary records.