Topographic Evolution Involving Co-Seismic Landslide, Deformation, Long-Term Folding and Isostatic Rebound: A Case Study on the 2004 Chuetsu Earthquake

Co-seismic landslide volume information is critical to understanding the role of strong earthquakes in topographic and geological evolution. The availability of both preand post-earthquake high-resolution digital elevation models (DEMs) provides us with the opportunity to develop a new approach to obtain robust landslide volume information. Here, we propose a method for landslide volume estimation and test it in the Chuetsu region, where a Mw 6.6 earthquake occurred in 2004. First, we align the DEMs by reconstructing the horizontal difference. Then, we quantitatively obtain the landslide volume in the epicentral area by differencing the preand post-earthquake DEMs. We convert the landslide volume into the distribution of average catchment-scale denudation and the resulting long-term crustal rebound. Our findings reveal that the Chuetsu earthquake mainly roughens the topography in the low-elevation Chuetsu region. Our results indicate that the preserved topography not only is due to the uplift caused by fault-related folding on the hanging wall of the Muikamachi fault but also undergoes erosion caused by seismically induced landslides and crustal rebound also modifies the topography in the long term. This study confirms that the differential DEM method is a valuable approach for quantitative analysis of topographic and geological evolution.


Introduction
The role of strong earthquakes and earthquake induced landslides is increasingly recognized as critical to understanding topographic and geological evolution. The volume of the co-seismic landslides contains important information about the landscape evolution. Previous studies use area-volume scaling relationships to calculate landslide volumes (e.g., [1][2][3][4]) but the results may contain large uncertainties [5]. Different co-seismic landslide volume results have been reported for the 2008 Wenchuan earthquake [1,[4][5][6][7]. Recently, high-resolution and multi-temporal light detection and ranging (LiDAR) digital elevation models (DEMs) or DEMs generated from stereo pairs of remote sensing images have been proven valuable in monitoring geomorphic and co-seismic surface deformations [5,[8][9][10][11][12][13]. Differencing of digital elevation models (DEMs) from different epochs has been widely used to study changes in surface topography duo to various natural (e.g., [14][15][16][17]) and anthropogenic (e.g., [18]) causes. The data source used to build DEMs are easy to acquire by satellite imaging or aerial scanning of the epicentre area immediately after the earthquake. Researchers even do not have to worry about road blockage by the co-seismic landslides. The high precision of the differential DEM method also made it a widely used method in many cases of earthquake-induced landslides. At the Longmen Shan region, seismic induced landslide volume was estimated by differencing pre-and post-Wenchuan Earthquake DEMs [5]. The landslides triggered by the 2008 Iwate-Miyagi Nairiku Earthquake were detected by using RGB color composite image of SRTM and ALOS/PALSAR InSAR DEMs [19]. Based on DEMs built by stereo SPOT6/7 data, the 2015 Gorkha (Nepal) earthquake-induced landslide volume was estimated [20]. By differences in pre-and post-earthquake DEMs, we could quantitatively analyse topographic changes and evaluate landslide volumes.
Thrust earthquakes are the main cause of surface uplift, and hence mountain building. Strong thrust earthquakes and landslides have been realized to play important roles in topographic and geological evolution in regions of steep slope and high elevation, such as the marginal zones of high plateaus of the Himalaya [21][22][23][24], Taiwan [25] Longmen Shan [4,5,25] and Andes [26]. A recent study in the foothills of Peru suggests earthquake is the primary trigger of landslides [26]. Another study [27] on a strong earthquake at the margin of the Tibet plateau (Nepal) also suggests a strong relation between earthquakes and landslides. Previous studies have demonstrated that landslides limit the slope [28,29] and the height of mountain peaks above adjacent river valleys in steep orogenic regions [22,26,30]. Quantifying the erosion rate is critical to understanding the role of tectonic events in mountain building. However, due to longterm mountain-building and topographic and geological evolution, previous studies have mainly been regional studies based on sparse thermochronological dating [31,32], cosmogenic dating [33][34][35] or modern hydrological observations [36] in regions of high mountains. These regional studies cannot show the details of how tectonic events affect the geological and topographic evolution of low mountain regions. Strong historic earthquakes are the most recent tectonic events that provide us with a valuable opportunity to study the role of such events in the current topographic evolution of low mountain regions. However, in regions with low mountains, the role of strong earthquakes in geological and topographic evolution is rarely reported. As our approach to obtain the landslide volumes is by differencing the pre-and post-earthquake DEMs, it is independent of factors such as local relief and hill slope, therefore it is also applicable to low mountains.
The 2004 Mw 6.6 Chuetsu earthquake occurred in Niigata Prefecture in Japan, where the local relief of the epicentral area is low, with a maximum elevation of 765 m ( Figure 1). In this study, we use high-resolution pre-and post-earthquake DEMs [37] to study the topographic changes due to the Chuetsu earthquake by comparing the slope angle, slope aspect, relief and roughness pre-and post-earthquake. The co-seismic denudation distribution pattern is also analysed using the co-seismic landslide volume with the availability of multi-temporal high-resolution topographic data. We can further analyse the long-term influence of the isostatic response to the landslide unloading on the surface height. We finally discuss the role of earthquakes in the geological and topographic evolution of the Chuetsu area with low mountains.

Tectonic Setting
The 2004 Mw 6.6 Chuetsu earthquake occurred at Chuetsu, Niigata Prefecture, Japan, where the convergent plate boundary between the Amur and Okhotsk plates is located ( Figure 1 [40][41][42]). The epicentral area has low elevation with a maximum of 765 m and is composed of Holocene to Miocene sedimentary and volcanic rocks. The sediments were mainly formed in the early Miocene, concurrently with the opening of the Japan Sea ( Figure 1). The sediments have been folded under an E-W to WNW-ESE compressional stress field since~2-3 Ma ago (e.g., [43,44]). The continued compression has deformed the strata and landforms and caused repeated seismicity in the Chuetsu area. The Shinano River is the main river flowing through the Chuetsu area, and the flood plain is mainly composed of Holocene to late Pleistocene sediments ( Figure 2). The mountain area is mainly composed of Pleistocene to Pliocene sediments, accompanied by a tectonic window revealing Jurassic sedimentary rocks ( Figure 2). The fault track constructed from fold geometries is consistent with the fault geometries revealed by aftershocks of the 2004 earthquake [40] and the inclined antithetic shear with a dip of 85 • in the hanging wall above a single reverse fault indicates it is a fault-related fold [40]. Geological studies indicate that the surface folding was developed since~3.5 Ma ago, following the extension stage of the Japan Sea [45]. The uplift of the mountain is proposed to be due to faultrelated folding caused by thrusting along the N-S-trending Muikamachi-Bonchi-Seien fault ( Figure 2 [40,46,47]). The 2004 Chuetsu earthquake was a thrust-dominated earthquake with minor lateral motion [48]. A co-seismic surface rupture zone 1 km in length has been reported, with 20 cm of vertical co-seismic offset and a lateral offset of less than 20 cm on a previously unmapped fault [48], which lies along the northward extension of the Muikamachi fault ( Figure 2 [38,39]). Hence, the most probable causative fault of the Chuetsu earthquake is the Muikamachi fault, according to the focal mechanism and locations of surface ruptures [48]. Previous studies have mapped the subsurface fault with a detachment depth of~10-13 km, which agrees well with the distribution of aftershocks [40,47,48]. These studies have found that the fault-related folding on the hanging wall of the Muikamachi fault was responsible for the growth of the geological structures [40,[47][48][49]. Palaeo-seismology studies reveal at least two strong earthquakes that occurred in the past 9000 years prior to the occurrence of the 2004 Chuetsu earthquake [48]. The co-seismic displacements of the two palaeoearthquakes were almost identical at~1.5 m and almost 15 times that of the 2004 event (~10 cm). The 2004 Chuetsu earthquake triggered thousands of co-seismic landslides, which dramatically modified the local topography [50][51][52][53]. Hence, the geological and topographic evolution in the epicentral area should be closely related to the co-seismic landslides caused by strong earthquakes.

Data
The pre-earthquake DEM is 10 m in resolution with absolute vertical precision within 2.5 m. The 10-m-resolution DEM is generated from stereo pairs of aerial photographs or topographic maps that cover the whole area of Japan and are published by the Geospatial Information Authority (GSI) of Japan (freely available at http://fgd.gsi.go.jp/download (accessed on 15 December 2020)). The post-earthquake DEM has a 2-m resolution with a root mean square (RMS) error within 0.12 m and is generated from airborne LiDAR data surveyed in 2005 with a point density larger than 1 pt/m 2 , released by the GSI of Japan [37]. These DEMs are of higher precision than those used in our previous studies in the Wenchuan area [5] (Figure 3). The landslide inventory map is interpreted based on high-resolution aerial photographs by the National Research Institute for Earth Science and Disaster Prevention (NIED), Japan (Figure 4 [50][51][52][53]). The geological information is derived from the 1:200,000 geological maps provided by the Geological Survey of Japan (GSJ) (Figures 2 and 4a, freely available at https://gbank.gsj.jp/seamless/index_en.html (accessed on 15 December 2020)) and the active fault map of Japan ( Figure 1 [38,39]).

Differential DEM
With the availability of high-resolution pre-and post-earthquake DEM data, especially the LiDAR DEM, the differential DEM method is widely used in detecting topographic changes [5,9,54], co-seismic deformations [8,10,55] and sediment budgets [12,13]. Previous studies have shown that the differential DEM method using multiple-scale and multiple-source DEMs is effective in detecting topographic changes caused by co-seismic landslides [5,6,54]. The available pre-and post-earthquake DEMs in the Chuetsu area provide us with the opportunity to study the topographic changes caused by co-seismic landslides due to the 2004 Chuetsu earthquake. In the differential DEM method, the precise geo-reference and correlation between the multi-temporal DEMs is one of the key issues before subtracting [56]. To analyse the topographic changes under a compressional environment, we are mainly interested in the elevation changes. As the DEMs are data of matrix, when conducting correlating and differencing between the pre-and post-seismic DEMs, the matrix dimensions must agree with each other. The airborne LiDAR DEM was down-sampled to 10 m to match the pre-earthquake DEM, under the rule that the volume in each 10 × 10 m area stays unchanged before and after the down-sampling. Example of down-sampling work can also be found in the research of near-field deformation of the El Mayor-Cucapah Earthquake , in which study high-resolution (≥9 returns per square meter) post-earthquake LiDAR data and lower-resolution (5-m-pixel) LiDAR data were used. At first, the horizontal differences between the pre-and post-earthquake DEMs were calculated and then reconstructed by back-slipping the horizontal differences. The coregistration of optically sensed images and correlation (COSI-Corr) software was developed to measure sub-pixel ground deformation using optical satellite and aerial images and has an accuracy of 1/10 of the input pixel size [10,[57][58][59]. In this study, we estimated the horizontal differences between the pre-and post-earthquake DEMs using the COSI-Corr software (freely available at www.tectonics.caltech.edu/slip_history/spot_coseis/index. html (accessed on 15 December 2020)), following Zhou et al.'s method [10]. We used a correlation window of 64 pixels followed by 32 pixels with a step of 4 pixels (40 m). The subpixel matching procedure was performed on the frequency content, which is more accurate than the statistical correlator [10,57]. Consequently, we obtained the NS (Figure 3a) and EW (Figure 3b) components of the horizontal differences between the pre-and post-earthquake DEMs ( Figure 3). Then, by reconstructing the mean horizontal differences in both directions across the whole DEM, we obtained precisely geo-referenced and correlated DEMs. Finally, by differencing the pre-and post-earthquake DEMs, we obtained the elevation changes caused by the Chuetsu earthquake (Figures 3c and 4b). In order to better avoid canopy effect of the forest area, we derive the landslide volumes using the elevation change data within each landslide area, meanwhile, the data within the densely forested region is not used (Figure 4). The canopy effect is partially avoidable because most of the landslides did not occur in the forested region, because the seismic landslides mainly occurred on steep slopes which are favourable for landslides. If there are strong canopy effect, then the positive and negative values of the DEM difference should include the canopy, which should be apparently different, hence, could not be comparable. However, our results also show the summing results of the negative and positive value is almost zero. Meanwhile, the largest landslide clearly showed the source and deposit areas, which occurred in the low mountains composed of late Miocene to Pliocene non-marine sediments (Figure 4a). Meanwhile, the derived landslide volumes also showed consistent results that deep-seated landslides are the main contributor to the landslide volume ( Figure 5). The above results indicate that the canopy effects do not seriously affect our results of the landslide volume and the distribution pattern of denudation. An earthquake-triggered landslide inventory map (Figure 4a) is made by the Geographical Survey Institute of Japan, based on aerial photos, photographed at October 24 after the earthquake, and field survey within 2 days after the earthquake [53], in scale of 1:10,000 and 1:25,000. The elevation changes (Figure 4b) revealed by the differential DEM show a similar distribution with the earthquake-triggered landslides on the inventory map ( Figure 4a). Hence, it can be inferred the elevation changes are mainly due to co-seismic landslides. Additionally, the differential DEMs method only detects elevation changes during the time period between the generations of the pre-and post-earthquake DEMs, usually less than 10 years. Therefore, we suppose the landslides independent of the earthquake is small in number comparing with the seismic triggered ones. As a result, elevation changes ( Figure 4b) are inferred to be mainly due to earthquake-induced landslides.

Topographic Analyses
Steep slopes are prone to landslides [29,60,61], such as the co-seismic landslides triggered by the 2008 Wenchuan Mw 7.9 earthquake, which mainly occurred on slopes with angles larger than 30 • [62]. Previous studies have widely applied slope angle, slope aspect, relief and roughness as the four main topographic features for geomorphological research, which can be used to analyse the topographic changes due to co-seismic landslides. Statistical comparison of pre-and post-earthquake topographic features has been proven useful in analysing co-seismic topographic changes [5]. In this study, based on the downsampled 10-m-resolution pre-and post-earthquake DEMs, we compare the pre-and postearthquake slope angle, slope aspect, relief and roughness ( Figure 6). The co-seismic displacement is reported to be less than 20 cm in the epicentral area [48]; The results obtained by differencing the pre-and post-DEMs show up to 29 m of vertical displacements (Figure 4b), much larger than the reported maximum co-seismic displacements [48]. Thus, in this study, we statistically compare the topographic changes within each landslide polygon ( Figure 6).

Catchment-Scale Denudation Depth
Landslides are the dominant mass wasting process in humid uplands [63]; thus, quantifying the denudation using landslide volumes is reasonable. The topographic changes within a landslide area should be much larger than the~10 cm co-seismic displacements. Thus, deriving landslide volume using the subtracted DEM is reliable. We calculated the co-seismic landslide volumes in the epicentral area by summing the elevation changes within the landslide area and multiplying the summed value by the area of one pixel (100 m 2 ). There are positive values and negative values for each single landslide because there are source and deposit areas, respectively (Figure 4b,c). Hence, we first summed both the positive and negative values for each landslide, then summed the absolute value together and finally took the average values as the landslide volumes for these landslides. The catchments were derived using the flow accumulation data of the streams and outlet data for each catchment from the pre-earthquake DEM with stream lengths as short as 5 km (Figure 7). To obtain the distribution pattern of the catchment-scale denudation depth, we summed the landslide volumes within each catchment. The average denudation depth was obtained by dividing the summed landslide volume by the catchment area (Figures 7 and 8).

Isostatic Response
We analysed the isostatic response to the removal of the landslide mass produced by the 2004 Chuetsu earthquake. As the sum of the negative and positive value of the DEM difference is nearly zero, the denudation volume is almost equal to the volume deposited. The Shinano River had not removed a significant part of the landslide mass until the LiDAR scan. Although we do not precisely know the duration for the mass removal, it may take much shorter time comparing with a seismic cycle [25]. Isotatic rebound will compensate the mass unloading in the long term, hence, further modify the topography [5,64]. Here, we implement the flexural isostatic rebound (FIR) model by [65].
By this model, flexural rigidity of the lithosphere distributes the isotatic response, leading to long-wavelength, low-amplitude tilting of the earth's surface [65,66]. We smoothed the co-seismic landslide volumes to 500 × 500 m pixels and divide the total pixel volume by an area of 0.25 km 2 . Then, we got the denudation depth distribution in 500-m pixels (Figure 8a). The pixel averaged denudation depth pattern has a higher resolution than the catchment-scale depth when calculating the isostatic response. Another input is the effective elastic thickness (T e ) of the lithosphere in the epicentral region. Using topography and observed gravity anomalies, Takeshi et al. [67] estimated the effective elastic thickness across the Japanese island arc lithosphere. For the region where Chuetsu locates, the T e is estimated to be 7 km. Other inputs: Young's modulus = 10 GPa; Poisson's ratio = 0.25; density of crust = 2700 kg/m 3 ; density of mantle = 3300 kg/m 3 [65,66]. By Equation (1), the flexural rigidity is~3.05 × 10 21 N m.
E is Young's modulus, T e is the effective elastic thickness, and v is Poisson's ration. Finally, we have the predicted isotatic rebound in response to the landslide mass removal (Figure 8b). Since the FIR model do not predict crustal rebounding rate, we assume the isostatic rebound is a long-term process comparing with the time needed for landslide unloading.

Results
Using the COSI-Corr software, we obtained a NS component difference of −0.27 m with a standard deviation (STD) of 0.47 m and an EW component difference of 0.21 m with a STD of 0.51 m (Figure 3b,c). Before differencing the DEMs, we first reconstructed the corresponding differences in the NS and EW components. Then, we obtained the elevation changes by subtracting the pre-earthquake DEM from the post-earthquake DEM. Flat ground surfaces located far from the epicentral area should be stable during an earthquake; i.e., the real elevation changes should be nearly zero. Hence, the obtained elevation differences in such regions represent the accuracy of the differential DEM in this study (Figure 3c). The elevation differences have a very small range at the nondeformed flat region (Figure 3c). We then obtained a mean elevation change for the whole region of −0.24 m with a STD of 0.26 m (Figure 3c) by eliminating the elevation changes between at the flat region, i.e., setting the delta-z values within this range to 0. In the landslide region, the elevation changes can reach tens of metres, which are much greater than at the flat region; hence, analysing the topographic changes using the preand post-earthquake DEMs is reliable. Therefore, we directly derive the landslide volume without using the centimetre scale (within the error level) co-seismic deformation data. As shown in Figure 4c, the scarp and toe of the largest landslide are clearly shown on the differential DEM map. By the denudation depth distribution revealed by Figure 8a, the denudation volume in each 500-m pixel can be calculated. Then we summed up the denudation volumes across the landslide region and estimated the total landslide volume as~22.5 million m 3 . Using high-accurate DEM from SPOT5 satellite images, Tsutsui et al. [68] reported a~25 million m 3 total landslide volume induced by the 2004 Chuetsu earthquake. These two estimated total landslide volumes are comparable with each other. The catchment-scale average denudation depth distribution shows a maximum catchment scale average denudation of 419 mm, which is not located immediately above the surface ruptures (Figures 7 and 8 [48]). The maximum denudation correlates with the uplift pattern in the epicentral area suggested by the fault-related folding on the hanging wall of the Muikamachi fault (Figure 9 [40,46,47]). The co-seismic topographic changes in the epicentral area show consistent increases in slope angle, relief and roughness ( Figure 6). The comparison of the pre-and postearthquake topographic features within each catchment also indicates that the average hillslope, relief and roughness increased after the Chuetsu earthquake ( Figure 6). The slope aspect decreased from 0 • -135 • and 270 • -360 • and increased from 135 • -270 • (Figure 6). The observed slope aspect changes might be related to the co-seismic lateral displacement of the 2004 Chuetsu earthquake, which was reported as~10-20 cm [48].
Although the duration for the landslide removal is unknown, the isotatic rebound will compensate the mass unloading and further modify the topography in the long term. Using the flexural isostatic rebound calculation by [65], we predict the crustal rebound distribution in the epicentral area (Figure 8b). The predicted rebound maximum, 290 mm, also correlates with the uplift pattern (Figure 9 [40,46,47]), lying~8 km from the surface ruptures (Figures 8 and 9).

Effectiveness of the Differential DEM
The occurrence of the Chuetsu earthquake provides a valuable opportunity to test the effectiveness of the differential DEM in detecting the co-seismic topographic changes, denudation caused by co-seismic landslides and the isostatic response to the mass removal, hence revealing the tectonic process and topographic growth in regions with low mountains. This study might be the first to investigate topographic changes using the differential DEM method in regions with low mountains. The key problem in differential DEM study occurs when the pre-and post-earthquake DEMs are from different sources, which might represent different surfaces such as the topographic or bare-earth surface. Usually, a DEM from stereo pairs of images should represent the topographic surface, including the forest canopy, and the LiDAR-derived DEM should represent a bare-earth DEM. However, in this study, no systematic elevation errors are observed. Meanwhile, the blank area in Figure 4c shows the elevation differences are small, hence canopy effect will not seriously affect the calculation result of the landslide volume. The differential DEM results should represent the canopy of the forest, which indicates that the pre-and post-earthquake DEMs refer to the same surface; i.e., our results represent the real co-seismic elevation differences. However, due to the precision and resolution of the DEMs, they cannot clearly show co-seismic elevation changes less than~20 cm. In this study, based on the pre-and post-earthquake DEMs, we statistically analyse the topographic changes caused by the Chuetsu earthquake in terms of the slope angle (Figure 7), slope aspect (Figure 7), relief (Figure 7), roughness (Figure 7), catchment-scale denudation (Figures 7 and 9) and we also predicted the isostatic response to the landslide removal in the long term (Figures 8 and 9). The slope angle, relief and roughness are all coseismically increased after the Chuetsu earthquake on the landslide scale ( Figure 6) and catchment scale (Figure 9) in the epicentral area ( Figure 6). The slope aspect changes show decreases from 0 • -135 • and 270 • -360 • and increases from 135 • -270 • , which might be associated with the co-seismic lateral displacement along the N-S-trending rupture [48]. The comparisons of pre-and post-earthquake data suggest that the Chuetsu earthquake mainly increased the topographic relief ( Figure 6), which is consistent with the role of long-term seismic landslides [22,26,28,30].The statistical analyses prove that the differential DEMs method is an effective approach to detecting the co-seismic topographic changes in regions with low mountains.

Surface Deformation
The earthquake on a thrust fault modifies the topography in multiple ways. The fault related folding increase the topography on the hanging wall. The earthquake induced landslides roughen the landscape. Erosion reduces the relief and will eventually remove the landslides from the epicentre area. The unloading of the landslides will lead to isostatic rebound. The present topography is the net result of these topo-modifying factors. The epicentre region of the 2004 Chuetsu earthquake provides us an excellent opportunity to study how these factors modify the topography.
Using the differential DEM method, we obtained the denudation depth of the landslide. Then the denudation depth was converted into landslide volumes. We also predicted the long-term elastic crustal rebound caused by the removal of the landslides. We used a swath profile (Figures 7 and 8) normal to the fault to extract the topographic relief, denudation depth and predicted isostatic rebound patterns (Figure 9a-c) to find their relations with the faulting event. Formerly reported post-earthquake deformation pattern [69] (Figure 9d) and fold deformation pattern [40,46,47] (Figure 9e) are also projected on Figure 9.
The modern relief is interpreted as the net result of fault related faulting, denudation and isostaic rebound. However, we found that in the Chuetsu region, the denudation depth distribution (Figure 9b), the predicted isotatic rebound (Figure 9c) and the fold deformation pattern (Figure 9e) do not show a direct correlation with the topographic profile (Figure 9a). The denudation and isostatic rebound patterns show a correlation with the deformation pattern (Figure 9e) suggested by fault-related folding on the hanging wall of the Muikamachi Fault (Figure 9e, [40]). The fault-related folding suggests that the main uplift area lies~8-10 km from the Muikamachi fault on the hanging wall [40,46,47,49]. Duo to the folding, the deformation pattern also shows subsidence near the surface rupture zone (Figure 9e), a similar pattern with the post-seismic displacement [70] profile (Figure 9d) extracted from the interferogram generated by synthetic aperture radar (SAR) data [69]. The subsidence zone corresponds to the maximum topographic relief (Figure 9a). Instead, in the steep Wenchuan area, previous studies have found that long-term high denudations are concentrated in a narrow zone along the Longmen Shan thrust belt, revealed by erosion rates from kyr-scale cosmogenic 10 Be and Myr-scale low-temperature thermochronology dating methods [31,34,35]. The highest co-seismic denudation is also mainly concentrated in the narrow corridor between the co-seismic surface ruptures produced by the 2008 Wenchuan earthquake [5].
Two possibilities might explain why the topographic relief profile (Figure 9a) does not directly correlate with the uplifting pattern in the Chuetsu region ( Figure 9). The first possible reason might be the erosional differences. High erosion occurred in the high-uplift area due to fault-related folding, but the area near the fault experienced low erosion and uplift ( Figure 9). Thus, almost homogeneous topographic relief was preserved in the epicentral area under the coupled processes of co-seismic uplift and denudation. The correlation between the denudation, isostatic rebounds and uplift also indicates that the topographic growth in the orogenic belt is closely related to the deformation associated with major faults, especially in regions where uplift is dominated by fault-related folding [40,49]. Hence, strong earthquakes should have been playing an important role in mountain growth in the epicentral area since the initiation of compression and folding at~2-3 Ma ago [43,44]. The second possible reason is that the high topographic relief near the Muikamachi fault might be due to strong palaeo-earthquakes associated with large co-seismic vertical offsets, as revealed by trenching [71]. Palaeo-seismological studies reveal that at least two strong earthquakes occurred during the past~9000 years, which produced co-seismic vertical displacements almost 15 times greater than that produced by the Chuetsu earthquake [48,71]. A vertical slip rate of 0.4 mm/yr was calculated using the dating results of the displaced units in a trench [71]. The Chuetsu area is composed of young sediments mainly formed since the Miocene, which have been under continuous thrusting, folding and isostatic response to wast unloading since~2-3 Ma ago [44,48]. The total uplift can be estimated as~800-1200 m, which is consistent with the current mountain height of~700-800 m considering that erosion also occurs caused by landslides.

Conclusions
By quantitatively comparing pre-and post-earthquake high-resolution DEMs, we detected denudation depth of the earthquake induced landslides, calculated the total volume of co-seismic landslide and predicted long-term isostatic rebound in response to the landslide unloading. Our statistical analyses reveal that the slope angle, relief and roughness ( Figure 6) are increased after the 2004 Chuetsu earthquake in the epicentre area. The analysis results indicate the differential DEM method is a robust approach to detecting the topographic changes induced by earthquake in low mountainous regions.
We report the role that the Mw 6.6 Chuetsu earthquake played in the geological and topographic evolution of a young and low mountain region. The preserved mountain peaks not only are uplifted by thrusting and folding but also undergo erosion caused by seismically induced landslides. In response to the following landslide removal, the epicentral area will undergo a long-term isostatic rebound, further modifying the topography. The denudation depth, predicted isostitic rebound and fold deformation distributions (Figure 9b,c,e) yielded by the swath profile normal to the fault show a different pattern with the topographic relief profile (Figure 9a). The main deformations lay~4-8 km from the fault track. That pattern is different from changes observed in the old and steep Longmen Shan orogenic region, with the co-seismic denudation of the 2008 Wenchuan earthquake concentrated at the surface rupture zones. we suggest that strong earthquakes might play different roles in geological and topographic evolution in gentle and steep mountain regions.
Author Contributions: J.L. and Z.R. contributed to conceiving the idea, prepared the figures, data compilation, and wrote the manuscript text. T.O. provides the LiDAR data. Z.R., T.O., P.Z. and S.U. discussed the contents of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the Institute of Geology, China Earthquake Administration.