Evaluating variability in coseismic slips of paleoearthquakes from an incomplete slip history: an example from displaced terrace flights across the Kamishiro fault, central Japan

Examining the regularity in slip over seismic cycles leads to an understanding of earthquake recurrence and provides the basis for probabilistic seismic hazard assessment. Systematic analysis of three-dimensional paleoseismic trenches and analysis of offset markers along faults reveal slip history. Flights of displaced terraces have also been used to study slips of paleoearthquakes when the number of earthquakes contributing to the observed displacement of a terrace is known. This study presents a Monte Carlo-based approach to estimating slip variability using displaced terraces when a detailed paleoseismic record is not available. First, we mapped fluvial terraces across the Kamishiro fault, which is an intra-plate reverse fault in central Japan, and systematically measured the cumulative dip slip of the mapped terraces. By combining these measurements with the age of the paleoearthquakes, we estimated the amount of dip slip for the penultimate event (PE) and antepenultimate event (APE) to be 1.6 and 3.4 m, respectively. The APE slip was nearly three times larger than the most recent event of 2014 (Mw 6.2): 1.2 m. This suggests that the rupture length of the APE was much longer than that of the 2014 event and the entire Kamishiro fault ruptured with adjacent faults during the APE. Thereafter, we performed the Monte Carlo simulations to explore the possible range of the coefficient of variation for slip per event (COVs). The simulation considered all the possible rupture histories in terms of the number of events and their slip amounts. The resulting COVs typically ranged between 0.3 and 0.54, indicating a large variation in the slip per event of the Kamishiro fault during the last few thousand years. To test the accuracy of our approach, we performed the same simulation to a fault whose slip per event was well constrained. The result showed that the error in the COVs estimate was less than 0.15 in 86% of realizations, which was comparable to the uncertainty in COVs derived from a paleoseismic trenching. Based on the accuracy test, we conclude that the Monte Carlo-based approach should help assess the regularity of earthquakes using an incomplete paleoseismic record.


Introduction
Understanding the spatiotemporal patterns of past earthquakes helps to estimate the size and timing of future earthquakes (e.g., Zielke et al., 2015). As historical records cover only a short time period compared with the recurrence interval of large earthquakes (In this article, large earthquakes mean those with surface displacement.) (e.g., McCalpin, 2009), geological records significantly contribute to examining whether there is a certain regularity of recurrence intervals and magnitudes of large earthquakes (e.g., Shimazaki and Nakata, 1980;Schwartz and Coppersmith, 1984;Grant, 1996;Weldon et al., 2004;Zielke et al., 2015). Although many studies have presented evidence that faults appear to behave regularly with regard to the size and recurrence intervals of large earthquakes (e.g., Klinger et al., 2011;Berryman et al., 2012), others have reported fluctuations (e.g., Chen et al., 2007;Schlagenhauf et al., 2011;Rockwell et al., 2015;Komori et al., 2017;Scharer et al., 2017;Mechernich et al., 2018;Wechsler et al., 2018). This raises the question of how variable earthquake recurrences are in terms of their size and recurrence interval.
The coefficient of variation (COV; a ratio of standard deviation to mean) is a statistical index of variability that is useful when discussing the recurrent behavior of large earthquakes. The COV for coseismic slips (COVs) represents slip variability over seismic cycles. When COVs is 0, it implies the repetition of the same amount of slip, and as COVs increases, the inter-event difference in slip amount also increases (see Fig. 14 in Zielke et al., 2015 for a graphical representation). One of the advantages of calculating COVs is that it facilitates the discussion about which earthquake magnitude-frequency distribution best describes the observation (e.g., Hecker et al., 2013;Zielke 2018), which helps to develop a realistic seismic hazard assessment (e.g., Hecker et al., 2013;Field et al., 2014;Nicol et al., 2016).
One of the problems in estimating COVs is that at least three or more slips of paleoearthquakes are necessary to calculate COVs, which is difficult to achieve. Paleoseismic trenching is the most popular way to reveal a slip history of a fault; however, it requires an ideal site where past earthquake displacements are preserved by well-stratified layers that are rich in datable materials. Ascertaining the number of earthquakes that contributed to the observed displacement is often quite difficult. This occurs owning to the poor stratigraphy (e.g., a massive unit) or limits of geologic and geomorphic records preserving past earthquake slips (e.g., Zielke et al., 2015), which is also applicable to other methods to reconstruct a slip history, such as the systematic analysis of offset features (e.g., Sieh and Jahns, 1984;Zielke et al., 2010) or the analysis of the 36 Cl concentration of a limestone fault scarp (Zreda and Noller, 1998;Schlagenhauf et al., 2011;Mechernich et al., 2018).
Flights of displaced terraces have also been used to reconstruct the timing and displacement of paleoearthquakes (e.g., McCalpin et al., 2009;Bollinger et al., 2014;Berryman et al., 2018). The basic assumption is that the difference in cumulative displacement between successive terraces resulted from displacement of a single earthquake (McCalpin et al., 2009). Only when this assumption holds does the displaced terrace sequence reveal an accurate slip history. This approach is useful in cases where the vertical displacement is too large to identify multiple events from a paleoseismic trench or an outcrop (e.g., Bollinger et al., 2014). Another advantage is that it can be applied to areas where excavation of deep trenches is hindered by local environmental conditions, such as high groundwater levels (McCalpin, 2009). However, similar to other approaches mentioned in the previous paragraph, it is often difficult to ascertain the number of paleoearthquakes that contributed to the cumulative displacements observed on each terrace, which prevents us from reconstructing the history of large earthquakes from a succession of displaced terraces.
This study presents a method for exploring inter-event variability in coseismic slips at a site especially when a detailed paleoseismic record is not available. As a case study, we used a displaced terrace sequence across the Kamishiro fault in central Japan and tested the accuracy of our approach by performing the same analysis for a fault whose slip history is well constrained. We first mapped fluvial terraces and measured the cumulative dip slip following a semi-automated method developed by Wolfe et al. (2020). We then estimated the coseismic slip of the penultimate event (PE) and the antepenultimate event (APE) of the Kamishiro fault based on the ages of the paleoearthquakes derived from paleoseismic trenches and historical documents. To reveal plausible slip variability of the fault, we performed a simple Monte Carlo simulation that considered the uncertainty of cumulative slip and age constraints. This approach allowed us to identify the slip regularity of the fault from terrace flights whose age constraints were limited and not accurate enough to determine the number of earthquakes that the terraces experienced since their formation.

Study area
The Kamishiro fault is a 26 km long, north to northeast trending reverse fault located at the northernmost part of the Itoigawa-Shizuoka Tectonic Line active fault system (ISTL, Shimokawa et al., 1995;Okumura, 2001, Fig. 1). The northern ISTL including the Kamishiro fault is one of the three major segments of the ISTL and is composed of east dipping reverse faults. Based on the relocated aftershocks of the 2014 Nagano-ken-hokubu earthquake with a three-dimensional velocity structure, the subsurface portion of the Kamishiro fault dips 30-45°SE at a depth of 0-4 km and 50-65°SE at a depth of greater than 4 km (Panayotopoulos et al., 2016). The Kamishiro fault merges with the Otari-Nakayama fault at approximately 4 km in depth (black lines in Fig. 1b), which initially formed as a normal fault that occurred with the opening of the Sea of Japan during the Miocene age and reactivated as a reverse fault (e.g., Okada and Ikeda, 2012;Panayotopoulos et al., 2016). The regional stress regime changed from an east-west extension to an east-west compression in the late Pliocene age (Sato, 1994;Sato et al., 2004;Ikeda et al., 2004), and continuous GPS observation indicates that the contraction strain axis is N 110°E (Sagiya et al., 2004). As the steeply dipping Otari-Nakayama fault is unfavorably oriented to the current east-west compression, it has been inactive since the early Pleistocene age (Kato et al., 1989;Ueki, 2008), and instead the Kamishiro fault formed as a footwall shortcut thrust (Panayotopoulos et al., 2016).
Timings of paleoseismic events are constrained by both paleoseismic trenches and historical documents. Katsube et al. (2017a) excavated a trench at Iimori (Fig. 1b) and found that the PE of the Kamishiro fault occurred after C.E. 1645 (2 σ) and produced vertical displacement comparable to that of the 2014 event (~50 cm). They also argued that the instrumental seismic intensity of the 2014 Fig. 1 Active fault map around the ISTL active fault system. a Active faults in central Japan and the ISTL composing of three major segments, each separated by change in slip direction and surface geometry of the fault trace (Nakata and Imaizumi, 2002). b The 2014 surface rupture, denoted as blue lines, occurred on approximately one-third of the previously known Kamishiro fault. The trench sites correspond with those in Fig. 10. The Otari-Nakayama fault is from Kato et al. (1989) and Nakano et al. (2002). ISTL Itoigawa-Shizuoka Tectonic Line, JMA Japan Meteorological Agency earthquake was similar to the estimated intensity of the 1714 Otari earthquake (Fig. 2, Tsuji, 2003) and correlated the 1714 earthquake of the PE on the Kamishiro fault. Toda et al. (2016) conducted trenching at Oide and Iida (Fig. 1b) and estimated the timing of PE was after 1695 and 1659 (2 σ), respectively, which is consistent with the result of Katsube et al. (2017a). The APE of the Kamishiro fault is thought to be a historical earthquake either in C.E. 762 or 841 (e.g., Katsube et al., 2017a). At the Oide site on the Kamishiro fault, Toda et al. (2016) reported a paleoseismic event dated between C.E. 420-865 (2 σ). Paleoseismic events of similar ages (around C.E. 700-900) were identified at multiple sites between Kamishiro and the Gofuku river ( Fig. 1) (e.g., Okumura, 2001;Maruyama et al., 2010). This overlap of the event ages led Maruyama et al. (2010) to conclude that the northern and a part of the central ISTL ruptured together during that period. Based on historical documents (Usami et al., 2013), Maruyama et al. (2010 correlated the event around C.E. 700-900 with a historical earthquake either in C.E. 762 or 841. Although several studies suggest additional events might have occurred at Oide between the ninth and the eighteenth century (e.g., Sugito et al., 2015;Toda et al., 2016), clear evidence of such events have not been found.
In this study, we focused on the Oide site, located in the northern part of the 2014 rupture zone, 2 km west of the epicenter (Figs. 1b and 3). There are two major rivers in this area: the Matsukawa River and the Himekawa River (Fig. 3a). The Matsukawa River flows eastward and formed a massive fan and flights of terraces. The Kamishiro fault runs through the middle of the fan, creating west-facing scarps, which dams the Matsukawa River to form a swamp along the foot of the scarp (Fig.  3c). During the 2014 earthquake, surface ruptures appeared mainly along the pre-existing fault scarp and were accompanied by minor secondary faulting, such as flexure deformation and rupture on branch faults (Fig. 3a: Okada et al., 2015;Ishimura et al., 2019). The total amount of dip slip at Oide was 1.2 ± 0.1 m, which corresponds with the maximum value over the entire rupture area (Ishimura et al., 2019). We mapped fluvial terraces at Oide using a 1 m meshed digital elevation model (DEM) taken after the 2014 earthquake and aerial photographs taken in the 1940s by the US Army. As the topography has been significantly altered from its original form by human activity, several terrace risers are hardly recognizable even with the highresolution modern DEM. It was only the historical photographs that allowed us to map the original terrace boundaries. Therefore, we mapped the terraces based on aerial photography interpretation and georeferenced the results using ArcGIS. This process enabled us to accurately digitize the mapping results and measure the cumulative displacement of each terrace. We used a radiocarbon-dating technique to determine the age of the terraces. Immediately after the 2014 earthquake, the upslope-facing fault scarp stopped the river flow, which locally submerged and formed a swamp on the footwall side (Figs. 3 and 4). Thus, radiocarbon ages obtained from the upper part of the terrace deposits should predate the terrace formation age, and the ages from the bottom of the swamp deposits overlying the terrace deposits on the downthrown side should postdate the terrace formation age (Fig. 4). In practice, because it is difficult to confirm if the sample is taken from the bottom of swamp deposits or swamp deposits over colluvium at the foot of a scarp (Fig. 4a), a younger limit on a terrace age tends to be overestimated. We collected several radiocarbon samples from the swamp deposits, and together with the radiocarbon ages reported in previous studies (Sugito et al., 2015;Toda et al., 2016), we estimated the ages of terrace formation using Oxcal v.4.4 (Bronk Ramsey, 2008) and IntCal20 (Reimer et al., 2020). Slip measurements often include errors due to the observer's lack of experience or knowledge. To minimize the effects of human error, we systematically measured the cumulative dip slip of a terrace using a semi-automated tool called the Monte Carlo Slip Statistics Toolkit (MCSST; Wolfe et al., 2020). Measuring fault offset requires the management of various uncertainties that can affect the measurement and even lead to Fig. 4 Typical stratigraphy along an uphill-facing scarp. a Schematic illustration to explain our strategy for determining terrace age. The age for each terrace is bracketed by a 14 C age from the swamp deposit and that from terrace-forming deposit. b An exposure of paleoseismic trench by Toda et al. (2016) at Oide. Both Loc. 4 and Loc. 5 are located within this trench. c A drilled core obtained at Loc. 4. Swamp deposit continues from the surface to the depth of 2.4 m and covers a gravel layer that corresponds to that shown in b misinterpretation (e.g., McGill and Sieh, 1991;Klinger et al., 2011;Zielke et al., 2012Zielke et al., , 2015Zielke, 2018). The MCSST helps us find a plausible dip slip by iteratively calculating the slip using several key fault parameters (e.g., extent and average slope of hanging/footwall and fault dip) and their uncertainties (Thompson et al., 2002). According to Thompson et al. (2002), the dip slip (S) is given by where x p is the position at which the fault intersects the scarp, δ is the fault dip, and m h, f, s and b h, f, s are the slope and intercept, respectively, of the linear regression lines for the hanging wall, footwall, and scarp surface. The calculation process is summarized in Fig. 5. First, we defined the extent of the footwall, hanging wall, and fault scarp, based on topography and the 2014 surface displacement (Ishimura et al., 2019) (Fig. 5a). As an example, a significant 85 cm vertical separation of terrace T1 appeared at x = 320 m, and a minor 25 cm vertical separation occurred at x = 450 m (Fig. 6b). Therefore, we considered an area of x > 450 m to be the hanging wall. Then, we calculated the slope and interception of linear regression lines for the profile of each component. The probability density function of a slope and an intercept follows a normal distribution with a standard deviation equal to the standard error of linear regression. We then assumed the fault dip based on the ratio of the vertical and horizontal components of the 2014 displacement (Ishimura et al., 2019) (Fig. 5a). The PDF of the fault dip is given by a uniform distribution. In addition, we created 30 transects on each terrace and calculated the cumulative dip slip on each transect (Fig. S1). There were two reasons for doing this. The first reason was that the measurement depends on the relative orientation between the topographic transect and the fault strike. This caveat primarily relates to the measurements for terraces T2 and T5 (Fig. 6a). The second reason was to average the along-strike Step 1: For each of the 30 profiles on each terrace (Fig. S1), define extents of footwall, hanging wall, and fault scarp on each topographic profile. The fault dip is determined from a ratio of the vertical and horizontal displacement of the 2014 earthquake. The PDF of the fault dip is assumed to be uniform. Slope and interception of topographic profiles (footwall, hanging wall, and fault scarp) are determined from linear regression, and their probability densities are assumed to be normally distributed with mean and standard deviation being the best fit and standard error of the regression. The fault is assumed to intersect with the scarp within the lower 5% of the scarp, and its PDF is uniform. b Step 2: Calculate the dip slip using the PDF defined in step 1. We performed 10,000 realizations for each of the 30 profiles and calculated the mean (denoted as S in the figure). c Step 3: Calculate the mean and standard deviation of S in step 2 as a plausible dip slip and its uncertainty. Results are shown in Fig. 9 Takahashi and Toda Progress in Earth and Planetary Science (2021) 8:15 variation of the surface displacements, as T1 extended 240 m along the fault (Fig. 6a).

Results and discussion
Previous studies have revealed that the PE resembled the most recent event in terms of the amount of slip (Katsube et al., 2017a) and seismic intensity distribution (Tsuji, 2003); however, the APE slip is still unknown at Oide. In this section, we first estimate the amount of PE and APE slip based on the results of dip-slip measurements and paleoseismic history. Then, we perform the Monte Carlo simulation to estimate unknown pre-historic events and assess the slip regularity of the Kamishiro fault from several past surface-rupturing earthquakes.

Terrace mapping and age determination
We mapped terraces around Oide based on the relative height of each terrace tread and labeled terraces as T0 to T5, from the oldest to the most recent (Fig. 6a). Considering the meander wavelength of the Matsukawa River (> 1 km), original forms of the terrace risers should have been almost straight around the fault (approximately several hundred meters from the fault). Therefore, we correlated terraces on the hanging and footwall side of the fault based on the geometry of terrace risers. We also mapped terraces Th1 and Th2 along the Himekawa River, which flowed down to the northeast (Fig. 6a). However, because these terraces were subparallel to the fault traces and did not cross the fault, they were not used. Our result is generally consistent with previous studies that also used aerial photographs taken in 1940s for terrace mapping (Togo et al., 1996;Matsuta et al., 2006;Sugito et al., 2015), and the difference is that we classified lower terraces (T3 to T5, Fig. 6a) in more detail. The relative heights of all the mapped terraces from the current Matsukawa River bed were lower than 10 m, and the vertical separations between successive terraces were from 1 to 5 m ( Fig. 6b-d), which may indicate that these terraces were formed over a short period of time.
Although the terrace risers are clearly marked by the paddy field boundaries indicated in aerial photographs taken in the 1940s (Fig. S2), they are now hardly preserved because of artificial alterations. Thus, some terrace risers are not visible in the DEM (Fig. 6a, S2). On the hanging wall, terrace deposits consisting of clastsupported gravel were exposed, whereas, on the footwall, the terrace gravel was covered by alternating units of sand and silt layers that were deposited after the abandonment of the terraces (Sugito et al., 2015;Toda et al., 2016). This stratigraphic sequence of footwalls suggests that erosion occurred after terrace abandonment was minimal. Instead, a downthrown side of the fault scarp was locally submerged, trapping fine sediments along the foot of the scarp (Fig. 3c). Therefore, we assume that the vertical separation that was observed across the fault primarily reflects single or multiple coseismic slips in past surface-rupturing earthquakes. During the 2014 event, secondary surface ruptures appeared to the east of the prominent fault scarp (Fig. 3a). However, there was no discernible scarp along these secondary faults in old aerial photographs taken in 1940s. We determined the terrace ages of T1, T2, and T5 using Oxcal (Bronk Ramsey, 2009). Ages of T3 and T4 could not be determined because we lacked radiocarbon samples to constrain their ages. The model output is summarized in Fig. 7 and the original output is in Fig.  S3. Two radiocarbon samples used to establish the age of the T2 surface were taken from the top of the terrace deposits and the bottom of the swamp deposits over the terrace deposits (Fig. 4, Table 1). As explained in the "Methods" section, the former sample should yield the age before and the latter after the abandonment of T2 (Fig. 4). The modeled age of the T2 terrace was in the range of 1150-1480 cal BP (2 σ) (Fig. 7). For ages of T1 and T5, because we lack ages from terrace deposits, we used ages of swamp deposits over terraces older than T1 and T5 as older limits. There are several fluvial terraces on the left bank of the Matsukawa River and one of them (T L ) is at approximately 5 m above the T1 surface suggesting terrace T L is older than T1. Same as terraces at the Oide site, the Kamishiro fault created a westfacing scarp on T L and swamp had developed at the foot of the scarp soon after the abandonment of T L , which is evidenced by a humic layer covering a fluvial gravel layer (Loc. 1, Fig. S4, Sugito et al., 2015). Thus, the radiocarbon age from the humic layer over T L should predate the formation of the T1 surface. We estimated the age of the T1 surface to be 2590-5010 yBP (2 σ) using samples from Loc. 1 and Loc. 6 ( Fig. 7, S4). Similarly, we used ages of swamp deposit over T3 (or T2) (Loc. 3, Fig.  S4) and over T5 (Loc. 2, Fig. 8) to determine the age of T5 as 200-1110 yBP (2 σ) (Fig. 7, Table 1).

Cumulative dip slip of mapped terraces
We modeled the cumulative dip slip of the T1, T2, and T5 surfaces iteratively using the MCSST (Wolfe et al., 2020). We estimated the fault dip on each profile based on the ratio of vertical to horizontal (fault-normal) component of the surface displacement of the 2014 earthquake (Ishimura et al., 2019). The estimated dip angles were 54-62°for T1, 48-56°for T2, and 44-52°for T5. For convenience, we used the LiDAR DEM taken before the 2014 event to calculate the cumulative dip slip. Using the measured slips from 30 profiles for each offset terrace (Fig. S1), we produced composite histograms and determined the mean and standard deviation of the cumulative dip slip (Fig. 9). The calculated cumulative dip slip was 7.9 ± 0.5 m for T1, 5.0 ± 0.5 m for T2, and 1.6 ± 0.6 m for T5 (1 σ). The dip slip for each terrace with an error (1 σ range) is plotted in Fig. 10. The dip slip shown in Fig. 10 is the sum of the dip slip mentioned above and the amount of coseismic slip at the 2014 earthquake derived from differential LiDAR (Ishimura et al., 2019).
Based on these dip-slip measurements and terrace ages, we were able to estimate the dip-slip rate of the Kamishiro fault. Dividing cumulative dip slip of T1 (7.4-8.4 m) by its age yields, the average dip-slip rate of 1.5-3.1 mm/year in the last 2.7-5 ky at Oide. Assuming that the fault dip is 54-62°, the equivalent vertical slip rate is 1.2-2.7 mm/ year, which is consistent with the rates obtained at Loc. 1 (> 1.6 mm/year) (Fig. 1, Sugito et al., 2015).

Dip slip of the PE and the APE
Identifying the amount of slip per event from a terrace sequence requires a complete catalog of paleoearthquakes to ensure that the cumulative slip difference between successive terraces was produced by a single earthquake. If we assume the PE occurred in C.E. 1714 and the APE occurred around C.E. 700-900 (e.g., Maruyama et al., 2010;Toda et al., 2016;Katsube et al., 2017a), it follows T2 experienced the PE and the APE, and T5 experienced the PE. However, as suggested by several studies (e.g., Sugito et al., 2015;Toda et al., 2016), there might have been additional events between the ninth and the eighteenth century. If there were any missing events, the cumulative dip slip of T3 should be between T2 and T5. Thus, to examine the possibility of such missing events, we also measured the cumulative dip slip of the T3 using the MCSST. The fault dip angle was assumed to be 44-51°, and the topographic profiles  used were shown in Figure S1. The calculated dip slip of T3 was 1.8 ± 0.7 m (1 σ) and was almost equal to that of T5 (1.6 ± 0.6 m). This similarity indicates that there was no surface-rupturing earthquakes between the formation of T3 and T5 or that surface displacement that occurred in that period was small. The similarity also indicates formation or abandonment of terraces at Oide can be associated with lateral migration of a channel and is not necessarily linked to faulting. Moreover, the post-2014 InSAR analysis did not detect any significant slip on the surface rupture (Omata et al., 2017). Therefore, it should be reasonable to assume that the cumulative slip of T2 resulted from coseismic displacements of the C.E. 1714 and the 700-900 earthquakes the slip of T5 from the displacement of the 1714 earthquake. Based on the calculated dip slips, terrace ages, and the dates of the paleoearthquakes, we were able to calculate the dip slip at the PE and the APE. Before the 2014 earthquake, the T5 terrace only experienced the PE, whereas the T2 surface was deformed by the PE and the APE (Fig. 10). The dip slip of the T5 terrace represents that of the PE, 1.6 ± 0.6 m. The difference between the dip slip of T2 and T5 should be equivalent to that of the APE, 3.4 ± 0.8 m. Since the dip slip of the 2014 event at Oide was 1.2 ± 0.1 m (Ishimura et al., 2019), the slip of the PE was slightly larger than that of the 2014 event.
Our result is consistent with that of Katsube et al.  Fig. 3c. To determine the depth of the gravel layer of T5, we inserted a geoslicer and a soil sampler until they hit gravel and stopped; we repeated this measurement many times. We confirmed that the thickness of swamp deposits was 30-50 cm from the surface, which suggested that radiocarbon samples we obtained were from the bottom of the swamp deposits and, thus, were suitable for estimating the age of T5 (2017a) who found that the slip of the PE was similar to that of the 2014 event at the center of the 2014 surface rupture zone (Iimori site in Fig. 1b). Two independent results showing the similarity of coseismic slips at both the PE and the 2014 earthquake, supported by the PE's similar seismic intensity (Fig. 2, Tsuji, 2003), indicate that the events are similar. In contrast, according to the compilation of the paleoseismic trenches along the northern and central ISTL (Okumura, 2001;Maruyama et al., 2010), the event in the year 841 or 762 ruptured the entire northern segments of the ISTL and a part of the central ISTL. Our results show that the APE produced a much larger surface displacement than the two recent events, which is in agreement with previous studies (Okumura, 2001;Maruyama et al., 2010).

Inter-event variability in the coseismic slip at Oide
To explore the variability of slip per event at Oide, we performed a simple Monte Carlo simulation to determine additional events and their possible slip amounts during the time period between the APE and the abandonment of terrace T1. According to the result of a paleoseismic trenching at Kamishiro 5 km apart from Oide (Okumura, 2001), one to three surface-rupturing earthquakes occurred on the Kamishiro fault during the simulated time period. Therefore, we considered three cases where one, two, or three earthquakes ruptured both Oide and Kamishiro during that time. We considered that all the cases assume (1) no more than three earthquakes ruptured Oide and (2) at least one earthquake ruptured both Oide and Kamishiro, which will be discussed in the next paragraph. Figure 11 shows the Monte Carlo simulation procedure. First, we modeled the dip slips of the 2014 event and the terraces by sampling them from a normal distribution. The mean and standard deviation of probability distributions used in this step are based on the results produced by the MCSS T (Fig. 9). Next, we calculated the dip slip of the recent three events (S 1 , S 2 , S 3 in Fig. 11: S i represents the dip slip of the ith most recent event at the Oide site). We then calculated the dip slip of the simulated events (S 4 , S 5 , S 6 ), which occurred between the APE and the abandonment of terrace T1, such that We will discuss the minimum (S min ) and maximum (S max ) threshold later. Because we assume at least one paleoearthquake out of three identified at Kamishiro (Okumura, 2001) ruptured Oide, any of the simulated events at Oide can be correlated to one of the second, third, and fourth most recent events at the Kamishiro site (Fig. 10). For example, in case 3, when only the sixth youngest event at Oide ruptured Kamishiro, S 6 can be associated with any events except for the first one (occurred C.E.~800) identified at Kamishiro (Fig. 10, Okumura, 2001). We repeated the calculations until we obtained 10,000 sets of realizations, which satisfied the necessary conditions for the slip, and then calculated the COVs to assess the plausible variability of each slip per event on the Kamishiro fault.
Here we justify two assumptions that we made in the simulation: (1) no more than three earthquakes ruptured Oide and (2) at least one earthquake ruptured both Oide and Kamishiro. The first assumption is based on the smaller dip slip rate at Oide than at Kamishiro. Our result revealed that the dip slip rate at Oide in the last 2.7-5 ky was 1.5-3.1 mm/year whereas the rate at Kamishiro since 10 ka was 2.8-4.1 mm/year (Niwa et al., 2018). Given this difference in slip rates, it should be reasonable to assume that the number of earthquakes that ruptured at the Oide site during our modeled time period is equal to or less than that ruptured at the Kamishiro site. Regarding the second assumption, Kondo et al. (2019) identified an earthquake dated to~4400 yBP at the southern end of the Kamishiro fault. Based on an empirical relationship between surface rupture length and maximum surface displacement (Matsuda et al., 1980), they concluded that this event ruptured the entire Kamishiro fault and an adjacent fault to the south. Furthermore, although a 1 km-wide step and 10°deflection of surface fault trace between Oide and Kamishiro played a role in stopping the rupture of the 2014 event (Katsube et al., 2017b), an area with these minor geometrical complexities can often be traversed by earthquakes of Mw > 6 ( Wesnousky, 2016, 2017). These observations support the second assumption of our current model.
To eliminate unrealistic estimates, we defined the 2014 slip as the minimum and the APE slip as the maximum limits on the amount of slip based on available paleoseismic records. The dip slip of the 2014 event was chosen as the minimum value because the paleoearthquakes in the modeled time window should have been larger than the dip slip in the 2014 event. This is because Okumura (2001) identified events based on angular unconformity, which is probably insensitive to small displacements near the surface. The fact that Okumura's trench sites experienced a broad uplift (~20-30 cm) without any apparent surface break (Ishimura et al., 2019) also justifies our choice of the minimum value. The maximum  (Toda et al., 2016), Iimori (Katsube et al., 2017a), and Kamishiro (Okumura, 2001)  limit is based on a compilation of paleoearthquake records at 42 sites over the entire ISTL (Maruyama et al., 2010). The results show that an earthquake that took place either in the year 841 or 762 was the greatest earthquake on the Kamishiro fault in the last 12 ky, which supports our choice of maximum threshold. We also performed the same Monte Carlo simulation without the maximum or minimum threshold on the dip-slip amount to see the effect of slip limits on the results. The resultant COVs are generally in the range between 0.30 and 0.54 and decrease with an increase in the number of simulated events when using upper and lower limits on slip per event (Table 2, Fig. 12). Cases without the maximum threshold yielded the similar results. This occurs probably because the maximum threshold (3.4 ± 0.8 m) is close to the total amount of slip allocated to simulated events (2.9 ± 0.7 m). Moreover, the removal of the minimum threshold considerably changed the results. This result emphasizes the importance of properly choosing a minimum threshold of slip in estimating COVs. Notably, the influence of a maximum/minimum threshold can vary depending on the number of simulated events and their total slip. Nevertheless, establishing an upper or lower bound of slip per event is often difficult owning to the lack of paleoseismic records covering a sufficiently long period compared to an average recurrence interval. While an upper bound can be constrained from earthquake scaling relationships among maximum slip and fault dimensions (e.g., Wells and Coppersmith, 1994;Stirling et al., 2013), a lower bound should depend on the way the resulting COVs will be used. When one uses COVs for seismic hazard assessment of a fault that is somewhat distant from populated areas where an earthquake of relatively small magnitude (e.g., Mw < 6) does not cause  severe damage, slips associated with such insignificant earthquakes can be ignored in estimating COVs. When one wants to compare a resulting COVs with those of other faults, a lower bound of slip may be constrained from local erosion rate or how COVs of other faults are derived. The minimum amount of slip that can be geologically or geomorphologically reconstructed depends on the local erosion rate and ability of the method resolving paleo-slips (e.g., Hecker et al., 2013;Zielke, 2018).

Testing the accuracy of the Monte Carlo approach
To test the accuracy of our approach, we performed the same simulation using a well-documented paleoseismic record of the Velino-Magunora fault (VMF) in central Italy (Schlagenhauf et al., 2011). The VMF is a normal fault that is composed of 4 major segments; each one is 10 km long, and the total length is approximately 45 km. Large earthquakes exhumed limestone fault plane along VMF, and an in situ 36 Cl profile of these fault planes records the paleoseismic activity of VMF (e.g., Schlagenhauf et al., 2011;Benedetti et al., 2013). Schlagenhauf et al. (2010Schlagenhauf et al. ( , 2011 have reconstructed the timings and slips of nine earthquakes by measuring 36 Cl concentration at five sites on VMF. These researchers suggested that some of these events might have ruptured one or two segments of VMF or entire VMF, which is similar to the rupture pattern of the Kamishiro fault where the fault ruptured on its own or together with adjacent faults (e.g., Katsube et al., 2017a). Given the quality of the slip history and similarity of the rupture pattern, the slip history of VMF is suitable to examine the predictive power of our method.
We focus on the recent five events identified at site MA3 (see Schlagenhauf et al., 2011 for details) because these events are identified at multiple sites; thus, their slip amounts are expected to be most reliable. Table 3 summarizes the best estimates and uncertainties of slip per event shown in Fig. 5 of Schlagenhauf et al. (2011). Analogous to the Kamishiro fault (step 3, case 2 in Fig. 11), slips of the fourth and fifth events are assumed to be unknown and computed. The simulation steps are as follows.
(1) Sample slip amounts of five events from pre-defined PDFs and calculate actual COVs, which will be compared Fig. 12 Probability distributions of COVs. a Case 1 assumes that one earthquake occurred at Oide during the period between T1 formation and the APE. b Case 2 assumes that two earthquakes occurred at Oide in the modeled period. c Case 3 assumes that three earthquakes occurred at Oide in the modeled period. d A composite histogram of all cases. All cases include the maximum and minimum thresholds of slip per event Takahashi and Toda Progress in Earth and Planetary Science (2021) 8:15 to simulated COVs.
(2) Compute slips of the fourth and fifth youngest events. These are randomly chosen between the maximum and minimum thresholds instead of PDFs. The sum of two simulated slips equals to that of actual events defined in the first step.
(3) Calculate simulated COVs using slips of three most recent events (from step 1) and two simulated events (from step 2), and compute the difference between the actual and simulated COVs. (4) Repeat steps (1) to (3) 10,000 times. We assigned triangular PDFs to each slip per event. The peaks of these PDFs occur at best estimates and upper/lower bounds correspond to max/min uncertainties. When sampling slips from PDFs, for simplicity, we assume that slip per event is independent of each other; although in reality, each slip varies depending on those of other events (Schlagenhauf et al., 2010). The maximum threshold of slip per event is 3.3 m, which is determined by the scaling relationship between surface rupture length (45 km) and maximum displacement (Wells and Coppersmith, 1994). The minimum threshold is 0.5 m, which is the detection limit of the 36 Cl approach (Schlagenhauf et al., 2011). The absolute difference between actual and simulated COVs is shown in Table 4. The numbers in the table represent a fraction of all realizations whose absolute difference is less than the specified value. Thus, 70% of realizations yielded an absolute difference of < 0.10, and the absolute difference of 86% of realizations < 0.15. It is very unlikely that the difference exceeds 0.20. These errors are almost half of the uncertainty of COVs estimates for the Kamishiro fault, where the one sigma uncertainty range is~0.2 (Table 2). This discrepancy is due to the large uncertainties in the terrace offsets we used in the simulation (± three-sigma range). If we consider only the most probable value (± one sigma range) for the terrace offsets, the uncertainty in COVs of the Kamishiro fault decreases to ± 0.05-0.08, which is consistent with the result of the accuracy test.
It is difficult to determine whether this accuracy is sufficiently good or not because the uncertainty of estimated COVs is not often presented (e.g., Zielke et al., 2018). However, COVs reported by Wechsler et al. (2018) is in the range between 0.50 and 0.65. Their estimate was based on the correlation of buried channels across a strike-slip fault, which is one of the most popular and straightforward ways to study a slip history of a strike-slip fault. The error in COVs of our accuracy test is comparable to the uncertainty reported by Wechsler et al. (2018); thus, we believe that our current approach should help explore a possible range of COVs when a detailed paleoseismic record is not available.

Conclusions
Accurate cumulative slip measurements and detailed chronologies are indispensable when revealing slip history from a flight of displaced terraces. We measured the cumulative dip slip based on key fault parameters and their uncertainties to obtain plausible estimates. We used these estimates, together with paleoseismic records, to reconstruct the amount of dip slip in the PE and APE. The results showed that the slip of the APE was three times as large as that of the 2014 event, emphasizing the need to assess the earthquake regularity of the Kamishiro fault. We also estimated the possible range of COVs at the Oide site based on the available paleoseismic records and the Monte Carlo approach. The resulting COVs is most likely to be~0.42.
We demonstrated that our Monte Carlo approach helps to estimate a possible range of COVs, even when it is uncertain how many earthquakes contributed to the observed cumulative displacement of a terrace sequence. Identifying multiple paleo-slips requires an ideal environment for past surface displacements to be preserved and detected as discrete events. Finding a site that fully satisfies such conditions is not always easy, which limits the number of detectable events at one site. We recognize that our modeling approach may be too simplistic to explain some fundamental issues in interpreting earthquake regularity. However, we believe that the concept itself should facilitate research into earthquake variability when using paleoseismic records that are missing some events or those that are conflicting.
Additional file 1: Figure S1. Profiles used to measure cumulative dip slip of the mapped terraces. Since the profiles do not intersect with the