A future scenario earthquake and ground motion hazards for Kathmandu, Nepal

Since Nepal is an earthquake-prone country due to the collision of the Indian and Eurasian plates, it is crucial for its capital of Kathmandu to evaluate ground motion hazards from a future large earthquake. For this purpose, we constructed a realistic scenario earthquake with realistic rupture parameters in a likely location. To obtain the location, a new distribution of the rupture zones of large historical earthquakes along the Main Himalayan Thrust (MHT) was obtained by integrating the distribution from a previous study and the results of trench surveys. Global Navigation Satellite System observations indicated that the plate boundary is strongly coupled from the southern boundary of the MHT to a depth of approximately 10 km and there is almost no lateral change in the coupling. This implies that all regions along the MHT have similar rates of strain increase. Therefore, it is most probable that the rupture zone of the oldest previous event will rupture as a scenario earthquake. In the new distribution, the 1255 earthquake is the oldest. However, large earthquakes occurred in 1934 and 2015 within its rupture zone. Thus, we adopted the area obtained by removing the 1934 and 2015 rupture zones from the western part of the 1255 rupture zone. The relationship between the rupture area size and seismic moment of the 2015 earthquake lies between the scaling formulas for crustal earthquakes and plate-boundary earthquakes, but is closer to the former. Therefore, using this and the scheme for characterized source models, we determined realistic rupture parameters. We then simulated broadband ground motions in Kathmandu using these rupture parameters, our 3-D velocity structure models, and a hybrid method combining the finite-difference method and the stochastic Green’s function method. We obtained the peak ground accelerations (PGAs) of simulated ground motions, and calculated the seismic intensities in the Modified Mercalli intensity scale from the PGAs as indexes of hazards for Kathmandu. Intensities IX coincide with the center of the Kathmandu Valley, and intensities VIII and VII are found in the area surrounded by the sedimentary boundary and the southernmost part of the valley.


Introduction
The collision of the Indian and Eurasian plates makes the Himalayas a seismically active region.This collision causes the Indian lithosphere to underthrust the Himalayas along the Main Himalayan Thrust (MHT), as shown in the inset of Fig. 1.A large-scale branch fault of the MHT reaches the ground surface at the Main Frontal Thrust (MFT; red lines in the main map of Fig. 1), as also shown in the inset.Damaging earthquakes have been generated by this mechanism in the past.
Nepal is located in the center of the Himalayas (Fig. 1).Kathmandu is its capital and has a high population density of 5,169 people per square km (National Statistics Office 2023) and high urban vulnerability as discussed by Bhattarai and Conway (2010).Therefore, Nepal must anticipate such a large earthquake that could cause damage to Kathmandu in the future.In this study, we first construct a new distribution of the rupture zones of large historical earthquakes along the MHT and find the likely location of a future scenario earthquake.We then model its rupture parameters, simulate broadband ground motions from it to Kathmandu, and calculate PGAs and seismic intensities as ground motion hazards.
Constructing a new distribution of rupture zones Bilham (2019) made a chronological list of significant Himalayan earthquakes.From this list, we extracted earthquakes with magnitudes of 7.8 or greater and location of Nepal (Table 1).Bilham (2019) also determined the rupture zones of listed earthquakes using the rule of thumb that intensity VIII or higher shaking signified typical intensities above a rupture zone.His results for the 1934 and 1505 earthquakes are explained below and included in our new distribution, but since we have results of a seismic source inversion for the 2015 earthquake, we use them to derive our own rupture zone.We use the results of paleoseismic studies to determine the rupture zones of the 1225 and 1344 earthquakes.
For the 1934 earthquake, which is called the Bihar-Nepal earthquake, the rupture zone was determined using the intensities compiled by Ambraseys and Douglas (2004; numbers in green boxes in Fig. 1) and Martin and Szeliga (2010).The outline of the zone is shown in yellow with the number "1934" in Fig. 2 (1934 rupture zone).Although intensities VIII were also observed in Kathmandu and around the Nepal-India border, they were assumed to be due to sediment amplification in the Kathmandu Valley and Indo-Gangetic Plain.Therefore, the western boundary is located to the east of Kathmandu, and southern boundary is to the northern edge of the Gangetic Plain, which is near the MFT.The northern boundary is vague because of few intensity observations.Bilham (2019) adopted the locking line (orange triangles in Fig. 1), which is the lower end of the interseismic locking area along the MHT and almost follows the 3,500 m elevation contour (Avouac 2003).
For the 1505 earthquake, intensity data are not available; however, accounts of its shaking and damage are available in historical sources of Tibet.As summarized by Bilham (2019), the historical sources compiled by Jackson (2002) indicate that the rupture zone of this earthquake was located close to Guge (79.8°E, 31.5°N), Purangs (81.2°E, 30.3°N), Globo (84° E, 29° N), and Mangyul Gungthang (85° E, 29.2° N).They also indicate that the region from Globo westward and southward, including the Kali Gandaki region of central Nepal was more severely shaken than Mangyul Gungthang.In the villages north of Kathmandu, such as Gyirong, houses collapsed but there were fewer casualties.These locations are shown in Fig. 1 with white circles.The severe shaking in the Kali Gandaki region implies that the northern boundary of the rupture zone is located between the region and Globo.The locking line is also located between them; thus, Bilham (2019) readopted it for the northern boundary of the rupture zone.The MFT was readopted for the southern boundary.Assumedly, the eastern boundary was located between the Kali Gandaki region and Gyirong where fewer casualties were experienced.The western boundary was assumed to be in the south of the westernmost shaking location of Guge.The rupture zone surrounded by these boundaries is shown in yellow with the number "1505" in Fig. 2 (1505 rupture zone).Bilham (2019) noted that shaking and damage were documented up to Agra, which was located at least 250 km from the assumed rupture zone (Figs. 1 and 2).This may have been due to surface wave amplification, because Agra is located on the other side of the Gangetic Plain.
The latest damaging earthquake in Nepal is the 2015 earthquake, which is called the Gorkha earthquake.Several authors conducted source inversions for this earthquake and obtained similar results.Among them, we adopted the results of Kobayashi et al. (2016) with the corrections of Miyakoshi et al. (2017).The violet rectangle in Fig. 1 represents the corrected source fault.Their slip distribution is displayed in color tone with the number "2015" in Fig. 2. Assuming that a slip of approximately 3 m, which is 150% of the average slip, or more is significant, the part of the source fault bounded by the yellow dashed lines in Fig. 2 can be considered the rupture zone of this earthquake (2015 rupture zone).The northern boundary roughly coincides with the locking line at an altitude of 3,500 m, similar to those of the 1934 and 1505 earthquakes.
Next, we determined the rupture zones of the 1225 and 1344 earthquakes using the results of paleoseismic studies.Mugnier et al. (2011), Wesnousky et al. (2017a), Sapkota et al. (2013), andWesnousky et al. (2017b) conducted trenching surveys at Koilabas, Tiberi, Sir Khola,  boundaries is shown in green with the number "1255" in Fig. 2. Evidence of the 1344 earthquake was discovered at Koilabas and Butwal, but not at the other eastern sites.For the same reason as above, the rupture zone of the 1344 earthquake can be located coinciding with the easternmost portion of the 1505 rupture zone (green area marked "1344" in Fig. 2).

Constructing a scenario earthquake
It was also recently found from Global Navigation Satellite System (GNSS) observations that the boundary between the Indian and Eurasian plates is strongly coupled from the southern boundary of the MHT to a depth of approximately 10 km, with almost negligible lateral variations in the coupling (Tabei et al. 2021).This implies that all regions along the MHT have similar rates of strain increase.The MHT is covered by the rupture zones of large historical earthquakes as shown in Fig. 2. Therefore, assuming that during a large earthquake most of the accumulated strain within its rupture zone is released, it is most probable that the rupture zone of the oldest previous event will rupture as a future scenario earthquake.Among the earthquakes whose rupture zones are shown in Fig. 2, the 1255 earthquake is the oldest; however, large earthquakes have already occurred in 1934 and 2015 within its rupture zone, reducing the accumulated strain.Thus, we removed the 1934 and 2015 rupture zones from   Miyakoshi et al. 2017).We assumed the rupture zone of the latter was an area between the two dashed lines in yellow.The rupture zones of the 1255 and 1344 earthquakes (green areas) were determined based on the results of the paleoseismic studies (blue squares along the MFT in red; Mugnier et al. 2011, Sapkota et al. 2013, Wesnousky et al. 2017a, b, Okumura et al. 2021).A scenario earthquake for seismic hazard assessment in Kathmandu is assumed to occur beneath the orange rectangle the western part of the 1255 rupture zone, obtaining the area indicated by an orange rectangle in Fig. 2. As this area is close to Kathmandu, we assumed that a future scenario earthquake with seismic hazards for Kathmandu would occur there.
To construct a source model for the scenario earthquake, we first investigated the source characteristics of the 2015 Gorkha earthquake using the trimming procedure of Somerville et al. (1999).Its corrected source fault, comprising 17 × 10 subfaults (10 × 10 km), was trimmed using the following criterion.An edge row or column of a source fault is removed if the average slip per fault element in the entire row or column is less than 0.3 times the average slip of the whole fault.The edge row or column that has the lowest slip per fault element is removed first, and thereafter the process is repeated until all edge rows and columns have normalized average slip per fault element of 0.3 or greater.The procedure resulted in a source fault 120 km long and 80 km wide.We considered its area of 9,600 km 2 as the rupture area size of the 2015 Gorkha earthquake.Kobayashi et al. (2016) obtained a seismic moment of 7.4 × 10 20 N⋅m and a moment magni- tude M w (Kanamori 1977) of 7.8 for this earthquake.Somerville et al. (1999) determined a relationship between rupture area sizes S (km 2 ) and seismic moments M 0C in dyne⋅cm of crustal earthquakes.This resulted in 0C , which is equivalent to a formula with M 0 in N⋅m: For plate-boundary earthquakes in subduction zones, Somerville et al. (2002) obtained The rupture area size and seismic moment of the 2015 Gorkha earthquake are plotted with a red star symbol in Fig. 3.This falls between the relationships for crustal earthquakes and plate-boundary earthquakes in subduction zones (green and blue lines in Fig. 3), but is very close to the former relationship.This implies that an earthquake along the MHT has similar source characteristics to those of a crustal earthquake.Thus, we constructed a characterized source model (Miyake et al. 2003) for the scenario earthquake based on the formula (1).

Characterized source model
The characterized source model (Miyake et al. 2003) is based on the asperity model of Kanamori (1978), where an earthquake source is represented by a heterogeneous (2) distribution of stress drop on a fault (Madariaga 1979).In order to simulate realistic ground motions, we constructed this characterized source model for the scenario earthquake.
The fault related to an earthquake is called the rupture area in the characterized source model.As we considered it significant to avoid underestimating seismic hazards in Kathmandu, the rupture area of the scenario earthquake was assumed to extend under the whole of the orange rectangle in Fig. 2. As shown in the cross-sections of Fig. 4, beginning near the ground surface, the rupture area extends along the fault of MFT to the surface of the Indian lithosphere, and thereon, extends along the surface to the deeper MHT under the northern side of the rectangle.The upper part along the MFT is 8 km long and its dip angle is 35°, while the lower part in the MHT is 56 km long and its dip angle is equal to that of the Indian lithosphere, which is 7°.Therefore, the width of the rupture area is 8 km plus 56 km, which is 64 km.We can use the strike 289° and length 220 km of the orange rectangle for those of the rupture area.We assumed pure reverse faulting as the focal mechanism of the scenario earthquake, so that the rake was set at 90°.These geometrical parameters of the rupture area are summarized in the top rows of Table 2.
The size of the rupture area is 220 km × 64 km, which is 14,080 km 2 , and represented using subfaults 10 km long and 8 km wide.Using this rupture area size S and the formula (1), we obtained M 0 = 15.9 × 10 20 N⋅m.The M w for this M 0 is 8.1.As the cross-sections in Fig. 4 show that most of the rupture area is located in the darkgreen layer, we used the properties of this layer for the S-wave velocity β, density ρ, and rigidity μ in the source region.The average slip D defined through M 0 = µDS is 3.4 m.Brune (1970Brune ( , 1971) ) proposed the formula M 0 = 16/7 • a 3 �σ for a circular rupture area with the radius a and constant stress drop �σ .We calculated the average stress drop � σc assuming that this formula was valid for the scenario earthquake with a = (S/π) 1/2 and �σ = � σc .f max , which is a high-frequency limit of seis- mic source radiation (Hanks 1982), was set to 6 Hz based on Fujiwara et al. (2009).The parameters in this paragraph are listed in the lower part of Table 2.
The rupture area of a characterized source model comprises a few asperities with large slip and background region with less slip (Miyake et al. 2003).In addition to the formula (1), Somerville et al. (1999) proposed the constant ratio of the combined area of asperities S a to the rupture area size S to be 22%.According to this ratio, the S a for the scenario earthquake should be 3,098 km 2 ; however, it cannot be represented by the 80 km 2 subfaults mentioned above.Therefore, we slightly increased it to 3,200 km 2 .Somerville et al. (1999) also proposed the average number of asperities to be 2.6.We adopted the number of 2 by removing the fraction from their proposal.Somerville et al. (1999) further proposed a constant contrast of the average slip in the asperities ( D a ) to that in the whole rupture area ( D ) to be 2.01.Using this, we obtained D a of 6.8 m.Madariaga (1979) identified the relationship of average stress drops in the asperities ( �σ a ) and the whole rupture area ( � σc ) to be � σc = �σ a S a /S .This results in �σ a = 10 MPa.We assumed that the effective stress in the asperities ( σ ea ), which is necessary for computation of short-period ground motions, is equal to �σ a .Irikura and Miyake (2001) proposed the ratio of asperity areas in a two-asperity model to be 16:6.For the scenario earthquake, the area of the larger asperity S a1 should be 2,327 km 2 ; however, it cannot be represented by the 80 km 2 subfaults.We removed the fraction smaller than 1,000 km 2 so that S a1 is 2,000 km 2 comprising 5 × 5 sub- faults.The area of the smaller asperity S a2 is the remainder of S a , which is 1,200 km 2 comprising 5 × 3 subfaults.These asperities are located in the middle of the rupture area so that they are mostly in the V S 3.5 km/s layer of the MHT as shown in Fig. 4.Both asperities have the same stress drop �σ a .Dan et al. (2001) proposed that where D ai stands for a slip on the i-th asperity and M 0a = µD a S a .Introducing r a , r ai , and γ i such that S a = πr a 2 , S ai = π r ai 2 , and γ i = r ai /r a , we obtained from (3).Using (4), we identified D a1 = 7.4m and D a2 = 5.8m .Next, we identified M 0a1 = 4.9 × 10 20 N⋅m, M 0a2 = 2.3 × 10 20 N⋅m, and M 0a = 7.2 × 10 20 N⋅m using M 0ai = µD ai S ai and The parameters for the asperities are summarized in the upper part of Table 3.The area of the background region S b is S minus S a , which is 10,880 km 2 .Miyake et al. (2003) assumed no stress drop there, but some slip does occur because of �σ a in the asperities as shown by Madariaga (1979).This slip generates short-period ground motion and there is effective stress σ eb corresponding to it.The average of the slip ( D b ) was calculated to be 2.4 m using M 0b = µD b S b where M 0b = M 0 − M 0a = 8.7 × 10 20 N⋅m.We calculated σ eb to be 2.0 MPa using the relation σ eb = 0.2 • �σ a , which was obtained by Miyatake (2002) from the results of dynamic rupture simulations for the (3) asperity model.These parameters for the background region are summarized in the lower part of Table 3.
For the Gorkha earthquake, the hypocenter was located in the westernmost part of the rupture zone, and the rupture propagated eastward from the hypocenter with a velocity of 3.3 km/s (Kobayashi et al. 2016).Similarly for the scenario earthquake, we assumed that the rupture initiation point, which is equivalent to a hypocenter, was within the westernmost part of the asperities that are equivalent to a rupture zone.In addition, we assumed that the rupture propagated eastward with the same velocity, which is shown in the lowermost row of Table 2 as rupture velocity v r .The rupture initiation point was assumed to be located at the southwestern corner of the larger asperity.If the rupture propagates eastward from this location, it will head toward Kathmandu, so that the directivity effect that amplifies ground motion (Koketsu et al. 2016) is maximized there.Therefore, based on our policy of avoiding the underestimation of seismic hazards, we chose this location for the rupture initiation point (star symbol in Fig. 4).

Ground motion simulation
Irikura and Miyake (2001) suggested using the hybrid Green's function method for simulating ground motions from a future earthquake, such as the scenario earthquake, and cited Kamae et al. (1998).However, their method applies to a source model that comprises a few Fig. 3 Relationships between rupture area size and seismic moment for crustal earthquakes (green line; Somerville et al. 1999) and for plate-boundary earthquakes in subduction zones (blue line; Somerville et al. 2002).The 2015 Gorkha earthquake is plotted with the red star symbol, which falls between the two relationships, but it is located close to the relationship for crustal earthquakes.In addition, the relationship of Murotani et al. (2008) is also plotted with a blue dashed line subevents and does not include a background region.Fujiwara et al. (2009) extended it to such a source model comprising subfaults as our characterized source model.In this extended method, Green's functions refer to ground motions from subfaults, which are calculated using deterministic and stochastic approaches.
To calculate the long-period component of ground motion, we applied the finite-difference method (FDM) as a deterministic approach, following Kamae et al. (1998) and Fujiwara et al. (2009).Our FDM is based on a fourth-order difference scheme in the 3-D wavefield (Levander 1988), a staggered grid (Virieux 1986), the viscoelastic formulation of Robertsson et al. (1994), and the absorbing boundary layer of Cerjan et al. (1985).As a velocity structure model, we adopted the regional-scale model of Koketsu et al. (2024), where the V S 0.44 km/s and 2.9 km/s layers were considered the engineering bedrock and the basement, respectively.However, owing to computational limitations, the FDM was applied to the engineering bedrock and deeper parts.As the FDM does not support free-surface topography, the squashing method explained in Aagaard et al. (2008) was applied to deform the engineering bedrock and layers immediately below it in the vertical direction.This was performed so that the free surface, which is the top surface of the engineering bedrock, becomes flat and aligned at the sea level (Fig. 5 upper).We used a variable grid, in which the vertical interval varied from 100 to 600 m.The horizontal interval was fixed at 100 m.This grid was applied to the flattened model.Point sources with the focal mechanisms shown in Table 2 were located at the centers of subfaults in the asperities and background region of the characterized source model (Fig. 5 lower).Their seismic moments were calculated by dividing M 0a1 , M 0a2 , and M 0b by N a1 , N a2 , and N b , respectively.These parameters are listed in Table 3. Nakamura and Miyatake (2000) proposed the physically reasonable time function of a slip rate dD(t)/dt (Fig. 6): where ( 5) Day (1982) called t r the rise time and obtained the second and last empirical formulas in (6).Once the four parameters in (6) are determined, the remaining parameters can be determined to be the functions of t b , represented as Numerically solving the equation  6), (7), and ( 8) have different values depending on the asperities and background region, so three slip rate functions were obtained.It is noteworthy that W b was assumed to be the largest width of the background region, which is equal to the width W of the rupture area.We then performed the FDM computations to obtain ground velocities at the top of the engineering bedrock in the Kathmandu Valley using the obtained slip rate functions and the source representation of Graves (1996) for a point source.
For the short-period component of ground motion, we applied the stochastic Green's function method as a stochastic approach.This choice is again the same as that in Kamae et al. (1998) and Fujiwara et al. (2009).Combining the results of Aki (1967), Brune (1970), andHanks (1982), Boore (1983) proposed the spectrum of seismic ground acceleration at the basement from a point source: we used ρ and β in Table 2.For the seismic moment M 0 , M 0a1 /N a1 , M 0a2 /N a2 , and M 0b /N b were used as in the FDM computation.R θφ in the first term represents the radiation pattern and was calculated using the method of Kagawa (2004).f c in the second term is the corner frequency calculated as which was obtained from Brune (1970Brune ( , 1971)), as shown in Text S1. m in the third term was obtained to be 4.2 by Satoh et al. (1997) and we used m = 4 by removing the fraction.The whole path attenuation is accounted for by the fourth term with the distance R from a point source to a target point and ( 9) of Satoh et al. (1997).The last term represents the amplification due to the impedance contrast of a point source and the basement with ρ b = 2.5 g/cm 3 and β b = 2.9 km/s.Satoh et al. (1994) studied the envelopes of ground accelerations observed in boreholes using the function of Jennings et al. (1968): and obtained scaling relations for the time differences in (12).Based on the formulas of Sato (1989) and Ohsaki (1984), they were rewritten as ( 12) with ( 13) was adopted for the envelopes of ground accelerations at the basement.The M 0 and R in ( 13) are the same as for (9).Following Boore (1983), we first generated a time sequence of white Gaussian noise using a pseudo-random number generator and a seed for it, and windowed it with (12).We next Fourier transformed the windowed time sequence into the complex spectrum, which was then multiplied by (9).The result of the multiplication was inversely Fourier transformed into a time sequence of ground acceleration at the basement.This process was repeated using various seeds of the pseudorandom number generator until a time sequence with a reasonable waveform was obtained.
We extracted the horizontally layered structure between the basement and engineering bedrock beneath a target point, from the 3-D velocity structure of Koketsu et al. (2024).Using this structure and such a matrix method as that of Thomson (1950), we obtained the ground acceleration at the engineering bedrock from that at the basement.This is called the stochastic Green's function (SGF).

Ground motion synthesis
The FDM computation is so time-consuming that its repetition is not practical.Therefore, for the computation of the long-period component, we introduced all the point sources into the flattened velocity structure.We directly obtained ground velocities at the engineering bedrock including their contributions in one computation.Ground accelerations were obtained by numerically differentiating them, and shortly called FDM accelerations.For the short-period component, we applied the SGFs of the point sources.They were summed up using the method of Irikura et al. (1997), which was explained in Miyake et al. (2003).This summation resulted in ground accelerations at the engineering bedrock, which were called SGF accelerations.We synthesized broadband ground accelerations at the engineering bedrock using a pair of matched filters as proposed by Kamae et al. (1998).Our pair comprised a high-pass filter that allows 50% transmission at 0.5 Hz and 100% transmission above 0.6 Hz, and a low-pass filter that blocks 50% transmission at 0.5 Hz and 100% transmission above 0.6 Hz (Fig. 7a).The former and latter were applied to the SGF and FDM accelerations, respectively.The filtering results were summed in the frequency domain and inversely Fourier transformed to obtain broadband ground accelerations.Such a process is called hybrid synthesis, because the obtained accelerations are hybrids of SGF and FDM accelerations (e.g., Kamae et al. 1998).For the NS, EW, and UD components at the engineering bedrock beneath the point P indicated by the magenta dot in Fig. 7b, the accelerations before and after the hybrid synthesis are shown in Fig. 7c with the labels of SGF, FDM, and Hybrid as an example.The peak ground acceleration (PGA) was obtained from the amplitude of a Hybrid waveform.
The above was repeated for all points within the area surrounded by the basin boundary of Shrestha et al. (1998), which is shown in Figs.7b and 8a by the brown lines.In addition to the point P, the broadband ground accelerations at the points A and B are shown in Fig. 8b as further examples.The distributions of the obtained PGAs are illustrated in Fig. 9a.We then calculated the PGAs at the ground surface by adding the effects of the shallowest part of the velocity structure to those at the engineering bedrock.JICA (2018) examined this part and summarized the results in the distribution of AVS30 (average S-wave velocity of the upper 30 m; Borcherdt et al. 1979), whose values range from 158 to 597 m/s with an average of 262 m/s within the sedimentary boundary.Fujimoto and Midorikawa (2006) proposed the empirical relationship between AVS30 and amplification of the PGA in the shallowest part of a velocity structure: where F is an amplification factor for the PGA, and β bedrock is an S-wave velocity at the engineering bedrock, which is 440 m/s in this study.γ eff stands for an effec- tive strain representing the degree of nonlinear response of the shallowest part.This was assumed to be 1 × 10 −3 in this study.If AVS30 was greater than β bedrock , which implied that the engineering bedrock or a deeper part was exposed, F is equal to unity.The PGA at the engi- neering bedrock was multiplied by F to obtain the PGA at the ground surface.The PGA distribution obtained is shown in Fig. 9b.Finally, we calculated seismic intensities at the ground surface based on the results obtained above adopting the Modified Mercalli (MM) intensity scale and the method of Wald et al. (1999).In this method, the PGAs are used with the formula  9a, the highest intensities of IX were first caused by large ground motions from the scenario earthquake and their amplification by the intermediate velocity structure between the basement and engineering bedrock.The second cause was the amplification by the shallowest velocity structure, as shown in Fig. 9b.Since the waveforms in Figs.7c and 8b show relatively short wavepackets rather than long wavetrains, we can ascertain that the effect of the intermediate velocity structure is not caused by the development of long-period ground motions in a basin (e.g., Koketsu (15) MM intensity = 3.66 log (PGA) − 1.66 .
7 The pair of matched filters for the SGF and FDM accelerations (black and gray lines) and the topographic map for the Kathmandu Valley (a and b).c shows the SGF and FDM accelerations at the engineering bedrock (black and gray traces labeled "SGF" and "FDM") at the point P indicated by a magenta dot in b.The blue traces labeled "Hybrid" represent the ground accelerations obtained by their hybrid synthesis using the filters in a.The upper, middle, and lower sets of traces are for the NS, EW, and UD components, respectively.SGF and FDM stand for the stochastic Green's function method and the finite-difference method, respectively and Miyake 2008), but by the reverberation amplification of body waves.
Furthermore, we found intensities VIII and VII in the southernmost part of the Kathmandu Valley, which is outside the sedimentary boundary (Fig. 9d).This part is located at the periphery of the valley (Fig. 9c), so that the effect of the velocity structure is considered insignificant.Therefore, the large intensities should have been mainly caused by larger ground motions at shorter distances from the asperities of the scenario earthquake.

Variation in results
Among the parameters included in our method, the ones that have the greatest influence on the results are the magnitude of the seismic moment, the location of the rupture initiation point, and the location of the asperities.We here discuss the variation in results due to the former two parameters.Regarding the variation due to the third parameter, Iwaki et al. (2017) discussed it using crustal earthquakes in Japan.
Our source modeling started with the location and size S of the rupture area shown in Fig. 2 and Table 2.We then determined M 0 to be 15.9 × 10 20 N⋅m from S = 14,080 km 2 and the formula (1), but M 0 can be brought closer to the formula (2) considering the variation.Murotani et al. (2008) proposed a relationship that lies between the formulas (1) and ( 2), so we used it to assume another M 0 of 9.3 × 10 20 N⋅m.We con- structed a characterized source model based on this M 0 and called it Model M, while the original model was called Model S. Figure 10a, which is the same as Fig. 9b, shows the PGA distribution for Model S, and Fig. 10b shows that for Model M. Comparing them, the distribution patterns are similar, but the PGA of Model M is generally smaller than that of Model S.
We next varied the rupture initiation point.As shown in Fig. 10d, the original one is in the west, while the varied one is in the east.Figure 10c shows the PGA distribution for Model S with the eastern rupture initiation point.When we compared this with Fig. 10a, we found that the PGA for eastern initiation was generally smaller and the distribution pattern was different from that for western initiation.The reason for these is that Kathmandu is not in the forward direction with respect to the eastern initiation point, so the directivity effect is small and the pattern is different from that for the western initiation point.

Conclusions
We constructed the new distribution of the rupture zones of great historical earthquakes along the MHT.Based on the finding of Tabei et al. (2021) that there is almost no lateral change in the plate coupling of the MHT, we considered it most probable that the rupture zone of the oldest previous event in the distribution would rupture as a future scenario earthquake.Based on this principle, we obtained the rupture zone of a scenario earthquake for the assessment of the seismic hazards in Kathmandu, and constructed the characterized source model for the scenario earthquake.We simulated broadband ground motions in Kathmandu applying this source model, our 3-D velocity structure models, and a hybrid method combining the finite-difference and stochastic Green's function methods.In the simulation results, intensities IX in the MM intensity scale coincide with the center of the Kathmandu Valley.We identified intensities VIII and VII in the area surrounded by the sedimentary boundary and in the southernmost part of the Kathmandu Valley.The intensities outside these areas are VI and V.There are many houses and buildings of weak and ordinary masonry in Kathmandu.Regarding intensity IX, Richter (1958) stated that weak masonry is destroyed, ordinary masonry is heavily damaged, and general panic occurs.Ground motion of intensity VIII damages ordinary masonry, which can lead to partial collapse.Similarly, ground motion of intensity VII damages weak masonry including cracks.It is crucial for Kathmandu to prepare for such ground motion hazards.

Discussion
The method used here, which is a combination of the characterized source modeling and the hybrid Green's function method, has been verified by Subcommittee for Evaluation of Strong Ground Motion (2002Motion ( , 2008) ) and Iwaki et al. (2016).The advantage of the hybrid Green's function method is that it can generate not only PGAs and seismic intensities, but also realistic ground motion time histories as shown in Figs.7 and 8. Furthermore, since response spectra used to evaluate the impact on houses and structures can be calculated from the time histories as shown in Fig. 8c, analyses of ground motion hazard can be linked to analyses of ground motion risk.
Figure S1 compares the simulated intensity distribution for the scenario earthquake (Fig. 9d) with the observed intensity distribution of the Gorkha earthquake (Martin et al. 2015).The maximum intensities for the Gorkha and scenario earthquakes are VIII and IX in the MM intensity scale, respectively.This difference may be due to the difference in seismic moments (7.4 and 15.9 × 10 20 Nm).The largest intensities of VII and VIII for the Gorkha earthquake are mainly distributed in the northern part of the basin from west to east, probably because large slips of the earthquake were adjacent to this part as shown in Fig. 2. On the other hand, those of VIII and IX for the scenario earthquake are distributed in the center of the basin.Since large slips of the earthquake were apart from Kathmandu as shown in the lower panel of Fig. 10d, the large intensities may have been mainly caused by the distribution of thick sediments.

Fig. 1
Fig. 1 Index map for historical seismicity in Nepal and neighboring regions.Kathmandu is identified by the brown boundary of the Kathmandu Valley (Shrestha et al. 1998).The inset shows a cross-sectional view of the velocity structure by Koketsu et al. (2024) along the section line XY in cyan.The orange triangle represents the intersection point of the section line and the locking line at an altitude of 3,500 m.Paleoseismic trenches were conducted at the locations marked by the blue squares.The red lines in the map and dashed black line in the inset indicate the Main Frontal Thrust.The white circles indicate the locations of heavy damage by the 1505 earthquake.Seismic intensities of VIII or higher observed during the 1934 earthquake are plotted with Arabic numbers in green boxes.A source inversion was performed for the 2015 Gorkha earthquake using the corrected source fault shown with the violet rectangle

Fig. 2
Fig.2Rupture zones of great historical earthquakes in1505 and 1934 (yellow areas;Bilham 2019) and the corrected slip distribution of the 2015 Gorkha earthquake (color tones;Kobayashi et al. 2016;Miyakoshi et al. 2017).We assumed the rupture zone of the latter was an area between the two dashed lines in yellow.The rupture zones of the 1255 and 1344 earthquakes (green areas) were determined based on the results of the paleoseismic studies (blue squares along the MFT in red;Mugnier et al. 2011, Sapkota et al. 2013, Wesnousky et al. 2017a, b, Okumura et al.  2021).A scenario earthquake for seismic hazard assessment in Kathmandu is assumed to occur beneath the orange rectangle

Fig. 4
Fig. 4 Left: Bird's eye view of the characterized source model for the scenario earthquake.The white rectangle and gray squares represent the rupture area and asperities of the model, respectively.The rupture area comprises an upper part with high dip angle and a lower part with low dip angle.They are connected at the light-green line.The star symbol represents the rupture initiation point.The red lines indicate surface faults of the MFT.The brown boundary indicates the Kathmandu Valley (Shrestha et al. 1998).Right: cross-sectional views of the model along the sky-blue lines A and B on the left.The background shows the velocity structure byKoketsu et al. (2024)

Fig. 5
Fig.5Lower: geometrical configuration of the characterized source model comprising the asperities (gray squares) and background region (white part).They are divided by blue lines into 22 × 8 subfaults 10 km long and 8 km wide with point sources at their centers (small black dots).The star symbol represents the rupture initiation point.The brown boundary surrounds the Kathmandu Valley(Shrestha et al. 1998).Upper: cross-sectional view of the rupture area with an aspect ratio of 1:1 along the cross-section line (sky blue) in the lower.The background is the flattened velocity structure model.We show the part between 10 and 90 km from the southern end of the cross-section line.The blue dots represent the intersections of the cross-section line and dividing lines, indicating the upper and lower sides of the subfaults.According to the aspect ratio, it is shown that a subfault along the MFT has the same width as those of subfaults in the MHT Figure9dpresents the distribution of the MM intensities calculated based on the PGAs at the ground surface.Intensities IX were obtained in the center of the Kathmandu Valley.We identified intensities VIII and VII in the area surrounded by the sedimentary boundary ofShrestha et al. (1998), which is indicated in Fig.9cby the black line.As the largest PGAs at the engineering bedrock were also found in the valley center as shown in Fig.9a, the highest intensities of IX were first caused by large ground motions from the scenario earthquake and their amplification by the intermediate velocity structure between the basement and engineering bedrock.The second cause was the amplification by the shallowest velocity structure, as shown in Fig.9b.Since the waveforms in Figs.7c and 8bshow relatively short wavepackets rather than long wavetrains, we can ascertain that the effect of the intermediate velocity structure is not caused by the development of long-period ground motions in a basin (e.g., Koketsu

Fig. 8 a
Fig. 8 a Topographic map for the Kathmandu Valley including the representative points (magenta dots) in downtown Thamel (A) and the basin center (B).b Three components of ground acceleration simulated at A and B in a. c Pseudo-velocity response spectra of the NS and EW components (Left and Right) of ground acceleration at A, B, and P in a. P is the same as in Fig. 7.The basin and sedimentary boundaries defined by Shrestha et al. (1998) are shown in a with brown and black lines, respectively

Fig. 9
Fig. 9 The distributions of the PGAs at the engineering bedrock and ground surface (a and b).d is the distribution of MM intensities at the ground surface.The brown boundaries in all the distributions are the basin boundary defined by Shrestha et al. (1998) as in the topographic map c, where the sedimentary boundary is also shown in black

Fig. 10
Fig. 10 The distributions of the PGAs for Models S and M with the original rupture initiation point in the west (a and b).c is that for Model S with the eastern rupture initiation point.Models S and M, and the both rupture initiation points are shown in d.The brown boundaries in all the panels are the basin boundary defined by Shrestha et al. (1998)

Table 1
Damaging earthquakes in Nepal with magnitudes of 7.8 or larger (from Table 1 of Bilham 2019)

Table 2
Parameters for the whole rupture area of the characterized source model for the scenario earthquake

Table 3
Parameters for the asperities and background region of the characterized source model for the scenario earthquake finally gives t b .As shown in Table3, W i , σ ei , and D i in ( Model of slip rate function (based on Nakamura and Miyatake 2000, reprinted from Koketsu 2018 with permission of Kindai Kagaku)