Paleotsunami history of Hachinohe, northern Japan: a multiproxy analysis and numerical modeling approach

Paleotsunami studies along the Pacific coast of Tohoku, northern Japan, have been considerably developed recently, particularly after the massive impact of the 2011 Tohoku-oki tsunami. Nevertheless, in the southernmost Shimokita Peninsula, studies pertaining to paleotsunami are underdeveloped, leading to a vague understanding of the tsunamigenic sources northward of the Tohoku region, along with incomplete hazard evaluation. Paleotsunami deposits in Shimokita can be related not only to the Japan Trench along the Sanriku coast but also to the Kuril trench along the Pacific coast of Hokkaido. In this study, we unveiled the paleotsunami history of Hachinohe in northern Tohoku. Using a combination of sedimentological, geochemical, paleontological, and mineralogical proxies, we characterized seven sand layers that dated from ca. 2700 to ca. 5500 yr BP based on radiocarbon (14C) ages as event deposits of marine origin. Sedimentological and paleontological evidence coupled with ground-penetrating radar imagery revealed a marsh environment comprising successive extinct ponds, controlling the depositional environment. Numerical modeling ruled out the possibility of storms as genetic sources, leading to the conclusion that the presence of event deposits with marine sediments in the study area would be associated with tsunami inundation episodes. Based on 14C dating, the mean frequency of recurrence of tsunamis is estimated as 384 years (320–450 yr, 95% confidence interval) and a coefficient of variation of 0.78 (0.68–0.99, 95% confidence interval). The previously recorded limited paleotsunami evidence and absence of an estimated recurrence interval in the Shimokita Peninsula reaffirm the importance of Hachinohe as a tsunami record site for the activity of both trenches.


Paleotsunami studies in the Pacific coast of northern Japan
The Pacific coast of northern Japan is established as one of the regions worldwide with large and frequent earthquakes. Earthquake events occur not only in relation to the Japan Trench, facing the Tohoku region but also to the Kuril Trench, facing Hokkaido Island (Seno 1978;Ide and Aochi 2013). Historical events demonstrated the remarkable tsunamigenic capacity of the Japan Trench in 869 A.D., 1611A.D., , 1896A.D., , 1933A.D., , 1968A.D., , and 2011A.D., (Hatori 1975Tanioka and Satake 1996;Nagai et al. 2001;Ozawa et al. 2011;Uchida et al. 2016) (Fig. 1a). In the seventeenth century, a massive tsunami associated with the activity of the Kuril Trench exhibited a high capacity of inundation along the coast of Hokkaido, leaving deposits far inland and in high cliff areas (Ioki and Tanioka 2016). The magnitudes of these earthquakes have attracted the attention of the scientific community, particularly after the catastrophic 2011 Tohoku-oki tsunami (TOT), which clearly demonstrated that understanding tsunami hazard was not adequate to develop an effective disaster risk reduction strategy and prevent casualties due to tsunami episodes (Goto et al. 2021). To increase the knowledge regarding the frequency and magnitude of events from historical records, several studies evaluating the sedimentological evidence of tsunamis have been conducted along the Sanriku Coast (Minoura et al. 1994;Sugawara et al. 2012;Ishimura and Miyauchi 2015;Takada et al. 2016;Inoue et al. 2017;Ishizawa et al. 2018;Ishizawa et al. 2022) and the Pacific coast of Hokkaido (Nanayama et al. 2002(Nanayama et al. , 2003(Nanayama et al. , 2007Hirakawa et al. 2005;Ishizawa et al. 2017).
Regarding geological studies of paleotsunami deposits, clarifying the origin of tsunami, regional stratigraphic correlation, age, and source estimation is still a matter of discussion, not only in Japan but worldwide Costa et al. 2016;Jaffe et al. 2016;Ramírez-Herrera et al. 2016;Takada et al. 2016;Moreira et al. 2017;Ishizawa et al. 2018;Castillo-Aja et al. 2019). To the best of our knowledge, the only technique to obtain information regarding the spatial extent of tsunami deposits is consistent survey of coastal areas, with maximum possible resolution. The sedimentological characteristics of the tsunami deposits largely depend on the distribution and composition of sediments in the shore zone and topography on the backshore, thereby varying the geometry and lithology over short distances (Jagodziński et al. 2009;Abe et al. 2012;Moreira et al. 2017). Between Misawa in the Shimokita Peninsula (Tanigawa et al. 2014a) up to Hirono Town (Harashinai site), northern Sanriku Coast (Takada et al. 2016), there are no paleotsunami studies, rendering the region record incomplete (Fig. 1e). Southward, along with the Sanriku coast, several studies have been done (e.g., Imaizumi et al. 2010;Goto et al. 2015;Takada et al. 2016;Inoue et al. 2017;Ishimura 2017;Goto et al. 2019) (Fig. 1e), and the paleotsunami record is thorough. Due to the peculiar location and distance from the two neighboring subduction zones, the southern  1611, 1896, 1933. (Hatori 1975Tanioka and Satake 1996;Nagai et al. 2001;Ozawa et al. 2011;Ioki and Tanioka 2016;Uchida et al. 2016); To-a and To-Cu correspond to the isopachs of the respective tephra layers (Machida and Arai 2003). b Digital elevation model of the study area -5 m grid size provided by the Geospatial Information Authority of Japan. c Topographic map of the study area. Colored areas evidenced the levels of uplifted marine terraces. d Cross-section of the elevation model along the study area, from the A to A' . e Paleotsunami studies developed along the Sanriku Coast and the Shimokita Peninsula. The circles represent the approximate survey location for each study Page 3 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 Shimokita Peninsula helps in not only to recover the data on the maximum northward extent of tsunami events from the Japan Trench but also to obtain paleotsunami data that can elucidate the seismic activity of the Kuril Trench and the flexion between the two trenches.

Objective
This study contributes toward the understanding of tsunami hazard assessment in northern Japan by identifying, dating, and correlating paleotsunami deposits in Hachinohe. Our methodology was based on a combination of sedimentological, geochemical, geophysical, and paleontological tools, along with numerical simulations for inundation and sediment transport.

Brief review of approaches for paleotsunami deposit identification
Previous paleotsunami studies have demonstrated that characteristics such as upward fining, landward thinning and finning, erosional sedimentary structures, and paleocurrent directions are relevant indicators of the nature of the flow and key points to link deposit events to a tsunami source (Sugawara et al. 2008;Richmond et al. 2011;Cisternas et al. 2018). However, in several cases, such characteristics can also be the result of river floods (Inoue et al. 2017;Shyu et al. 2019). By performing grain size analysis, it is possible to infer some of these textural characteristics and contribute toward the understanding of the sedimentary environment, energy, and transport media involved during its deposition (López 2017). In addition, the characterization of the sedimentary fabric can be complemented by identifying the relative density changes generated by lithological variations (Kain et al. 2015). X-ray computed tomography (CT) can illustrate the compositional variations in the sedimentary sequence because it records the attenuation of X-rays on materials due to the Compton scattering effect, the magnitude of which depends primarily on density and consequently on the mineralogy (Mees et al. 2003). For instance, Paris et al. (2019) utilized high-resolution micro-X-ray CT to characterize the sedimentary fabric, grading, and imbrication angles, inferring tsunami as the transport agent. In addition, geochemical characterization of paleotsunami deposits can be performed by determining relative quantities of major elements, trace elements, and carbonates (Chagué- . Because of the leaching effect of rainwater and groundwater on high-mobility elements, only major and trace elements must be utilized to characterize deposits with long-term diagenesis, that is, paleotsunami deposits Shinozaki et al. 2015;Chagué-Goff et al. 2017;Judd et al. 2017), However, because the characterized group of elements can possess different mineralogical arrangements and hence different geological meanings, it is imperative to reinforce element interpretation with mineralogical identification techniques such as X-ray diffraction (XRD) (Croudace and Rothwell 2015). Additionally, XRF is beneficial for relative quantification of the organic matter content (Croudace et al. 2006;Chawchai et al. 2016), as demonstrated in previous tsunami deposit studies (Chagué-Goff et al. 2016;Judd et al. 2017). The obtained data have been analyzed using multivariate statistical methods, such as principal component analysis (PCA), to find associations or linear dependencies among the measured elements and clarify the sedimentary source (Chagué- Goff et al. 2016). To support geological data, few studies have used diatom analysis, and shell fragment recognition to clarify the terrestrial or marine origin of event deposits (Sawai et al. 2009;Tanigawa et al. 2014b;Chagué-Goff et al. 2015;Judd et al. 2017;Cisternas et al. 2018).
Regarding the differentiation of the origin of the deposit between storms or tsunamis (Morton et al. 2007;Sugawara et al. 2008;, recent studies have demonstrated that numerical modeling can be suitable for solving this problem by comparing the inundation capacity for both phenomena (Inoue et al. 2017;Watanabe et al. 2018).

Study area
The study area is in the northernmost part of the coastal zone of Hachinohe City, Aomori Prefecture, Japan, surrounded by the Oirase and Gonohe rivers. A succession of uplifted both Pleistocene and Holocene marine and river terraces resulting from the constant tectonic uplift related to the subduction in the Japan Trench (Miyauchi 1987;Niwa and Sugai 2021) controls the stepwise morphology of the area. Based on digital elevation models (DEMs), it is possible to infer the presence of the four different surfaces (Fig. 1c). The lowest surface (1) is about 3 m high; a following ~ 4 m high surface (2) is edged by the past lateral migration of the surrounding rivers; next, a ~ 5 m high surface (3) (Fig. 1d). These three surfaces have been named "Holocene Lowland" by Miyauchi (1987), and include river terraces, backshore deposits, and poorly developed marine terraces. Lastly, a 35 m high surface (4) corresponds to the marine terrace named Takadate Surface by Miyauchi (1987). The survey area is located on the surface (3).
Land use mostly includes industrial activities and housing on surfaces 1 and 2, paddy fields on surface 3, and housing on surface 4. The presence of 5 m high marine-surge protection structures on the beach berm and river dikes is remarkable in the current topography. In the northern part of the 4 m surface, next to the Oirase river, a j-shaped depression is present that can be Page 4 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 morphologically related both to a cutoff meander and a paleo-pond (Fig. 1c). As a first approach to paleotsunami research in Hachinohe, Velasco et al. (2019) conducted numerical modeling for storm and tsunami surges to compare their inundation propagation. To reproduce natural coastal settings, artificial features, such as protection walls and breakwaters, were removed from the DEM used in this study. Extreme and improbable storm parameters were assumed for the storm simulations-physical conditions of the 1979 Typhoon Tip, with a pressure of 870 hPa, and constant winds of 140 knots for 10 min; and the propagation speed was taken from the 2013 Haiyan Typhoon (Watanabe et al. 2018). For the tsunami simulation, the source parameters of the 2011 TOT were applied. The results showed that even under the most severe storm conditions, erosion and sedimentation are limited to areas below terrace 2, under the assumption of the present sea level (Fig. 1). In addition, tsunami surges derived from the 2011 TOT can inundate further inland and reach the edge of terrace 3 (Fig. 2, Velasco et al. (2019)). In practical terms, this means that up to a height of 3 m, the sedimentary environment can be affected by the occasional occurrence of extreme events such as river floods, storms, and tsunamis. Moreover, in areas higher than 3 m, the sedimentary environment can solely be affected by either tsunamis or river floods, restraining the differentiation of event deposits to a marine or terrestrial origin. As the modeled tsunami surge reached the terrace 3 edge, we selected a survey transect along such a surface (Fig. 2c).

Ground penetrating radar (GPR) survey and excavation using Handy Geoslicer
As reported by Velasco et al. (2019), who performed a preliminary site survey, GPR lines were used to select the best possible survey points by identifying paleosurfaces with the maximum sediment preservation potential. This study used processed GPR imagery for paleotopographic interpretation.
We used a GPR controller GSSI SIR ® -4000 (Geophysical Survey Systems, Inc., USA) with an antenna set up to a frequency of 400 MHz to achieve a maximum effective   0  200  400  600  800  1000  1200  1400  1600  1800 Distance ( Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 vertical penetration of 3 m and a maximum resolution of 18 cm (1/4 of wavelength). GPR data were processed using the RADAN 7 software by GSSI (Geophysical Survey Systems, Inc., USA). A Handy Geoslicer (Takada et al. 2002) was used for core sampling. Sediment core samples were extracted along a transect perpendicular to the shoreline.
3.2 Non-destructive core scanning: X-ray CT X-ray CT was used to identify and clarify lithological changes, sedimentary structures, and sediment grading. X-ray CT was conducted using the scanner Toshiba Aquilion Prime at the Center for Advanced Marine Core Research facilities at Kochi University. Horos software (open source licensed under the GNU Lesser General Public License, Version 3.0 (LGPL 3.0)) was used to generate 3D reconstructions, from axial-perpendicular slides obtained every 0.5 mm. In addition, we modified the colors depending on the CT units to enhance the visualization of lithological changes. Image J software (public domain Java image processing and analysis program) extracted a CT value profile from the front-axial images to be included in the log worksheet along with the X-ray fluorescence scanning results.
In the 3D reconstruction, the original grayscale was modified to blue, green, and yellow green scales to facilitate the recognition of the lithological features. These colors were manually selected on Horos for convenience to differentiate the lithological changes. Threshold values were determined by first identifying sand, peat, tephra layers, and roots and further manually assigning such values, thereby allowing faster identification of composite materials.

Non-destructive core scanning: magnetic susceptibility
Magnetic susceptibility (MS) was measured to identify lateral mineralogical changes linked to the tsunami inundation trajectory. As MS values are absolute quantities, we utilized them for lateral comparison among the measured cores. The MS of the sediment core samples was measured using the core logging sensor (Model MS2C, Bartington, USA) with a resolution of 10 mm.
The obtained values were plotted as a vertical profile and included in the log worksheet along with the XRF profiles.

Non-destructive core scanning: micro XRF (ITRAX)
We utilized XRF core-scanning by ITRAX (Löwemark et al. 2019) to obtain the content of the major elements and identify geochemical signatures that contribute to the source discrimination. The equipment was facilitated by the Center for Advanced Marine Core Research at Kochi University. The measurements were conducted using a molybdenum tube with acquisition parameters of 30 kV, 55 mA, exposure time of 10 s, and step size of 1 mm. XRF raw data were further processed to identify points with unreliable measured values due to the irregularities in the sample or misleading measurements. This processing used the "validity profile" delivered along with the raw data by the ITRAX equipment and discarded unreliable measurements. Then, elemental spectra were normalized by the total count rate (kcps) to correct the effect related to sample quality and measurement parameters (Chagué- Goff et al. 2016). Also, we calculated centeredlog ratios for each element measurement (clr) to improve the precision in correlations (Löwemark et al. 2011;Mondal et al. 2021) using the following equation: where I jj is the intensity of element i after normalization at depth j, and gm j is the geometric mean at depth j of the raw data after normalization by kcps.
PCA and linear discriminant analysis (LDA) were utilized to identify geochemical changes among the sand layers and their provenance. PCA has been primarily used in tsunami sedimentology to identify relationships or discriminate among sand layers (Trauth 2015;Chagué-Goff et al. 2016). Herein, PCA was applied to identify geochemical changes in the same sediment core. Based on the XRF-element detection rates (Bishop 2016) and pattern profiles, we used Sr, Ca, Mn, Ti, Fe, Cr, Ta, Zr, Si, K, Al, Rb, and inc/coh as inputs for the PCA. As our interest lies in the geochemical characterization of the sand layers, we extracted the data corresponding to every sand layer. For the data processing we used PCA code provided by Nell (2017) LDA was utilized for sand correlation among the core samples. LDA and Fischer discriminants are valuable statistical techniques for classification and reducing components when inferring a relationship among them. Such a linear function maximizes the distance between two classes and minimizes the variation between each class (Govindaraju and Rao 2014). Herein, LDA determined the relationships among the sand layers in every core sample. To perform such analysis, we selected XRF-element detection rates restricted to the sand layers in every core and compared them using the data visualization tool DisplayR ® . In XRF core scanning, because the quantification of major elements in counts per second (cps) largely depends on the measuring parameters and sample conditions (Chagué- Goff et al. 2016;Löwemark et al. 2019), it is not possible to conduct direct associations of elemental values among different cores, even after kcps normalization. As a ratio is a quantitative relation between two elements of the same core (Croudace and Rothwell 2015), it does not depend on analytical factors and can be utilized to compare with ratios among other samples. To conduct such a comparison along with every sand-event sheet, we applied LDA through the predominant elements in PCA, normalized by Ta (for the use of Ta as normalization element, see the discussion in the body text).

Grain size analysis
Grain size analysis was performed on the sand layers to identify vertical and lateral grading that could be used to identify marine-related surges. Each sand layer was divided into slices of 1 cm thick. A cubic sample of 1.5 cm wide and 2 cm depth was extracted from every slice, with 3.5 g of dry sample on average. Prior to the grain size analysis, samples were treated to remove organic matter by dissolution in 10% hydrogen peroxide (H 2 O 2 ), desiccation at 60 °C for 12 h, and homogenization by mixing with a spatula for 1 min. Subsequently, we utilized a sample splitter to obtain a suitable amount of sediment for the laser equipment. The sample was divided in two half-lots three times to obtain approximately 0.5 g of dry sample. Subsequently, we utilized the sediment sieving method to separate particles > 2000 μm to adequately characterize smaller particles with laser diffraction (Gerlach et al. 2002). We utilized a laser diffraction analyzer (Shimadzu SALD-2300, Japan) at Tohoku University. The resolution of the equipment ranges from 17 nm to 2500 μm. Grain size results were processed in GRADISTAT v8 (Blott and Pye 2001) to obtain the grain size distribution, mean, sorting, kurtosis, and skewness. GRADISTAT implements the graphical classification methods of moment and Folk and Ward.

XRD analysis
XRD analysis was conducted to identify the mineralogical source of key elements, such as Ca, identified in the XRF analysis and the presence of glauconite in the deposit. A total of 5 g was extracted from each sand layer for homogenization. Sampling was conducted throughout the entire layer to obtain a representative lithological composition. XRD analyses were conducted on the powder and oriented samples extracted from the sand layers. To separate and orientate the clay-size portion from the sediment sand layers, we followed the sample preparation technique proposed by Moore and Reynolds (1999). First, we disaggregated the unconsolidated minerals and homogenized the samples with a bead cell disrupter. Afterward, 2 g of the sample was dispersed by ultrasound in pure water. After the decantation of the heavier minerals, the suspended material was centrifuged and filtered using Millipore ™ . Finally, the clay-portion was transferred to the glass slide.
The oriented samples were treated and analyzed in three stages: air-dried, ethylene glycol solvated for smectite identification, and heated to 550 °C for 3 h for kaolinite identification. We utilized the diffractometer from the Earth Science Department of Tohoku University (Philips X'pert Pro MPD, USA). The measuring resolution was 0.5 min every degree, from 0° to 35° 2θ. The resulting patterns were processed and interpreted using the Match! software (Crystal Impact GbR, Germany).
Shell fragment recognition was conducted with two objectives: evidence of marine influence and high-energy surges. The shell fragments were separated and identified using a trinocular microscope.

Carbonate quantification
Carbonate quantification aimed to quantify biogenic carbonate and complement elemental and mineralogical interpretations based on XRF core scanning and XRD, respectively. Carbonate content was calculated by subtracting the initial and final weights after controlling pH dissolution. To fix the pH to 5.2, we used a solution of 500 ml of 2% acetic acid (CH 3 COOH) and 6.5 g of ammonium acetate (CH 3 COONH 4 ). Each sample (200 g) was poured into 100 ml of solution for 24 h. The sample was further dehydrated using a 0.8 µm Millipore ™ filter and desiccated at 24 °C for 24 h.

Tephra analysis
Tephra analysis aimed to identify the layers T1 and T6 to constrain the radiocarbon dating modeling. Samples were extracted from the tephra layers for analysis. Kyoto Fission-Track Co., Ltd. conducted the analyses. The petrological characteristics of the tephra layers were obtained via mineralogy analysis by microscope, selecting up to 200 grains of volcanic glass, light minerals, heavy minerals, and rock fragments. Characteristics such as heavy mineral content, volcanic glass shape, and refractive index of glass shards, hornblende, and orthopyroxene crystals were analyzed. For refractive index in T1, 81 volcanic glass shards and 60 orthopyroxene shards were measured. For T6, 66 volcanic glass shards and 60 orthopyroxene shards were measured. The obtained data were further correlated with standardized values of known tephra layers in Japan.
Herein, the term "tephra" is defined as any deposit event conformed mainly by volcanic materials, from clay up to pebbles, that are either possible fallout tephra layers or reworked tephra layers. The source and age of the layers can be either known or unknown.

Radiocarbon dating
Bulk samples were obtained from the interbedding peatymud layers for 14 C dating. Slices of 1 cm wide and 5 mm thick were extracted. The samples were treated with an acid wash to remove the bone collagen and cellulose. Furthermore, 14 C dating was conducted using accelerator mass spectrometry (AMS) in Beta Analytic Inc.
We employed OxCal 4.4.4 (Bronk Ramsey 2021) to calibrate dating results to calendar ages and applied the calibration curve dataset IntCal20 (Reimer et al. 2020). The sequence model (Lienkaemper and Ramsey 2009) was used to calculate the estimated deposition age of the sand layers (Lienkaemper and Ramsey 2009). In addition, using OxCal, it was possible to obtain probabilistic recurrence intervals for deposition ages for every sand unit using the "Difference ()" function.

GPR, core sampling, and sedimentary features
Thirteen GPR lines were collected from the survey area (Fig. 1b). Thirteen core samples (H01-H13) were obtained. Core sample Ha8 was obtained in a preliminary recognition of the area at the same point as H02 (Shinohara et al. 2017). These core samples were described in the field and analyzed in the laboratory.
GPR profile images revealed two major features regarding the stratigraphic configuration ( Fig. 3a): first, small basins associated with large depths of the paleosurfaces and thus locations with better preservation potential of sediments; second, the presence of two strong reflectors in the sedimentary sequence.
The relative positions of the core samples are illustrated in Fig. 3b. The general sedimentological characteristics of the sand layers found in core samples Ha8, H02, H03, H04, and H06 are described in Table 1. As stated in the Sect. 3, these cores were later used for the laboratory analyses. Core samples Ha8 and H02 contained the largest number of sand layers among the core samples taken (Fig. 3c), with ten layers. Core samples H03, H04, and H06 contained seven, six, and two sand layers, respectively. The sedimentary record generally displayed a succession of two or three tephra layers at the top; the following section comprised intercalation among sand layers and peaty-mud, with occasional ashy tephras and a thick tephra layer at the bottom.
As core sample H02 contained the most representative sedimentary succession ( Fig. 3b and c), it can be described as follows: the upper part of the core sample H02 was characterized by cultivated soil, followed by a succession of three thick tephra layers (T1, T2, and T3) along with two interbedding thin mud beds. The underlying sequence comprised the interbedding of sand and peaty-mud, along with two tephra laminae (T4 and T5). Within this section, it was possible to recognize eight sand layers with the naked eye, named as units S1, S2, S4, S5, S7, S8, S9, and S10 from top to bottom. Units S3 and S6 in H02 were recognized by a magnifying glass, and their boundaries were later constrained by XRF core scanning and CT images (see the Sect. 5 in the body text). The bottommost section was formed by a thick tephra layer (T6).

CT, MS, geochemistry, and grain size analysis
X-ray CT, MS, XRF, and grain size analyses were conducted on core samples H02, H03, H04, and H06.
The CT image reconstruction of core samples H02, H03, and H04 is shown in Fig. 4.
In CT reconstruction, it is possible to recognize the sand units (yellow green in CT color), interbedded nonbioturbated mud (blue in CT color), and tephra layers (see color palette for CT in Fig. 4). Sedimentological characteristics such as geological contacts (planar, erosive, or transitional), intraclasts, and bioturbation, of core H02 have also been delineated and shown in Fig. 4. Sand layers S3 and S6 were highlighted in the CT image, which allowed the boundaries to be constrained. The mud layer that separate units S8 and S7 shows high mixing with sand aggregates indicating bioturbation. Also, it is noticeable there is a considerable increase in the intensity of bioturbation, as evidenced by the heterogeneity of the sand layers. Units S2 and S3 are cut by roots (Fig. 4).
Page 8 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 High values of MS in core H02 can be linked to the presence of sand layers S4 to S10 (Figs. 5 and 6).
XRF elemental patterns obtained from core H02 exhibited conspicuous changes related to the interface between the peaty-mud and the sand layers. Si, Ca, Sr, K, Ti, and Fe had higher values on sands, whereas Ta exhibited lower values (Fig. 6). Elemental profiles also contributed to differentiating S3 and S6 from the neighboring mud layers. Fluctuation pattern of Sr was relatively like that of Ca. Ti exhibited the most remarkable changes among the elements measured; every sand layer can be characterized by an abrupt change in Ti across, over, and underlying mud beds. However, high Ta values were related to peaty-mud layers. A ratio of Ca to Ta was introduced to clearly delineate the interface between sand and peaty mud. High Mo inc/coh ratio values represented the high relative content of organic matter in the peaty-mud layers. Although sand layers were characterized by low Mo inc/coh intensity, the pattern behaved oppositely to Ca, Ti, and Ca/Ta. We applied PCA for the XRF data obtained on the cores H02, H03, H04, and H06.
The CT value profile extracted from the CT images ( Fig. 6) displayed high values related to the sand layers. As stated before, the CT images primarily reflect density changes; thus, this profile illustrates those related to lithological changes among sand, mud, and tephra.
In a first-order difference, vertical variations in MS (Fig. 5), elemental composition ( Fig. 6), and CT color intensity ( Fig. 4) were linked to the changes in mineral composition and density; at the second-order difference to grain size and compaction (Fig. 7).
Grain size analysis was conducted in cores H02, H03, and H04 to identify vertical and landward changes in each sample. Some sand layers (e.g., S1, S2, and S3 in H02) were not sliced and analyzed by grain size owing to the nominal thickness of the layer. In addition, sand layers in core sample H06 were not analyzed because of the high bioturbation and mixing with peaty mud. Regarding the sample processing, although it is expected that the sample splitter produces a minor overall inaccuracy level among the splitting methods, some segregation could have occurred because of the number of times the sample   Fig. 3 a Location of core sample H02/Ha8 and relative positions of core samples H01 and H13on the GPR profile. Note that core samples H01 and H13 are displaced laterally from the GPR line (Fig. 1b). b Relative position of horizontal core sample and sand events correlation, based on sedimentary features and stratigraphic relationships. c Stratigraphic column of core samples Ha8 and H02, separated laterally by ± 2 m Page 9 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19   was divided; nonetheless, it is expected that the error level is not greater than 2% (Gerlach et al. 2002). Grain size analysis revealed that the sand layers are mostly composed of sand with grain sizes of 1.5-2.0 phi, which frequently exhibit upward fining (Fig. 7). Thicker sand layers, such as S8b and S5 at H02, typically exhibited normal or inverse grading. In addition, the sand layers exhibited vertical changes in mode values (Fig. 7). Some curves exhibited multimodal grain-size distribution, primarily related to sand-mud mixing or sand-granule mixing. As observed in S4, S5, and S8, in addition to grading, sorting and skewness exhibited a vertical variation in the sand layers (Fig. 7).
In Table 2, the specific characteristics of each sand layer with respect to geochemistry, MS, CT, and grain size are described. Not all features are observed in every proxy.

PCA and linear component analysis
The first three principal components represented 59% of the total variability, with distributions of 39%, 11%, and 9%, respectively. After using several combinations of elements, it was not possible to separate the sand layers in the clusters (Fig. 8). Nonetheless, although the cloud of data is averaged to the center, the data illustrated the tendency of unit S8, S9 and S10 to move from  Page 13 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 Sr, Al, and Ca. Units S1, S2, and S3 were gently inclined to be associated with Rb, Cr, Ta, and Zr. We applied LDA to the predominant elements in PCA (Si, Ca, Cr, Zr, Cr, Ti, and Fe), normalized by Ta (see the discussion in the body text), which behaved opposite to sand-related elements in the XRF results. In essence, the LDA scatterplot shows every sand layer's mean value and its relative position to the elemental ratios. It helps to see how chemically "similar" one sand layer to another is and how correlatable the sand layers are among the core samples (Fig. 9). LDA results show the data is mainly distributed from ~ − 3.0 to ~ 2.0 in the first linear discriminant and from ~ − 2.0 to 4.0 in the second linear discriminant. H02_S1 and H02_S2 show no correlation with the other sand layers. Independent and explicit clustering of sand layers is not possible among the remaining sand layers. However, although it is possible to make some close associations among S3, S4, S5, S6, and S7, there are no clear and unique geochemical signatures for these.
Even after the separation of clay portion, plagioclase (varying compositionally from albite to andesine) and quartz were highly representative in every sample.

Paleontological analyses
Diatom analysis was conducted on 24 sub-samples obtained from the core sample Ha8. Diatom species related to freshwater (34 species), freshwater-brackish (3 species), and brackish-marine (5 species) environments were identified in core sample Ha8 (Fig. 10). Both freshwater and freshwater-brackish species were found in sand layers, peaty mud, and tephra layers. Brackishmarine (B-M) species were identified in S5, S7, S8, S9 and S10 (Table 2 and Fig. 10). The proportional amount of the B-M species with respect to the total found species was < 10% in every sample of every sand layer.
Using a trinocular microscope, a shell fragment of the bivalve was recognized in unit S5 (Fig. 11). In the other sand layers, it was not possible to find any mollusk fragments.

Tephra analysis
Tephra layers at the top and base of the core samples were sampled in Ha8 for identification.
T1 is conformed in general by two layers. The basal unit corresponds to a 1 cm thick ash layer overlayed by a chaotic massive tephra deposit containing volcanic material ranging from clay to lapilli. T6 is formed by a chaotic deposit of volcanic material that ranges from clay to lapilli.
Tephra analysis revealed that both tephra layers are related to the eruptive events of the Towada Volcano. Based on the petrological analysis (mineralogy, heavy mineral content, volcanic glass shape, and refractive index of glass shards, hornblende, and orthopyroxene crystals), T1 was identified as To-a, whose age is 915 CE    (Machida and Arai 2003). In addition, basal tephra was associated with To-Cu, with a deposition age of 5986-5899 cal yr BP (McLean et al. 2018). Reference values were taken from Machida and Arai (2003) (Table 3).

14 C dating
A total of 11 samples were obtained from core sample Ha8 for the 14 C dating. Three samples from H04 and two from H02 were utilized for the lateral correlation. The conventional radiocarbon ages obtained from core samples H02, Ha8, and H04 are listed in Table 4. The calibrated results are shown in Fig. 12.

Sand events correlation
Based on the GPR imagethe two strong reflectors identified were interpreted from top to base as the contacts between T3 and the underlying mud layer, and between T6 and the overlaying mud layer (Fig. 3a). Further, the initial lateral correlation of the sand-peaty mud sequence was conducted based on sedimentary features, color, and stratigraphic position (Fig. 3c). The persistent presence of tephra layer four (T4) among several cores served as a marker horizon (Fig. 3b, c; Table 5). Table 6 presents the lateral correlations among the core samples. The correlation was supported with CT, geochemical analysis, and radiocarbon dating.

Correlation between core samples H02 and Ha8
As the core sample Ha8 was the first core to be extracted, diatom analysis, tephra analysis, and radiocarbon dating for event recurrence were conducted on this core sample. When conducting the field survey, the extraction of the core sample was intended to be in a 2 m radius from the Ha8 sampling site, which corresponded to the GPS device horizontal error: ± 2 m, to assure a good correlation between them.
The correlation between Ha8 and H02 was determined based on the stratigraphic position of the sand layers, lithology, and sedimentary characteristics (Figs. 3, 13 Fig. 6 Log pane of the core sample H02. Log distribution: depth, stratigraphic units, CT image, validity, CT value profile, magnetic susceptibility, and micro X-ray fluorescence elements result (Si, Ca, Sr, K, Ti, Fe, Ta, Ca/Ta, and Inc/Coh) vs depth. Black arrows indicate trend changes in pattern intensity Page 15 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 • T1 is present in all core samples and was correlated using stratigraphic and lithological characteristics. It shows a basal ash layer (T1p) and an overlying chaotic deposit of volcanic material. (Fig. 3c). T1p is considered as primary tephra based on the following results: tephra identification results recognized this layer as the To-a event (915 CE); no diatom species were found that could show reworking; 14C results show the age below T1 corresponds to 1184 to 1271 cal BP (Fig. 13), which means this mud layer was deposited before the To-a; and, it fits the reported isopach map by Machida and Arai (2003), which reports a thickness < 5 cm (Fig. 1) for the To-a event. Consequently, the overlying chaotic deposit is considered the possible secondary deposit resulting from a lahar pyroclastic flow derived from the Oirase River, which directly drains towards the Towada caldera. On the other hand, freshwater diatom species are mixed only at the top part, which supports the interpretation of the deposit source (Fig. 10). • T2 and T3 were correlated between the two core samples using the lithological characteristics and stratigraphic position. These are considered secondary deposits by flooding events related to the Oirase river, supported by the layering with mud in T3 (Table 3, Fig. 3c), and the inclusion of freshwater diatom species into the deposit (Fig. 10) which points out the reworking. • Sand layers S1, S2, S3, and S4, were correlated using the relative stratigraphic position because the lithological characteristics were identical (Fig. 13). • Based on dating results, the mud layer overlying S5 in Ha8 -sample C5 can be directly correlated to the one overlying S5 in H02 -sample C12, with ages ranging from 2973 to 3075 cal BP and 2800 to 3000 cal BP, respectively (Fig. 13). • T4 was used as a marker horizon among the core samples (Fig. 3b), based on its thickness and clayey texture. However, T4 is considered a possible fallout tephra layer or reworked tephra layer which source and age are unknown. • S5 was above T4 for both Ha8 and H02. S6 was below T4, and the sedimentary characteristics were identical.  Table 2 Summary of sedimentological, paleontological, and geochemical analysis of the event sand layers (S1 to S10) For lateral correlation of the sand layers S1 to S10, see Sect. 5 in the body text Sand unit Grain size analysis  Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 • The dating result of the mud layer separating S8 and S9 in Ha8 -sample C9 did not coincide with the same stratigraphic position in H02 -sample C13 (Fig. 13). Nonetheless, the age difference is small, which could be related to the close time frame or reworking of sedimentary events during sand event deposition. However, it supports the correlation of S8 as it was deposited after 4892 cal BP, and the sedimentary characteristics of the layers coincide. • T5 is also considered a possible fallout tephra layer or reworked tephra layer which source and age are unknown. • In both cores, both S9 and S10 appear as thin beds underlying S8 (Fig. 13).
• T6 is present in all core samples and was correlated using stratigraphic and lithological characteristics. There are no clear insights about the primary or secondary character of the tephra layer. According to the isopach map by Machida and Arai (2003), To-Cu should be < 10 cm; however, the thickness of the observed layer is > 40 cm. We consider this thick layer was generated by an initial ash fallout and an ensuing lahar or pyroclastic flow derived from the Oirase River.
In general, the thickness of the sand layers changed considerably between the two core samples. S3 was not recognizable in Ha8, and S6 was significantly reduced in thickness to a medium lamina (Fig. 3b). Such strong lateral changes in thickness can occur in tsunami deposits (Nakamura et al. 2012).

Lateral correlation of cores H02, H03, H04 and H06 using LDA and 14 C
The XRF results and its analysis by LDA contributed to clearly differentiating the sedimentary facies of S1 and S2 from the other sand layers. However, LDA evidenced a substantial geochemical similarity among the sands S4, S5, S6, S7, S8, S9, and S10 in different cores, suggesting a similar lithological composition, and thus the source of the sediments. However, although it was possible to cluster the sand layers, some samples are close to different sand layers, which means that in the case of performing clustering by using automatic algorithms or unsupervised techniques, wrong associations can be done. For instance, S5 in core H06 lacks a good correlation with the corresponding depth in H02, H03, and H04 (Fig. 4), implying that the furthest landward sand could correspond to a different sedimentary event and source. Therefore, the clustering must be made carefully to avoid the false correlation of sand layers (Fig. 9). As for the correlation between Ha8 and H02, the 14 C dating results supported the lateral correlation of the sand layers (Fig. 13). In addition, the mud layer overlying S10 in Ha8 can be correlated to the same stratigraphic layer in H04, with age ranges of 4892-5034 cal BP and 4854-4960 cal BP, respectively. Calibrated ages show that the peaty-mud layers underlying S5 in Ha8 and H04 are of similar age. On the other hand, the dating result from the peat-mud underlying S7 in H04 does not correspond to the age obtained on the same stratigraphic level on Ha8, which can be caused by modern contamination during core extraction or sampling. As stated previously, differences in dating ages among core samples H02 and Ha8 may be related to the erroneously sampling of reworked material (Ishizawa et al. 2018), or sampling resolution (Fig. 13) Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 We followed the same methodology using the XRF data and the LDA for the lateral correlation in the remaining cores (i.e., H13, H01, H07, H08, H09, H10, H11, H12, and H05) (Fig. 3b).

Identification of tsunami origin
We consider brackish-marine diatom species, shell fragments, glauconite occurrence, landward thinning, and landward reduction of MS value as strong evidence of marine provenance of the sand layer (Table 2). Regarding marine provenance, a factual imprint is the presence of fragments or complete structures of marine organisms, which is, in our case, brackish-marine diatom species and shell fragments. Complementarily, glauconite also supports the marine origin of the sand layers by usually being an authigenic mineral of shallow oceanic platform environments (Huggett 2013). However, it can also be related to reworked material from the Miocene or the Pleistocene deposits (Kamada et al. 1991;Phillips et al. 2017).
Sand units S4, S5, S6, S7, S8, S9, and S10 can be considered as clear marine-origin events, based on attributes such as landward thinning, erosive basal contacts (Fig. 4), landward decrease in magnetic susceptibility value -indicating landward mineral gradation during the wave progression (Fig. 5) Fig. 10 Diatom analysis result for the core sample Ha 8. In the horizontal axis are described the diatom species found along with the sediment core: freshwater, freshwater-brackish, and brackish-marine environments. On the right side, diatom-origin proportion is described as a percentage. A close view depicts the presence of marine species (from 0 to 30%) 0 1 mm H02_S4 Fig. 11 Shell fragment recovered from S5 Page 19 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19   Page 20 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 brackish-marine species (Fig. 10), upward fining, and glauconite occurrence (Table 2). Complementarily, the sand layers can be categorized as tsunami deposits when using the results from numerical modeling that suggest that only tsunami events can inundate the 5 m-uplifted terrace. As mentioned before, Velasco et al. (2019) compared the inundation capacity of both storms and tsunamis in the area of this study by means of numerical modeling using the present-day topography without coastal engineered dikes. For the storm calculation, extreme parameters were utilized for a cyclone, which is unrealistic for Hachinohe because the trajectory of the storm was NE to SW (counter-clockwise), opposite to the expected Coriolis-induced SW to NE direction. Even under such extreme boundary conditions, the storm surge could not inundate and deposit sediments on areas higher than 3 m in elevation (Fig. 2a). For the tsunami surge, under the source parameters of the 2011 Tohoku-oki event, the calculated wave could reach the 5 m-uplifted terrace (Fig. 2b). Therefore, closer, or more potent tsunami sources can inundate the study area further inland.
According to sea-level reconstructions, sand layers S4 to S10 were deposited under a relative sea-level highstand about 2 m higher than the present day's (i.e., from ~ 6000 cal BP to ~ 2800 cal BP) (Niwa and Sugai 2021). It means, the edge of the current five-meterterrace should be three-meter high for that time. Even though we did not perform reconstructions for the paleo shoreline, such maximum sea level implies that even in the worst storm scenario, its inundation capacity should still be limited by the topography. Hence, the sand units identified as a marine origin (i.e., S4, S5, S6, S7, S8, S9, and S10) were deposited on the exclusive tsunami domain zone (Figs. 3, 4).
Contrastingly, sand units S1, S2, and S3 lack any such evidence; they are not related to marine-origin events; hence cannot be considered tsunami events. Although the lack of marine evidence does not necessarily indicate a non-marine origin , its association with a tsunami event is weaker than that of the other sand units. Suppose a marine origin is assumed for sand units S1, S2, and S3. In this case, one possible explanation is that they were deposited close to the limit of the maximum sediment extent, where the flow conditions only allow reduced sediment content, typically represented in mica flakes and clay minerals (Yoshii et al. 2017). Consequently, the remaining thin and extentlimited sand layers were highly susceptible to leaching, weathering, and poor preservation of marine-related evidence. This hypothesis is supported by Niwa and Sugai (2021). They showed that there was a slight sea-level fall from ~ 2800 cal BP in the Hachinohe area, meaning that the inundation capacity of tsunami events was progressively reduced.
Another possible explanation is that the sand layers S1, S2, and S3 can be associated with possible river floods, as the survey area is bounded by the Oirase and Gonohe rivers.5.3 Sedimentation processes.

Implications from XRF, XRD and CT data
The CT profile displayed bulk density changes associated with landward thickness reduction, packing, and grading (Figs. 4 and 7). Likewise, XRF-pattern analysis is significant not only for geochemical characterization but also for sedimentological analysis by showing elemental changes related to the mud-sand ratio change ( Fig. 6 and Additional file 3: Appendix 3, Additional file 4: Appendix 4, and Additional file 5: Appendix 5). Such trend changes in the sand layers exhibit a close relationship between the elemental pattern and changes in packing, and composition (Figs. 4, 6, and 7 and Additional file 3: Appendix 3, Additional file 4: Appendix 4, and Additional file 5: Appendix 5). As the mud is mixed into the sand layer, the clastic signal (i.e., sand fraction) decreases with sorting and composition change (Fig. 7). Thus, the upward trend of the decrease or increase in Ca and Sr patterns in the sand layers can be associated with the sand-mud ratio change, respectively. As these two elements (Ca and Sr) have similar atomic properties, they can easily replace each other in the lattice (Nichols 2009) and can be related to the same origin. In addition, changes in Ti display variations in the content of heavy minerals (Croudace and Rothwell 2015) within sandy deposits, reflecting the density-gradation by a decrease in the particle settling velocity, due to the landward change in MS.
Carbonate dissolution and XRD demonstrated that the positive calcium signal from sand layers on XRF was related to the plagioclase content rather than the marine origin. It undoubtedly helps to avoid misinterpretations of the origin of the element (Additional file 1: Appendix 1). Even though small quantities of carbonate were obtained in some sand layers (< 0.3%) ( Table 2), such values can be related to analytical error instead of real mineral content. Ta has high values in the peaty-mud layers because it tends to be dissolved by weathering and concentrated in the latest phases of fluvial sediment transport, depositing as clay minerals on floodplains (Parker and Fleischer 1968;Fricke and Heilig 2006). Thus, Ca/Ta represents the subtraction of the mud content from the total Ca (i.e., plagioclase), which means that this ratio reflects the relative siliciclastic content along with the sediment core sample and its trend variation. As depicted Page 21 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 by the rows in Fig. 6 and Additional file 3: Appendix 3, Additional file 4: Appendix 4, when the Ca/Ta value decreases, the grain size also decrease ( Figs. 4 and 7). Consequently, the Ca/Ta pattern mimicked the grainsize change. Ta was also helpful for normalizing the elements used in the LDA. As mentioned above, it is closely related to mud (Fig. 6 and Additional file 3: Appendix 3, Additional file 4: Appendix 4, and Additional file 5: Appendix 5), and it can be considered relatively homogeneous and abundant in the core samples. Thus, it serves as a reference for the comparison of the elements related to the sand layers.
On the other hand, CT imagery suggested a rapid landward reduction in sand content and density via a reduction in color intensity in units S8 and S5.

Implications from grain-size and MS data
The grain size analysis was not entirely conclusive for being used as a key factor for tsunami discrimination. For instance, the grain size of the coarsest portion of the sand layers S4, S6, S7, S8, and S9 at H03 was finer than that at H02 (100 m seaward of H03). Horizontal changes in the grain size of the coarsest portion of the sand were not evident between H03 and H04, except for the sand layer S4 (Fig. 7). Thus, landward fining was not convincing among the core samples. This may be related to several factors, such as the paleotsunami sediment source was mostly homogeneous in grain size, the particle size analysis was performed with a small amount of sample, approximately 0.5 g each, the number of core samples is limited, and particularly because the changes in grain size are slight. For instance, when comparing the maximum extent of the sand layer S8,approximately 200 m (Figs. 1 and 3b), to the number of core samples analyzed by grain size, it is not precise to assert that the sand layer is reducing landward grain size. Although the laser diffraction method offers high resolution (Shimadzu 2012), reliability, and reproductivity, the number of samples limits the quality of the results and increases the expected error for interpretation.
On the other hand, regarding sample preparation, sample homogenization using a spatula after organic matter dissolution may not be enough for sample disaggregation, which could leave some particle aggregates that could mislead the measurement. Thus, to improve the precision of the results, it is desirable to increase the sample volume and the number of sample measurements. This reduces systematic errors.
The sand layers, S5, S6, S7, and S8 exhibit slight or no lateral change in the distributions of the grain size. On the other hand, S4, S9, and S10 show a remarkable variation from sand at H02 to mud at H04. Interestingly, although S5 and S8 are the thickest sand layers and exhibit the maximum landward extent, the landward decrease in the mode value is not significant. This can be explained by the homogeneity in the source's grain size or greater flow speed and depth for the distance between H04 and H02, which is only 200 m (Fig. 3). However, some sedimentary characteristics can be observed in the grain-size distributions, which fluctuate along with the sand layer showing the particles' behavior when they settle down. In sand layers, such as S5 in H02 and S8b in H02 from the middle to the upper part, and S8 in H03, in addition to upward fining, grain size skewness tends to rise owing to a lower grain selection (Fig. 7). Some distribution curves exhibit multimodal distribution due to the mixing of sand-mud or sand-granules mixing. Coarser particles exhibiting negative phi values in a secondary mode peak are composed of rounded pumice grains and are distributed mostly towards the base and top of the layer (S5 in H02 and S8 in H02), and further inland (S5 and S8 in H04) (Fig. 7, Table 2). This can be explained by vertical density stratification during sediment transport in the inundation flow and differences in particle settling velocity.
In the mud layer interbedded with S9 in H02, the ripup clast was derived from the mud layer below S9a. The front (Fig. 4c) and lateral (Fig. 4b) views expose differences in axial lengths, which points out the possible transport direction perpendicular to the a-axis; hence, perpendicular to the coastline. In the CT images from cores H02, H03, and H04 (Fig. 4), units S6, S9, and S10 exhibit a remarkable landward decrease in thickness. In addition to the landward decrease in magnetic susceptibility (Fig. 5), landward thinning reflects the progressive reduction in energy and flow capacity during inundation from the sea (Fig. 4).
The lateral variation of MS in S7 is also different, with a maximum in the middle section of the transect (Fig. 5). High magnetic susceptibility values can be associated with high-energy environments (Černý et al. 2016), suggesting a flow velocity peak reached that point is possible. Regardless inverse grading observed in S7 is linked or not to a marine origin, diatom evidence points it out.
The inverse grading appeared in S7 and at the base of S8, both in H02 (Fig. 7). Three different mechanisms can explain this phenomenon. Firstly, Sohn (1997) described the traction-carpet model as a result of a bivariant velocity profile, which results in upward coarsening along with an increase in velocity. Such a model has been developed based on turbidites and hyperconcentrated flow deposits (Sohn 1997), as well as previous paleotsunami studies (Moore et al. 2011;Minoura et al. 2013). Secondly, Yoshii et al. (2017) proposed that the inverse grading at the base of the deposit may be related to processes such as kinetic sieving, geometrical sieving, and spatial differences in Page 22 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 flow speed (Middleton 1970;Kneller and Branney 1995;Dasgupta and Manna 2011), rather than traction-carpet sedimentation. Thirdly, Naruse et al. (2012) assumed velocity changes with time, including acceleration and deceleration cycles, which induce upward coarsening and fining, respectively, and are divided by an internal erosion surface (IES). The IES can be defined as the sedimentary truncation generated by the maximum shear velocity. We interpreted the boundary between S8a and S8b as an IES (Fig. 4). Thus, both S8a and S8b corresponded to the same run-up event and represent a flow cycle of acceleration-deceleration, respectively. The presence of the two subunits in H02 can be explained by the reduction in internal sedimentary units along with landward thinning (Naruse et al. 2012).
All three explanations fit the inverse grading sedimentation style. Although Sohn's model employs a logarithmic profile of flow velocity and its increase with depth reduction and Naruse's model assumes velocity and acceleration changes with time, Yoshii's model includes geometrical properties. One does not necessarily exclude the other. In that sense, during the acceleration phase, a traction-carpet sedimentation style can be developed to generate inverse grading, such as in S8a.
Based on the sedimentological characteristics of the sand layers, it is possible to classify them into three types of sedimentation: (1) normal grading with fast landward sediment run-out, marked by a strong landward grain size change from sand to mud and landward thinning (S4, S9, and S10); (2) massive sands with slight vertical and lateral changes (S6 and S7), and (3) normal grading with thicker and broader sediment deposition (S5 and S8). Such settling can be related to three aspects: first, the sediment-source distribution and availability, which vary seasonally in the shoreface (Dean and Dalrymple 2004); second, the paleosurface formed by successive ponds that cause a rapid run out of the sediment by a rapid decay in the sediment transport capacity of the wave; and third, the surge's energy is controlled by the source magnitude and distance.

Sedimentary environment
GPR imagery displayed the basal and top boundaries that limit the peaty-mud-sand-event sequence, characterized as undulating paleosurfaces (Fig. 14). The freshwater-brackish diatom species that are present along with the core sample and are included mainly in the peat-mud sediment depict a sedimentary environment of wetlands controlled by the influence of freshwater-brackish water and anoxic regimes. These two features indicate that the sedimentary environment corresponds to a marsh complex formed by sequential ponds (Fig. 14).  Such geomorphological setting not only controls the transport and deposition of organic-rich mud sediments, but also the maximum extent of the tsunami sand layers while restricting the progression of inundation.
Bioturbation intensification and appearance/disappearance of diatom species after sand layers S5 and S8 indicate sudden environmental changes (Figs. 4 and 10) (Dixit et al. 1992). Such alterations can be related either to the impact of extreme events on the local ecosystem or the progressive reduction in the water sheet thickness, hence increasing oxygen content and, consequently, bioturbation (Dixit et al. 1992;Nichols 2009).
The presence of freshwater and freshwater-brackish species related to the sand layers, mixed with marine species, can elucidate the reworking of marine provenance sediments with terrestrial material caused by the tsunami wave. Due to the minor thickness of the sand layers, it was not possible to determine whether the position of the marine species was related to a specific stratigraphic level. On the other hand, S8 exhibits that the marine diatoms are in the base and top of the layer. This can be related to the increased mixing of marine and terrestrial materials during the progression and regression of the wave, in which the addition of in situ material is increased. Such behavior is not unusual because the sedimentary background environment is related to a marsh, where entrainment of freshwater species included in terrestrial sediments is expected. In addition, freshwater diatom species have been previously reported in the lithological units surrounding the study area (Kamada et al. 1991).

Recurrence interval between 5600 cal yr BP and 1600 cal yr BP, and regional correlation of tsunami events
The modeled depositional dates for sand layers, with a probability of 95.4% (2σ), are shown in Table 2 and Fig. 12. The model calculation in OxCal was constrained by the To-a event at the top of the sedimentary sequence.
Owing to the minor thickness of S3 and its adjacent peaty-mud layers, it was not possible to constrain its deposition age. Even though T6 was identified as To-Cu, Table 6 Correlation of the sand layers among the core samples

T1
It was defined as the Towada-a event (Table 3), see the discussion in the body text. It continues throughout the core samples (Fig. 3), just below the paddy field soil. It is also consistent that two more tephra events were found underlying it, interbedded with peaty mud, except in H13, H07, and H06, where there was only one, probably due to erosion S1 It only appears in the core sample H02 (the thickest sedimentary record). Therefore, it is difficult to identify a lateral continuity, essentially due to its thickness (Table 1, Fig. 3)

S2
It was traced among H02, H01, and H07 due to its relative stratigraphic position and lithological characteristics S3 It only appears in the core sample H02 (the thickest sedimentary record). Therefore, it is difficult to identify a lateral continuity, essentially due to its thickness S4 It was traceable among H02, H01, H07, H08, and H03 (Fig. 3). Such correlation was made on the grain size, color, and relative stratigraphic position

S5
It was traced throughout all the core samples, representing the sand with the furthest extent. It is the second thickest sand layer, after S8, in the sedimentary sequence. Its position above the peaty-mud layer and T4 (marking horizon), upward fining, and bottom erosive contact allowed its lateral correlation. The dating result obtained from the peaty-mud layer underlying S5 between Ha8 and H04 supported the correlation (see Sect. 5 in the body text)

T4
It is a consistent ash layer, defined as a white clay with homogeneous thickness among the sand layers. Due to its consistency, lateral continuity, and lithological characteristics, it was used as a marker horizon to define the upper boundary of S6 and the lower boundary of S5, considering the peaty mud interbedding between the sand and tephra layers (

S6
It is located below the peaty mud underlying T4; it was found in the core samples H02, H03 and H04. Such correlation was made based on the characteristic lenticular stratification, silt composition, and the thickness of the layers (Table 1) S7 It was traceable among H02, H07, H08, H03, H09, H10, and H04 (Fig. 3). Grain size, the bottom transitional or erosional contact, and the relative stratigraphic position led to its lateral correlation

S8
It was traceable throughout all the core samples between H02 and H06 and represented, in general, the thickest sand layer. Its thickness, upward fining, and relative stratigraphic position allowed its lateral correlation S9 It was traceable among the core samples H02, H07, H08, H03, H10, and H04. The two thin laminas interbedded with the thin peaty mud led to the correlation between H02 and H03. The stratigraphic position allowed the correlation among the other core samples. The dating result obtained from the peaty-mud layer underlying S9 and overlying S10 in Ha8 and H04 supported the correlation (see Sect. 5 in the body text) S10 It was traceable among the core samples H02, H07, H08, H03, H10, and H04. The stratigraphic position, overlying To-Cu and the peaty mud, led to its lateral correlation. The dating result obtained from the peaty-mud layer underlying S9 and overlying S10 in Ha8 and H04 supported the correlation (see Sect. 5 in the body text)

T6
It was defined as the Towada-Cu event (Table 3), see the discussion in the body text. It continues throughout all the core samples at the end of the sequence. In all the cases, it was found as a thick lapilli layer Page 24 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 the lack of insights about its primary or secondary origin limits its usefulness as a boundary in the radiocarbon modeling.
Regarding carbon dating sampling, it was impossible to find fragments of organic matter in the peaty-mud layers for more precise age estimation. Hence, it is expected that the actual depositional ages of the sand layers can be older than the estimated ones due to unknown events such as local erosion and reworking.
The calibrated 14 C results exhibited a mean recurrence interval of 456 yr (320-784 yr, 95% confidence interval) and a coefficient of variation (CoV) of 0.74 (0.57-0.99, 95% confidence interval) for all sand layers (Table 5, Additional file 6: Appendix 6), during the period between 1200 and 5500 yr BP. However, the lack of marine evidence indicates that sand layers S1, S2, and S3 cannot be directly correlated to tsunami events, and their inclusion in the recurrence estimation reduces the reliability  Fig. 13 Carbon dating results and ages correlation among core samples H02, Ha8, and H04. C# notation corresponds to the sample number described in the dating model on Fig. 12 Fig. 14 Paleosedimentary environment interpretation for the base of the peaty-mud sequence. a Ground penetrating radar image obtained around H04 sampling area. Strong reflectors display the top and base of the sand/peaty-mud sequence. b Interpreted paleo environment before the deposition of T1, T2 and T3. It was contructed by using diatom results, sedimentary characteristics, and GPR imagery Page 25 of 29 Velasco-Reyes et al. Progress in Earth and Planetary Science (2022) 9:19 of the calculation. Thus, a mean recurrence interval of 384 yr (320-450 yr, 95% confidence interval) and a CoV of 0.78 (0.68-0.99, 95% confidence interval) was obtained from S4 to S10, from ~ 2700 to ~ 5500 yr BP (Table 5, Additional file 6: Appendix 6). This does not imply that no tsunami affected the Hachinohe area between ~ 1200 and ~ 2700 yr BP. Instead, it suggests no tsunami evidence was recorded or the geological conditions for such time were different and left tsunami evidence without a clear marine imprint. Besides, reworking and tsunami erosion may have removed sedimentological evidence of additional tsunami events that may have occurred among the registered events.
On the coast of Misawa city, Tanigawa et al. (2014a) identified a tsunami deposit formed from 2900 to 4800 cal BP, which can be broadly correlated with events S4, S5, S6, S7, S8, and S9 in the present study (Fig. 15). Hence, making it difficult to correctly correlate the deposit found in Misawa City with any of the sand layers from Hachinohe. To the south, tsunami events from Hachinohe can be correlated with events found by (Takada et al. 2016) in Harashinai (Hirono Town) (Fig. 15). Nonetheless, such chronological and recurrence correlations lack sufficient accuracy owing to the uncertainty in the dating results in Harashinai Site (Takada et al. 2016) and error ranges of the calibrated 14 C ages. Further south, in Noda (Maita), S4 can be correlated to the sand layers "Tse5" and "Layer 4" identified by Ishizawa et al. (2022) and Inoue et al. (2017), respectively (Fig. 15).
Anyway, the northern Sanriku Coast (i.e., Noda and Harashinai) and Hachinohe face both trenches at different angles and at different latitudes. This means that evidence from Hachinohe can be related to the activity of the Kuril and Japan trenches and their junction, while the evidence recorded in Noda and Harashinai can be primarily related to the activity of the Japan Trench (e.g., Figure 7 in Tetsuka et al. (2020)). Such geological records and recurrence intervals reaffirm the importance of Hachinohe as a suitable site for paleotsunami research.
To improve temporal correlations and recurrence estimation, it is imperative to increase the precision of the deposition ages by conducting detailed studies at every site. In addition, numerical modeling of the tsunami sources will improve our understanding of the size of tsunami events that can inundate the tsunami deposit sites examined here (e.g., Figure 7 in Tetsuka et al. 2020).

Conclusions
In this study, we carried out a geological investigation of paleotsunami in Hachinohe, Aomori Prefecture, northern Japan, to reconstruct the regional tsunami history using sedimentary evidence on uplifted marine terraces. Numerical modeling allowed us to determine the best possible survey point and rule out the possibility of storms as a genetic factor. We identified the marine provenance of seven sand layers. Sedimentological analysis indicated that three types of depositional settings existed among the sand layers: (1) normal grading with fast landward sediment run-out, marked by a strong landward grain size change from sand to mud and landward thinning (S4, S9, and S10); (2) massive sands with slight vertical and lateral changes (S6 and S7), and (3) normal grading with thicker and broader sediment deposition (S5 and S8). Such sedimentary differences can be related to sediment-source distribution and availability, the paleosurface formed by successive ponds and vegetation that causes a quick run out of the sediment by a rapid decay in the sediment transport capacity of the wave, and the surge's energy controlled by the source