Evaluating variability in coseismic slips of paleoearthquakes using a flight of displaced terraces across the Kamishiro fault, Itoigawa-Shizuoka Tectonic Line, central Japan

Examining the regularity between events during an earthquake slip leads to an understanding of earthquake recurrence and provides the basis for probabilistic seismic hazard assessment. Generally, scientists use systematic analysis of three-dimensional paleoseismic trenches and numerous offset markers along fault zones to study slip history. Flights of displaced terraces have also been used, under the assumption that the number of earthquakes contributing to the observed cumulative slip is known. This study presents a Monte Carlo-based approach to estimating slip variability from a series of displaced terraces when such an assumption cannot be satisfied. First, we mapped fluvial terraces across the Kamishiro Fault, which is an intra-plate reverse fault in central Japan, and systematically measured the cumulative net slip in the mapped terraces. By combining these measurements with the age of the paleoearthquakes, we estimated the amount of net slip for the penultimate event (PE) and antepenultimate event (APE) to be 1.5 ± 0.2 and 2.7 ± 0.4 m, respectively. The APE slip was twice that of the PE slip and 2.5 times larger than the most recent event, the Nagano-ken-hokubu earthquake, and measured 1.2 ± 0.1 m. This suggests that the APE ruptured along the entire length of the 26 km-long Kamishiro Fault or that there were multiple faults involving adjacent segments. As we are unsure how many earthquakes had occurred since the oldest terrace was formed, we assumed three cases based on available paleoseismic records. In each case, we calculated the slip that could reproduce the cumulative slips within a reasonable range of observed terrace offsets and then estimated the coefficient of variation for coseismic slips (COVs) of paleoearthquakes. The resulting COVs typically fell into the range of 0.3 to 0.5, indicating that, over the last few thousand years, the Kamishiro Fault did not regularly behave as it had done before the 2014 event. Instead, there were large variations in the fault’s coseismic slip, as suggested by the global dataset. Although we acknowledge that our approach may be oversimplified, the Monte Carlo-based approach should help assess the regularity of earthquakes from displaced terraces where limited data are available. 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 twice as large as that of the 2014 event and the PE, 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


Introduction
Understanding the spatiotemporal patterns of past earthquakes provides a guide to estimating the size and timing of future earthquakes. As historical records cover only a short period of time compared with the recurrence interval of large earthquakes (e.g., McCalpin, 2009), geological and geomorphic records contribute greatly when examining whether there is a certain regularity of interevent times and magnitudes during large earthquakes that are caused by a single fault or fault system (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 (e.g., Klinger  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 of the recurrence interval (COVt) is equivalent to the aperiodicity widely used in seismic hazard assessment (e.g., Nishenko and Buland, 1987; Kumamoto and Hamada, 2005). When COVt is 0, it means that the recurrence is completely periodic. As COVt increases, this indicates reducing regularity, and COVt > 1 indicates temporal clustering (See Fig. 14  There are several issues regarding the interpretation of COVs. The first is that the method used to obtain the timing or slip of earthquakes to calculate the COVs needs to be validated. This is because each approach has its own distinct error or uncertainty that would affect the values used to calculate the COVs. For example, when measuring the offset of many geomorphic and geological markers that are distributed along a section of a fault to estimate the average offset of past earthquakes (e.g., Sieh and Jahns, 1984; Zielke et al., 2010Zielke et al., , 2012, the result largely depends on how to measure the offset and assign probability density to each measurement (Zielke et al., 2015). The second issue is that COVs that is derived from points on a fault or a segment might not represent the same thing. When calculating COVs using slips observed at a site on a fault, the results should reflect the characteristics of the fault segment to which the site belongs. In contrast, if COVs is derived from slips averaged over multiple segments, the study area may contain fault segments that have different rupture histories (e.g., Haddon et al., 2016;Kurtz et al., 2018). Therefore, great care must be taken when comparing COVs derived from points and segments. The ability to detect small displacements from stratigraphy or geomorphology also complicates the interpretation of COVs. Here, "small displacements" are those that are hardly preserved because of surface processes after a major earthquake and cannot be detectable as discrete events in the geological record. Flights of displaced terraces, as well as paleoseismic trenching and systematic analysis of offset features, have 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 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 that it can be applied even when the vertical deformation is so great that it requires an unrealistically large trench or an outcrop to identify multiple events (e.g., Bollinger et al., 2014). It can also be applied to areas where excavation of deep trenches is hindered by local conditions, such as high groundwater levels (McCalpin, 2009). However, it is often difficult to ascertain the number of paleoearthquakes that contributed to the cumulative displacement observed on each terrace, which inhibits 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 using a displaced terrace sequence across the Kamishiro Fault. A 9 km-long section to the north of the Kamishiro Fault was ruptured by the Nagano-ken-hokubu earthquake (Mw = 6.2) on November 22, 2014, which caused a vertical displacement of up to 1 m (e.g., Okada et al., 2015;Ishimura et al., 2019). We first mapped the fluvial terraces and measured the cumulative dip slip following a semiautomated 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 from paleoseismic trenches and historical accounts. To reveal the plausible slip variability and the recurrent behavior of the fault, we performed a simple Monte Carlo simulation, taking into account the uncertainty of cumulative slip and age constraints. This approach allowed us to identify the slip regularity of the fault from observable terrace flights whose age constraints were limited and not accurate enough to determine the number of earthquakes since terrace formation.

Study area
The Kamishiro Fault is a 26 km-long, north-northeast trending reverse fault located at the

Methods
We mapped fluvial terraces at Oide using a 1 m meshed digital elevation model (DEM) and aerial photographs taken in the 1940s by the United States Army. As the topography has been significantly altered from its original form by human activity, several terrace risers are barely visible, even with a high-resolution DEM. It was only the historical photographs that allowed us to map the original topography. Therefore, we mapped the terraces based on aerial photography interpretation and georeferenced the resulting images using ArcGIS. This process enabled us to accurately digitize the mapping results and measure the cumulative displacement of each terrace. We used a radiocarbondating technique to determine the age of the terraces. Immediately after the 2014 earthquake, the upslope-facing fault scarp stopped the river flow (  (Fig. 4a). 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. 5a). The second reason was to average the along-strike variation of the surface displacements, as T1 extended 240 m along the fault (Fig. 5a), to consider the effect of alongstrike fluctuation.

Results And Discussion
Previous studies have revealed that the PE resembled the most recent event in terms of slip volume (Katsube et al., 2017) and seismic intensity distribution (Tsuji, 2003); however, the APE slip is still unknown. In this section, we first estimate the amount of PE and APE slip based on the results of dipslip 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 labeled the mapped terraces T0 to T5, from the oldest to the most recent (Fig. 5a). We also mapped terraces Th1 and Th2 along the Himekawa River, which flowed down to the northeast.
However, because these terraces were subparallel to the fault traces and did not cross the fault zone, they were not used. 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 meters, 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 and change in land use indicated in aerial photographs taken in the 1940s, they are now largely unpreserved because of artificial alterations. On the hanging wall, terrace deposits consisting of clast-supported gravel were exposed, whereas, on the footwall, the terrace gravel was covered by alternating units of sand and mud layers This stratigraphic sequence of footwalls suggests that no erosion occurred after terrace abandonment. Instead, a small deposit was locally submerged, trapping fine sediment in close proximity to the descending fault zone. Therefore, we assume that the vertical separation that was observed across the fault zone primarily reflects single or multiple coseismic slips in past surfacerupturing earthquakes. During the 2014 event, a few subsidiary minor surface ruptures appeared to the east of the prominent fault zone. However, aerial photograph interpretation found no discernible scarp along these secondary faults.
We determined the terrace ages of T1, T2, and T5 using Oxcal (Bronk Ramsey, 2009). We chose these three terraces because the other surfaces were not as well preserved on either or both the hanging walls and the footwalls. Radiocarbon samples were used directly from the top of the terrace deposits and the bottom of the swamp deposits, as shown in Fig. 3, to establish the age of the T2 surface. The modeled emergent age of the T2 terrace was in the range of 1,290-1,460 cal BP (1 σ). For the age constraints on the T1 and T5 surfaces, the following approach was taken because of the lack of radiocarbon material. There are several fluvial terraces on the left bank of the Matsukawa River, one of which is approximately 5 m above the T1 surface. The radiocarbon age of a sample obtained from a post-emerged terrace deposit should predate the formation of the T1 surface. We estimated the age of the T1 surface to be 2,710-4,960 yBP (1 σ) using samples taken at Loc. 1 (Fig. 1b) and Loc. 6 ( Fig. 2). Similarly, we determined the age of T5 as 380-1,080 yBP (1 σ) based on the samples taken at Loc. 2 (Fig. 2a) and Loc. 3 (Fig. 3a). Table 1 Radiocarbon ages Loc. (Fig. 1, 2  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. 6). The calculated cumulative dip slip was 7.9 ± 0.5 m for T1, 4.6 ± 0.4 m for T2, and 1.7 ± 0.3 m for T5 (1 σ). The dip slip for each terrace with a 1 σ range of ages is plotted in Fig. 7. Note that dip slip shown in Fig. 7 includes the amount of coseismic slip at the 2014 earthquake derived from differential LiDAR (Ishimura et al., 2019).
On the basis of these dip-slip measurements and terrace ages, we were able to estimate the dip-slip rate in the last 3-5 ky at the Oide site to be 1.5-3.1 mm/yr. Assuming that the fault dip is 54°-62°, the equivalent vertical slip rate is 1.2-2.7 mm/yr, which is consistent with the rates obtained at Iimori (> On the basis of the dip-slip and age constraint for each terrace 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. 7). Simple subtraction allows us to estimate the dip slip of the PE to be 1.7 ± 0.3 m and that of the APE to be 2.9 ± 0.5 m. Since the dip slip of the most recent event at Oide was 1.2 ± 0.1 m (Ishimura et al., 2019), the slip of the PE was a little larger than that of the 2014 event. Similar conclusions are given by Katsube et al. (2017), that is, the PE showed a slightly smaller displacement than the 2014 slip on the paleoseismic trench walls in 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 iso-seismicity (Tsuji, 2003), support the assumption that the events are similar. In contrast, based on his compilation of the paleoseismic trenches along the northern and central ISTL, Okumura (2001) argued that 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 Okumura (2001).

Inter-event variability in coseismic slip at Oide
To explore the variability of each 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 emergence of terrace T1. According to Okumura (2001), one to three surface-rupturing earthquakes occurred during this time period. Therefore, we can assume that one, two, or three earthquakes occurred during that time. Figure 8 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 MCSST (Fig. 6). Next, we calculated the dip slip of the recent three events ( in Fig. 8: represents the dip slip of the ith most recent event). We then calculated the dip slip of the simulated events (S 4 , S 5 , S 6 ), which occurred between the APE and the emergence of terrace T1, such that We will discuss the minimum (S 1 ) and maximum (S 3 ) threshold in the next paragraph. We repeated the calculations until we obtained 10,000 sets of realizations, which satisfied the necessary conditions  (Table 2). Table 2 COVs with and without an upper limit on dip slip amount (1 σ There are slight differences between the three cases; however, each case shows a similar slip variation (Fig. 9). Our results show that the COVs generally falls between 0.3 and 0.5 (1 σ range; Our results, showing that COVt is generally greater than COVs, agree with those that have been reported in recent studies (e.g., Nicol et al., 2016;Zielke, 2018). This discrepancy may reflect the natural characteristics of earthquake recurrence but can be, in part, attributed to the varying abilities of geological and geomorphic markers to record an earthquake with small displacements or short recurrence intervals. Immediately after an earthquake, surface processes begin to degrade the surface ruptures. This degradation makes it difficult for the earthquake displacement to be preserved at the surface. Errors in topographical data, such as those recorded by a DEM or optical satellite images, may be too large to resolve small displacements from the terrain (e.g., Zielke et al., 2015).
These difficulties often prevent an observer from identifying small displacements, resulting in an apparent decrease in COVs. In addition, written records, archaeological evidence, and dating strategies (e.g., Lienkaemper and Bronk Ramsey, 2009) can distinguish temporally clustered earthquakes (e.g., Scharer et al., 2017;Wechsler et al., 2018), which tends to increase COVt. These characteristics contribute to the apparent discrepancy between slip regularity and recurrence interval.

Discrepancy between observed COVs and true COVs
We consider that the recent three paleoseismic records, MRE, PE, and APE are complete. However, we may have overlooked minor events before the APE that relied entirely on paleoseismic evidence at one of Okumura's sites. If this is the case, the modeled COVs may not represent the true variability in slips from event to event (e.g., Hecker et al., 2013). Here, we attempt to quantify the effect of such missing events on the discrepancy between observed COVs and true COVs by introducing additional small events to the simulation. The dip slips associated with those small events () are modeled as follows: where S min denotes the minimum amount of dip slip. Technically, to estimate the true variability in slips at the Oide site, S min must be the minimum surface displacement that can occur at the site, which is difficult to determine. Instead, we set S min to be 0.25 and 0.50 m based on a threshold of detection developed by Hecker et al. (2013). The threshold of detection is based on expert comments and represents the threshold displacement with a 50% chance of being detected as a discrete event in a geological record. The threshold values of 0.25 and 0.50 m are generally applied to paleoseismic trenches, with lower thresholds considered to be better stratified (see Figure B1 in Hecker et al., 2013 for the details). Except for these conditions, the calculation process was almost identical to the process of ignoring small events in the previous section. We repeated this calculation until we generated 10,000 realizations. For example, when n = 4 and m=1 in Eq. 6 ( Fig. 11), we first calculated true COVs from all events (three recent events and two simulated events). As for the observed COVs, we added the slip of Eq. 4* to that of Eq. 4 and then calculated the COVs. This is because, even when a small event is not identified as a discrete event, its slip is conserved and interpreted in conjunction with other events (Fig. 11).
The results show that COVs increases as more small events are added to the dataset ( Fig. 12 and Table 3). It should be noted that, although unlikely, COVs can also decrease (negative value in Fig. 12), which is the opposite of what would be expected (e.g., Zielke et al., 2015). This happens when the average slip of an observed event is close to the slip of a missing event. Given that the recurrence intervals of the recent three events are approximately 300 and 900 years (Fig. 7), true COVs is expected to be around 0.5-0.7 (Fig. S2). These results demonstrate the significant impact of ignoring small earthquakes when estimating COVs and reaffirm that the Kamishiro Fault exhibits a more irregular recurrent behavior than expected before the 2014 event (e.g., Okada et al., 2015).

Conclusions
Accurate cumulative slip measurements and detailed chronologies are indispensable components 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 twice as large as that of the 2014 event and the PE, 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 between 0.31 and 0.49, which is comparable with the ranges previously reported (e.g., Hecker et al., 2013;Zielke, 2018). We demonstrated that it is possible 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.    Schematic illustration to explain our strategy for determining terrace age. The emergent age for each terrace is bracketed by post-event 14C age from the swamp deposit and the pre-event terrace-forming deposit. Step 2 as a plausible dip slip and its uncertainty. Results are shown in Fig. 6.    Flowchart illustrating the procedure to calculate the COV of the slip. S_2014 is the dip slip of Figure 11 An example of the observed and true paleoseismic record. Eq4* is an event that cannot be identified as a discrete event. Therefore, the slip of Eq4* is interpreted together with the slip of Eq4.