Estimation of tsunami debris on seafloors towards future disaster preparedness_ Unveiling spatial varying effects of combined land use and oceanographic factors

A large amount of tsunami debris from the Great East Japan Earthquake in 2011 was sunk on the seafloor and threatened the marine ecosystem and local communities' economy, especially in fisheries. However, few studies estimated spatial accumulations of tsunami benthic debris, comparing to their flows on the ocean surface. Here, a spatially varying coefficient model was used to estimate tsunami debris accumulation considering the spatial structure of the data off the Tohoku region. Our model revealed the number of vessels nearest the coast at the tsunami event had the highest positive impact, whereas the distance from the coast and kinetic energy influenced negatively. However, the effect of the proximity to the coast wasn't detected in the Sendai bay, indicating spatial dependency of these effects. Our model estimation provides the fundamental information of tsunami debris accumulation on the seafloor, supporting early reconstruction and risk reduction in marine ecosystems and local communities.


Introduction
The tsunami generated by the Great East Japan Earthquake in 2011 washed 5 million tons of debris into the North Pacific (Ministry of the Environment, 2012). The Japanese government estimated 3.2 million tons of debris sank in the coastal area (Ministry of the Environment, 2012). Tsunami debris accumulated on the seafloor would have an impact on marine organisms (Brown and Macfadyen, 2007;Chiappone et al., 2005) and damage fishing gears including trawl nets (Katsanevakis et al., 2007;Nash, 1992), although some of the debris may act as an artificial reef for some species (Girard et al., 2004;Gooday, 2002). The seabed debris removal should be efficiently advanced immediately after the disaster not only for the conservation of coastal and offshore marine life, but also to reduce the possible impact on human activities such as fishery in coastal and offshore areas (Ramachandran et al., 2005).
An approach using statistical or mathematical models can estimate the extent of debris deposition at a particular location by corresponding environmental conditions where debris is likely to be deposited. However, while model estimation of drifting debris has been relatively prolific (Kawamura et al., 2014;Lebreton and Borrero, 2013;Maximenko and Hafner, 2010;Prasetya et al., 2012), few studies have modeled the spatial accumulation of tsunami debris on the seafloor, which means that detecting debris on the seafloor is often much more difficult than for drifting debris, resulting in a paucity of researches (Barnes et al., 2009;Goto and Shibata, 2015;Watters et al., 2010). Although observational studies have reported that debris deposited on the seafloor tends to be influenced by geographic features such as distance from coasts and bathymetry (Bauer et al., 2008;Debrot et al., 2014;Galgani et al., 2000), hydrodynamics (Galgani et al., 2000), and the degree of urbanization and anthropogenic activities in surrounding areas (Loulad et al., 2019;Pham et al., 2014), the influences of these factors are not always consistent across the studies. This suggests there are combined effects of inland land use, ocean currents, and seabed topography conditions at the subject site, and their spatial dependency. Exploratory cleaning operations in such an environment can be more inefficient. In fact, the debris removal work that was carried out after the Great East Japan Earthquake is still being carried out eight years after the disaster (JUDGIT, https://judgit.net/projects/2513. last accessed 8 May 2020). Previous disaster cases have provided simulations of debris generation on land and are beginning to provide knowledge for planned cleanup work (Srinivas and Nakagawa, 2008;Wood and Good, 2004), but for debris deposited on the seabed, reports on the distribution of recovered debris are limited (Goto and Shibata, 2015). If the model can be used to estimate the locations where debris is likely to accumulate based on environmental conditions, it will be possible to more efficiently select the priority locations for cleaning operations (Bauer et al., 2008).
This study aimed to develop a model to estimate the location of tsunami debris accumulation on the seafloor using a suite of environmental variables. Our study area was in the northeastern part of Japan, which was tremendously damaged by the Great East Japan Earthquake in 2011. After the disaster, Japanese and local governments conducted surveys investigating tsunami debris on the seafloor, allowing us to overcome the lack of tsunami debris data in this study. Here, we investigate the probability of tsunami debris accumulation on the seafloor using a spatially varying coefficient model that considers explanatory factors of both oceanographic (bathymetry, ocean velocity, tsunami impact) and land-based resources (land use coverage and number of vessels in the inundation area) in addition to the coast proximity.

Study area and tsunami debris survey
The survey area is off the coast of Miyagi and Iwate prefectures in Japan, which is part of the North Pacific, i.e., ca. 100 km from the epicenter of the earthquake (Fig. 1). The area severely suffered from the Great East Japan Earthquake in 2011, and a large amount of debris was washed to the near shore by the tsunami runup (NOAA Marine Debris Program, 2015).
We used debris data obtained from a survey of sediment and debris off the coast of Miyagi Prefecture in the process of cleanup operations. The survey was performed along the coast area ( Fig. 1; from 39°5′N; 141°7′E to 37°8′N; 140°9′E) using a SONICE2024 multibeam sonar system (frequency, 200-400 kHz; swath width, 1°× 1.6°; beams, 256) from October 2012 to March 2013 (Department of Agriculture, Forestry, and Fisheries, Miyagi Prefecture provided). The horizontal measurement accuracy was < 0.6 m, and debris larger than 2 m × 2 m was recorded based on seafloor topographic information acquired at a resolution of 1 m × 1 m. The survey recorded the volume (m 3 ), type (ship, fishing gear, container, house, timber, car, artificial reef, accumulation, unknown items), location (latitude and longitude), and depth (m) of each debris. In the following analysis, artificial reefs were omitted, and the total volume of all types summed up per 1 km grid square was used. It was due to classification uncertainty of debris types in the survey data. Additionally, in some cases, multiple pieces of debris were intertwined and drifted or settled, and there were limits to what could be determined by estimating based on the type of information, so the analysis was conducted without distinguishing between types.

Environmental data
We modeled the probability of tsunami debris accumulations on the seafloor using variables related to oceanographic (bathymetry, ocean velocity, tsunami impact) and land-based resource (land use coverage and number of vessels in the inundation area) factors. The bathymetric gradient (depth) at 30 m resolution was derived from the M7000 Digital Bathymetry Chart provided by Japan Hydrographic Association to calculate the range of depth (m) at 1 km grid (Fig. 2a). Ocean velocity components (u, v) were obtained at 1.7 km resolution from numerical model products (Tanaka et al., 2018) to calculate cumulative values using monthly averages of surface kinetic energy (cm 2 s −2 ) from March 2011 to 2013. Kinetic energy explains the flow intensity on that site, with a high value indicating strong intensity of ocean flow in the area (Fig. 2c). The variable regarding the bottom kinetic energy was not used for the following analysis because of small values in our study area (the average cumulative values at the bottom kinetic energy from March 2011 to 2013 were 0.001 cm 2 /s −2 ± 0.007 cm 2 /s −2 ). The volume of debris (m 3 ) generated by tsunami was considered to be affected by tsunami height (m) and land-use type in the inundation area ( Fig. 2c, d, f, g). The values in the tsunami height along the coast were derived from tsunami survey results (Mori et al., 2011), and the maximum value at the nearest coast from each grid was used in the analysis. The coverage of each land use type in the nearest inundation area was calculated using the land use map provided by the Geospatial Information Authority in Japan. The following types were included: the urbanized area (m 2 ), riverine area (m 2 ), forest area (m 2 ) and crop land (m 2 ) ( Fig. 2c, f, g; Only the variables included as explanatory variables are shown in the figure. Details of the selection process are given in Section 2.3). For each land use, the total area included in a buffer zone 10 km from the nearest shore, which covers the majority of the inundated area, was calculated. The number of fishing vessels was derived from the 2008 census report of fisheries by the Ministry of Agriculture, Forestry, and Fisheries, with the latest version before the disaster to calculate the average numbers in the same buffer for landuse data (Fig. 2e). The distance from the coastline (km), the range of depth, each area of land use type, and the number of vessels were calculated for each 1 km grid using the ArcGIS version 10.3 (ESRI Inc.).

Model description
A spatially varying coefficient model (Murakami et al., 2017) was applied to estimate the probability of the tsunami debris accumulation on the seafloor to account for spatial autocorrelation of the debris accumulation and environmental variables. This method is appropriate to apply for this study region, which is likely to be affected by different factors at different locations. It was originally developed to detect the difference in spatial coefficient structures using eigenvectors of the proximity matrix (Dray et al., 2006). Each of them is computed by correlating the spatial distributions of a response and the environmental variables with each other. By using the eigenvectors of the proximity matrix to model spatial correlations, it is possible to represent covariance structures that cannot be explained by distance alone, such as differences in spatial correlation structures by region (e.g., a debris is less influenced by the current in the bay). This assumption enabled us to represent the spatial autocorrelation that may affect the accuracy and uncertainty of the parameter estimation in models (Griffith and Peres-Neto, 2006). In this study, Moran I statistics (Moran, 1950), popular spatial correlation statistics, was used to calculate the eigenvector of the proximity matrix (Murakami et al., 2017(Murakami et al., , 2019. We used the total volume of each debris type log-transformed for the response variable (0.5 was used as added constant included in our data; see Yamamura, 1999). To reduce any potential multicollinearity, variance inflation factors (VIF) (Heiberger, 2019) were calculated using these variables. The final environmental variables included in the model were the distance from the coast, surface kinetic energy, range of depth, tsunami height, coverage of urbanized area, riverine area, and number of vessels as a spatially varying coefficient variable. The distance from the coast was assumed to follow an inverse square law. This inverse square law is one of the universal physical laws expressing an effect is proportional to the inverse square of the distance from its source. We assumed various effects are expected to work on the discharge of debris to the sea from the land, but these effects would be attenuated with distance from the shore.
Each parameter of the environmental factors was estimated by assuming that it would vary spatially. For example, tsunami height and land use factors can have greater impact nearer the coast, whereas the ocean velocity may have greater effect further from the coast. The model is specified in the following form where VOL log is the estimated volume of tsunami debris. The formula (1) consists of three terms explaining global variations, local variations, and observation errors, where β k is n × 1 vectors of this model on kth explanatory variables x k where n is the number of grid squares. β k is composed of β k, 0 1 as a spatially constant component (1 is a vector of ones, that is there is no heterogeneity of eigenvector) and Eγ k as a spatially varying component. Eγ k is expressed using Moran's eigenvectors, E, and their random coefficients of a given explanatory variable, γ k . Covariance matrices (n × k) are expressed using the form of a tensor product, ⊗ to explain the spatial filtering effects of each explanatory variable on the tsunami debris. Error terms, ε followed the gaussian process, N(0, σ 2 Ι), which underlined the hypothesis that the Gaussian process would compensate for the observed patterns of tsunami debris caused by a spatial autocorrelation using sonar investigation (Fig. 1). The number of iterations was 30 times (default setting) to gain robust results, which means that the model is not just a result of random probability. For comparison, a generalized linear model was constructed with the same dataset where all variables were set as a constant variable and a null model without considering spatial autocorrelation of the data following a Gaussian distribution. Model performance was evaluated using the root mean square error (RMSE) of all possible models, which is a popular indicator to show deviation between actual and predicted values (smaller value indicates better predictive power of a model). Model analyses were performed using the R package "spmoran" (Murakami et al., 2017) and "gdistance" (van Etten, 2017). We projected the model constructed to the fishery area outside Tohoku ( Fig. 1; from 37°9′N; 140°9′E to 40°5′N; 141°7′E), where the urgent need for debris removal. The target area was selected to include at least 95% of the recent bottom trawling areas. This is because the government has listed the fishing area as a priority area for seabed debris removal operations.
To account the potential for significant volume variation depending on the specific debris type, the estimated volume of the debris accumulation by the model was normalized to represent "the probability of debris accumulations".

Relationships between debris and environmental variables
The volume of debris was 658.7 m 3 on average, with 82,711.2 m 3 at maximum, 121 m 3 at median per a 1 km grid from the survey data (Appendix 1); however, spatial variation of the volume was observed (Fig. 1). The volume of tsunami debris tended to decrease the further away from the coast (Appendices 2a, 3a), whereas more debris did not tend to be always sunk in deeper locations (Appendix 2b). Additionally, from the correlation analysis between volume (log-transformed) and each variable, environmental factors were not clearly related with the volume (Appendix 3).

Spatially varying coefficient model and its projection
Both oceanographic and land-based factors significantly influenced the volume of the debris (Fig. 3). The number of vessels nearest to the coast had the highest positive impact, followed by riverine area, tsunami height, and complexity of bathymetry, whereas the distance from the coast and kinetic energy influenced negatively (Fig. 3). A statistically significant effect in the urban area was not detected.
In our models, explanatory variables were assumed to have spatially varying coefficients (SVC). Five environmental variables among the influential variables were estimated to have SVC, except for the kinetic energy (Fig. 3). The distance from the coast had less impact on those in the Sendai bay as well as the range of depth and tsunami height (Fig. 4a-c). The number of vessels positively influenced the volume of debris at the northeast and inner areas of the Sendai bay, negatively affected the southern area (Fig. 4d). The effect of riverine area was greater around several estuaries (Fig. 4e).
The evaluation of prediction accuracy of the best model by RMSE showed that our model had lower RMSEs than the spatially constant model and null model, thereby indicating that our models had a better accuracy (RMSE = 2.43 (best model) 2.74 (spatially constant model), and 2.83 (null model)). The best model of the constant coefficient model showed the number of vessels and kinetic energy were significant variables, not tsunami impact, the range of depth, and proximity to the river and coast.
Estimation of the probability of the tsunami debris volume in our projection area showed that the volume decreased based on the distance from the coast. The probability of debris accumulation would be higher along the coast of Miyagi Prefecture than that of Iwate Prefecture (Fig. 5). Particularly, the Sendai bay was estimated to have a relatively high probability of debris accumulation, whereas the lower probability of accumulated debris was estimated outside Miyako (St. M).

Spatially varying factors on tsunami debris accumulation on the seafloor
Both land-born (number of vessels and riverine area) and oceanographic (kinetic energy, tsunami height, range of depth) factors affected the tsunami debris accumulation on seafloors as well as the coast proximity (Fig. 3). Furthermore, some of these factors had different coefficients across the studied area (Figs. 3, 4). Although tsunami debris accumulation on the seafloor is a socioeconomic and ecological concern that would affect fishery activities (Katsanevakis et al., 2007;Nash, 1992) and benthic organisms (Brown and Macfadyen, 2007;Chiappone et al., 2005), few studies have demonstrated the spatial explicit estimation of tsunami debris accumulation on the seafloor (Goto and Fig. 3. Results of spatially varying coefficient model for the total volume of tsunami debris. We show z value, which the mean value of the estimate divided by standard error of each environmental variable with an asterisk to denote statistical significance (flagged "**" for p < 0.01; "*" for p < 0.05; " ▪ " for p < 0.1). Each spatially varying coefficient of the variables was showed in the grey plot. Shibata, 2015). Our modelling approach that considered spatial structures of the data provided insights that multiple factors according to site-specific conditions affected the tsunami debris accumulation at the seabed.
The higher the tsunami run-up, the nearer from the coast, and the steeper the bathymetry have been proposed to increase the probability of the debris accumulation in the surrounding seafloor, which is consistent with that of previous reports (Lebreton and Borrero, 2013 for tsunami impact; Goto and Shibata, 2015 for the proximity to the coast; and Galgani et al., 1996;Pham et al., 2014 for the bathymetry). Tsunami backwash drags the debris on to the sea, with the force often positively correlated with tsunami height and negatively with the distance from the coast (Prasetya et al., 2012;Sugawara et al., 2009). The steep topography spreading in the eastern area of Miyagi Prefecture caused more debris accumulation (Fig. 2a) because the foot of the steep terrain acts as a crevasse within which more debris is piled up (Galgani et al., 1996;Pham et al., 2014). Hence, it is reasonable that the volume of debris responded to those factors. By estimating each coefficient, the effect of proximity to the coast was not detected around the Sendai bay (Fig. 4a) because of the effects of the Kuroshio current into the bay from the south (Fig. 2b). This flow prevented the debris from being washout to the offshore and might have counteracted the effect of distance attenuation. The bays with little inflow in this disaster area had a large amount of rubble in their reef after the tsunami backwash (Yokoyama et al., 2012), which corroborated with our interpretation.   Pham et al. (2014) showed that accumulation of marine debris was generally high in adjacent urban areas. By contrast, our study showed the effect of urban area coverage in the inundation zone was not detected, whereas the number of vessels and riverine area were observed. The urban areas in this study (St. K, St. M) were often adjacent to fishing ports with many vessels, which may also have resulted in various fishing materials contributing to tsunami debris (Fig. 2c, e, f).
Another land use factor, the riverine area coverage positively affected the accumulation. The effect of riverine coverage was stronger near Natori (St. N) and Ishinomaki (St. I) (Fig. 4e), which were near to watersheds of several rivers such as the Natori, Naruse, and Kitakami. The tsunami concentrated into the rivers and created vast inundated areas in Tohoku area (Kayane et al., 2014;Richmond et al., 2012;Yalciner et al., 2011;Yoshimura et al., 2014). Some field surveys after the disaster reported many tsunami debris and sediments through those rivers (Mikami and Suzuki, 2015;Ohta and Yamanaka, 2013). Our results were consistent with these observational reports that showed the significant influence of rivers on the outflow of debris into the ocean.
Based on the result of our model projection outside the coast of Tohoku, we confirmed that the accumulation would be larger in Miyagi Prefecture and smaller in Iwate Prefecture (Fig. 5) although the relatively higher tsunami run-up was observed in Iwate Prefecture (Fig. 2d). There are several explanations as to why the southern part of Tohoku area had more debris accumulation than the northern part. First, the inundated area was wider off the coast in Miyagi Prefecture ( Fig. 2f-g; bars in grey). It resulted in regional differences in the amount of debris generated by the disaster; 15.84 million tons in Miyagi and 4.52 million tons in Iwate as estimated by the Japan Ministry of Environment (Ministry of the Environment, 2012). These amounts reasonably reflect the difference in outflows into the ocean. Another possible reason is the difficulty of debris retention due to the high kinetic energy off Iwate (Fig. 2b). A previous simulation model for floating debris revealed that they would be quickly removed by a stronger velocity to convergent zones controlled by an ocean current (Maximenko et al., 2012). This tendency closely agree with that of previous studies, which revealed the sedimentation did not tend to occur in strong ocean flow (Onink et al., 2019;Yoon et al., 2010). Hence, the northern part of Tohoku area with relatively high kinetic energy, particularly off Miyako (St. M; Fig. 2b), had fewer debris accumulation, whereas the Sendai bay in the southern part (St. N and St. I) had more accumulation (Fig. 5). Considering that the windage of a debris may affect the distance from the coast to flow , the effect of kinetic energy still needs to be tested further along the gradient of buoyancy for debris types. In particular, it should be noted that ships have a relatively large volume and high buoyancy, making them easy to travel far away , but natural sediments which were not addressed in this study may respond differently to environmental factors. Furthermore, in the coastal area targeted in this study, it is possible that the influence of the backwash wave played a large role, but it should be kept in mind when predicting a more distant range from the coast.
Generally, debris was estimated to move further over time (Murray et al., 2018). However, our model was projected for one time slice after the disaster as a static estimation due to the need for a short-term cleanup operation during the early recovery period from the disaster. Under such urgent requirements soon after a disaster, a request for detailed information on debris based on model estimation often surpasses that of medium-and long-term predictions of debris accumulations. Temporal simulations of debris movements must be further explored but may be more related to mid-and long-term periods after a disaster where monitoring is desired to assess its impact on the ecosystem Kawamura et al., 2014).
In conclusion, given the potential limitation of human resources and budgets for seabed debris removal, prioritizing cleanup work should be planned. Although many predictions about the tsunami hazard area, few hazard maps assess tsunami debris accumulations on seafloors nor are prepared in advance for the future disaster (Bauer et al., 2008). Our model estimation provides information where the tsunami debris tends to be accumulated on the seafloor, supporting early reconstruction in the disaster area. Because data for these model variables are relatively common in other regions (for Japan; Administrative Management Bureau, Ministry of Internal Affairs and Communications https://www. data.go.jp/. last accessed 26 November 2019), our model may be applicable to other regions; however, it should be carefully tested due to different land use types in a target region. Future generalization across regions must provide a fundamental hazard assessment information for tsunami risk reduction in marine ecosystems and local communities.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.