BULLETIN OF THE SEISMOLOGICAL SOCIETY OF AMERICA, VOL. ???, XXXX, DOI:10.1785/,

Three-Dimensional Fault Structure Inferred from a Refined
Aftershock Catalog for the 2015 Gorkha Earthquake, Nepal
Masumi Yamada,1 Thakur Kandel,2 Koji Tamaribuchi3 Abhijit Ghosh4
Abstract.

In this paper, we created a well-resolved aftershock catalog for the 2015 Gorkha earthquake in Nepal by processing 11 months of continuous data using an automatic onset
and hypocenter determination procedure. Aftershocks were detected by the NAMASTE
temporary seismic network that is densely distributed covering the rupture area and became fully operational about 50 days after the mainshock. The catalog was refined using a joint hypocenter determination technique, and an optimal one-dimensional (1D)
velocity model with station correction factors determined simultaneously. We found around
15,000 aftershocks with the magnitude of completeness of M2. Our catalog shows that
there are two large aftershock clusters along the north side of the Gorkha–Pokhara anticlinorium and smaller shallow aftershock clusters in the south. The patterns of aftershock distribution in the northern and southern clusters reflect the complex geometryof the Main Himalayan Thrust. The aftershocks are located both on the slip surface andthrough the entire hanging wall. The 1D velocity structure obtained from this study isalmost constant at a P-wave velocity (Vp) of 6.0 km/s for a depth of 0–20 km, similar
to Vp of the shallow continental crust.

Introduction

et al., 2017] of the MHT shear zone. Aftershocks identified
using far-field data are strongly influenced by the velocity
structure, and the data are limited to large earthquakes.
Variations in the local subsurface structure also result in
estimation errors.
Another uncertainty is the three-dimensional (3D) structure of the MHT. Nepal lies on the collision zone between
the Indian Plate and the Eurasian Plate, and the Indian
lithosphere underthrusts beneath the Himalayas along the
MHT. Three major north-dipping faults formed because ofthis collision—the Main Frontal Thrust (MFT), the Main
Boundary Thrust (MBT), and the Main Central Fault(MCT)—from the south to the north. Figure 1 shows thegeological map of Nepal modified after Dahal [2006]. The
north–south (NS) section, in the direction of motion of the
Indian Plate, has been thoroughly studied using microseismicity [e.g. Pandey et al., 1995, 1999] and geophysical tools
such as receiver function [e.g. Tilmann et al., 2003; Nábelek
et al., 2009]. Although it is believed that the variation in
the east–west (EW) direction is largely homogeneous, thestructure is not well studied.
In order to improve the 3D spatial resolution of the aftershock activity, we used the temporary seismic network,
called NAMASTE [Karplus et al., 2015; Pant et al., 2016;
Mendoza et al., 2016; Ghosh et al., 2017; Karplus et al.,2017; Bai et al., 2019; Mendoza et al., 2019]. 47 near-source
stations were deployed; in comparison, the permanent seismic network has only seven stations in the same area. Thisdataset enables us to resolve and locate aftershocks down tomagnitude 2.0. Based on this well-resolved aftershock catalog, we will discuss the relationship between the aftershocks
and the 3D fault structure in central Nepal.

On April 25, 2015, an earthquake (Mw 7.8) in central
Nepal (later referred to as the Gorkha earthquake) caused
significant damage, including 8979 fatalities. This was thelargest earthquake that occurred after the installation of
the national seismic network in Nepal. The earthquake ruptured the Main Himalayan Thrust (MHT), the plate boundary between India and the Tibetan Plateau [Galetzka et al.,2015; Avouac et al., 2015].
The aftershock activity was studied immediately after theearthquake. Adhikari et al. [2015] manually processed theseismic data of the permanent national seismic network in
Nepal and presented a temporal and spatial distribution of
the early aftershocks. Baillard et al. [2017] processed thesame data with automatic onset and hypocenter determination procedures and created an aftershock catalog for thethree months following the earthquake. Bai et al. [2016]used a seismic array near the China–Nepal border and located early aftershocks using the double-difference earthquake relocation technique. Aftershock activities were alsostudied using global seismic data [e.g. Letort et al., 2016;
Wang et al., 2017]
Although the horizontal distribution is roughly consistentamong these studies, the hypocenter depths exhibit a larger
uncertainty, making it difficult to clarify the relationshipbetween the aftershocks and the fault structure. The mainshock rupture surface is located on the MHT décollement
based on the hypocenter depth and moment tensor mechanism [Adhikari et al., 2015; McNamara et al., 2017; Baillard
et al., 2017]. However, the aftershocks may be shallower[Bai et al., 2016; Baillard et al., 2017], deeper [McNamara
et al., 2017], or on the boundary [Letort et al., 2016; Wang

Data

1 Disaster Prevention Research Institute, Kyoto

University, Uji, Gokasho, 611-0011, Japan
2 Department of Mines and Geology, Kapurdhara Marg-9,
Lainchour, Kathmandu, Nepal
3 Meteorological Research Institute, 1-1 Nagamine,
Tsukuba, Ibaraki 305-0052, Japan
4 University of California, Riverside, Riverside, CA, 92521,
US

About 50 days after the Gorkha earthquake a temporary seismic network called NAMASTE was deployed in theepicentral area of the mainshock to record the aftershocks[Karplus et al., 2015; Pant et al., 2016; Mendoza et al., 2016;
Ghosh et al., 2017; Karplus et al., 2017; Bai et al., 2019;
Mendoza et al., 2019]. This network included broadband
seismometers, short-period seismometers, and accelerometers. For our research, we used 42 broadband and shortperiod records covering the period from June 25, 2015, to
May 14, 2016. The average spacing of the stations was

Copyright 2020 by the American Geophysical Union.
0886-6236/20/$12.00

1

X - 2 YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE
about 20 km. The three component continuous waveforms
were downloaded from the IRIS Data Management Center
website.

Methods

where A was the maximum amplitude of the vertical displacement in that time window (nanometers) and R was
the distance from hypocenter to station in kilometers. The
median value of the five stations is chosen as the event magnitude.
Refinement of the Phase Detection

Phase Detection by the T pd Method
Seismic waveforms were processed by first correcting the
instrumental response and removing the DC offset. Then,
we applied a fourth-order one-pass band-pass filter with acorner frequency of 2–10 Hz. Phase arrivals were detectedusing the T pd method [Hildyard et al., 2008; Hildyard andRietbrock, 2010]. This method identified seismic phase arrivals by observing the change in frequency of the waveforms. After applying the T pd method, the arrival time wasrefined by following the methods of Yamada [2017]. Figure
S1 shows an example of the automatic phase detection.
Once a phase arrival was detected, the refined arrival
time and the maximum vertical velocity 5 s after the P-wave
onset were computed. They were used for the hypocenter
determination. The ratio between the maximum horizontal amplitude and the maximum vertical amplitude (H/V)
in 1 s was used to classify P- and S-waves. If H/V < 0.5
or H/V > 2, the phase was classified as a P-wave or Swave, respectively. Otherwise, the phase was treated as anindeterminate type of wave, and the automatic determination algorithm (explained in the next subsection) identified
the phase that best fits the travel-time curve. Note that
this picking program was tuned to detect P-wave arrival, so
there were more P-wave observations than S-wave observations.
Automatic Hypocenter Determination
We implemented an automatic hypocenter determination
method using Bayesian estimation [Tamaribuchi, 2018] tocreate an earthquake catalog. The program read a set of arrival times and maximum vertical amplitudes and createda subset of observations that were matched to a possible
earthquake location. The closest 10 stations from the first
trigger station were defined as a trigger group. Once threeor more stations in the group were triggered within the theoretical travel-time limit, the hypocenter calculation began.
In order to determine the location, virtual hypocenters(i.e., particles) were distributed in 3D space. The weight ofeach particle (i.e., the likelihood that the particle represents
the actual hypocenter) was computed by taking the inverseof the difference between the observed and the theoreticalvalues. We assumed that this error followed the normal distribution. Therefore, if the location was correct, the error
between the observation and the theoretical value would besmall and the weight would be large.
We then used an importance sampling algorithm to refine the hypocenter location [Wu et al., 2014]. A resampling
process was performed so that the particle distribution reflected the probability density function of the hypocenter
location. When the likelihood calculation and resampling
were repeated, virtual hypocenters gradually converged on
the true hypocenter. The program stopped the iteration
process when the virtual hypocenter with the maximum
weight did not change after three successive iterations. It
produced a final hypocenter location only if all the following conditions were satisfied: 1) the longitude and latitude
errors were both less than 10 min (about 18 km), 2) theorigin time error was less than 2.0 s, 3) the depth error wasless than 100 km, and 4) 5 or more phases were used for thelocation determination.
After determining the hypocenters, we calculated themagnitude using the Nepalese attenuation model that was
proposed by Baillard et al. [2017]. We selected the five closest stations for each event and fixed the time window fromthe theoretical P-wave arrival time to three times the difference between the arrivals of the theoretical S- and P-waves.
We used the following equation to obtain the magnitude[Baillard et al., 2017]:
M = 0.9 log 10(A) + 1.2 log 10(R) + 0.0003R − 0.9,

(1)

In order to improve the accuracy in the depth direction,
we searched the P- and S-wave arrivals at around the theoretical arrival times based on the catalog we constructed.
We used the phase detection method with the variance ratio[Tamaribuchi, 2018].
First, we applied a band-pass filter with a corner frequency of 5–10 Hz for the three component waveforms to
reduce noise from tidal or distant events. We used waveforms for which either the P- or S-phase was detected by theT pd method. The two horizontal waveforms were rotated toradial and transverse directions. We used the vertical andradial components to detect P- and S-phases respectively.
Next, we computed the variance ratio rvar of the filteredwaveforms for a few seconds before and after the theoretical arrival time. The rvar was defined as:

X

t+N−1

rvar (t) =i=t(xi − x(t))2 /t−1
X(xi − x(t − N ))2 ,

(2)

i=t−N

where xi was the filtered waveform dataat the time step i,
Pt+N−1x(t) was the average of the data, i.e., i=t
xi /N , and N
was the length of the data. We used the number of samples
in one second for N . If the variance ratio exceeds 9 in the
few seconds time window, we computed the time recording
the peak of rvar , and used it as an arrival time.
Finally, we refined the P-wave arrival time by the autoregressive Akaike Information Criterion (AR-AIC) method[Takanami and Kitagawa, 1988, 1991]. In the AR-AIC
method, a few second data before and after the estimatedarrival time were treated as stationary noise and signal, respectively, and modeled as an AR process. The time seriesof the AIC value was expected to be the minimum at thephase arrival time. We applied this process only for the
P-wave arrival, since the method did not work well for the
S-wave.
After this phase detection process using variance ratio,
the number of the S-wave picking increased from 61,266
to 87,621, and that of the P-wave slightly changed from
133,088 to 138,436.
Joint Hypocenter Determination
We used the VELEST program to perform thejoint hypocenter determination (JHD) [Kissling et al.,
1994, 1995]. VELEST was a FORTRAN77 routine that simultaneously determines the location of a group of events,
the 1D velocity model, and station corrections that accountfor lateral velocity variations.
First, we selected relatively well-located earthquakes
from our catalog. The selected events were the earthquakes
where the number of picks is ≥20 and the azimuthal gapof the stations is ≤240 (see Result Section on these criteria). We defined these 3,304 events as well-located earthquakes. The waveform alignment of these events were visually inspected and we confirmed that they were all realearthquakes.
The initial velocity model was, as described by Pandey
et al. [1995], a three-layer model with boundaries at 23 and
55 km depths. We subdivided this model into layers witha 2 km interval for 0–10 km, and a 5 km interval for 10–60
km. We chose station NA280 as the reference station (i.e.,YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE
no station correction) since it was near the center of theaftershock distribution.
After performing the joint hypocenter determination, theprogram was run in single event mode to relocate all earthquakes with a fixed velocity model and station corrections.
In total, 14,721 events were detected for the 11 month period.

Results
Spatiotemporal Distribution of Aftershocks
Figure 2 shows the temporal and spatial distribution ofaftershocks. Large clusters are present along the northernside of the Gorkha–Pokhara anticlinorium (GPA), divided
by the seismic gap north of Kathmandu. The eastern cluster is very dense and spans about 100 km parallel to the
MBT fault line. The western cluster is sparse and scatters in the EW direction. Smaller clusters exist about 20
km west of Kathmandu and about 40 km east–southeast ofKathmandu. These clusters are much shallower than thelarge clusters.
Figure 2a shows the temporal variation of the aftershock
sequence. Since our catalog data starts after June 25, 2015,
we used a catalog created by Baillard et al. [2017] for earlier
dates. Only well-located events (with epicentral errors of <7
km and depth errors of <10 km) in the Baillard et al. [2017]
are used. Compared with the previous study, our catalog
has more earthquakes and shallower events. The seismicity
west of Kathmandu decays faster whereas the eastern seismicity around the largest aftershock is still active after 1
year.
In order to investigate the depth variation, the NS and
EW cross sections are shown in Figure 3. Aftershocks within±5 km from each section line are shown in the figure. The
cross sections show that the aftershocks are distributed in
the shallow area, between 0 and 20 km. In the NS section(Sections A1 and A2), the depth of the aftershocks becomesdeeper from the south to the north, which is consistent withthe dip direction of the MHT. For example, there are twoclusters in Section A2: the southern cluster is about 5 km
shallower than the northern cluster. Another feature is thatmost of the aftershocks are distributed in the hanging wall
between the ground surface and MHT fault boundary. The
depth variation in the EW section (Sections B1) is not as
obvious as the NS section. It seems the depth is about
constant in this direction. There is a gap with few events
between 80 and 100 km from the west end.
We plot the histograms of the magnitude of the catalog in Figure 4a. The magnitude distribution follows
the Gutenberg-Richter law, and the minimum magnitude
of complete recording (magnitude of completeness) is estimated about M2 from Figure 4a, which is smaller than inprevious studies (M4.0 for Adhikari et al. [2015] and M2.5
for Baillard et al. [2017]).
Estimated 1D Velocity Model
VELEST program simultaneously determines the location of earthquakes, the 1D velocity model, and station corrections. Figure 5 and Table 1 show the optimal velocity
model and station corrections. Additionally, we have performed this joint hypocenter determination with a different
initial model (JMA2001, Ueno et al. [2002]), but the results
are similar. The P-wave velocity (Vp) at depths <20 km
is almost constant, and estimated to be 6.0 km/s. This isslightly larger than Vp obtained in a previous study; 5.6
km/s in Pandey et al. [1995], 5.7 km/s in Monsalve et al.
[2006] for eastern Nepal, and 5.8 km/s in Negi et al. [2017]
for west of Nepal. This value is similar to Vp of the shallow
continental crust [Christensen and Mooney, 1995]. From 20
km depth, Vp gradually increased as the depth increased,
to 60 km depth, and we did not observe any strong velocity contrasts. However, we may not have enough resolution
at depths >20 km because of the limited number of deep
earthquakes in our dataset.

X-3

Station corrections have spatial variations, and values arelarger in the south (see Figure 5b). This southern part of
Nepal is a lowland area formed by Himalayan rivers, called
Terai. The surface structure consists of alluvial deposits,
so we would expect to see a slower Vp at shallow depth.
Larger values for station corrections are likely to reflect thislow-velocity structure. Note that the numbers of P-wave observations at stations NA080, NA090, NA330, and NA400
are smaller (less than 100 picks), so they may have poor resolution and show a different pattern from the surrounding
stations.
Uncertainty of the Phase Detection
We have evaluated the error of the automatic phase detection method by using the waveforms of the 31 welllocated earthquakes on July 1, 2015. Figure S2 shows the
waveforms for radial and vertical components as a functionof the epicenter distance. They show a good alignment for
both P- and S-phases, so these events are not noises, andthe locations are mostly correct.
We manually pick the P-wave arrival of the 584 waveforms and the S-wave arrival of the 337 waveforms. The
histogram of the residuals between the manual and automatic detections is shown in the Figures 6a and 6b. The
rms (root mean square) error of the P-phases is 0.14 s. The
waveforms at the distant stations have multiple phases ataround the S-wave arrival, so it is difficult to detect the correct S-wave arrival time. Therefore, average rms error of
the S-phases is as large as 0.77 s. If we remove the stations
at the epicenter distance larger than 80 km, the rms error
decreases to 0.35 s. The error of the automatic method is
considered acceptable compared to the uncertainty of thevelocity structure and local site conditions, which is on theorder of a second.
Figure 6c shows the error between the observed and estimated arrival times after the relocation, as a function of
the epicenter distance. The rms error for the P-phase is 0.38
s, and that for the S-phase is 1.03 s. The P-phase error isabout constant for any distance, but S-phase error increases
as a function of the distance. As shown in Figure S2, theS-wave arrivals at the distant station have an emergent onset, which makes it difficult to detect S-phases accurately.
Note that the phase arrival times at the stations close to
the epicenter are more sensitive to the depth of the event.
For this dataset, the rms error of the S-phase at the stations
within 30 km from the epicenter is only 0.45. Therefore, we
think our phase detection has reasonable accuracy.
Uncertainty of the Aftershock Catalog
Figure 7 shows histograms of the rms, the hypocenter
location error in the X (longitude), Y (latitude), and Z(depth) directions before the relocation by VELEST. The
average rms is 0.66 s, and the depth error is larger thanthe horizontal error. The error in the X direction is slightly
smaller than that in the Y direction, due to the network
geometry. In order to evaluate the effect of the number ofphase detections, we plot the relationship between the number of phase detections and the above uncertainty parameters in Figure 8. The rms residuals increase as a functionof the number of phase detections, whereas the uncertaintyof the location decreases as the number of phase detections
increases. As the number of phase detections increases, ingeneral, the distance from the hypocenter to the station
increases. The uncertainty of the velocity structure accumulates for distant stations, which results in the larger rms.
If we select earthquakes determined by 20 picks or more,
we are able to filter out most of the earthquakes with horizontal errors larger than 3.5 km. In order to remove earthquakes far from the network, we used earthquakes with anazimuthal gap of stations less than 240 degree. By applying this condition on the network geometry, we are able
to filter out the earthquakes with relatively large errors in

X - 4 YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE
Figure 8 (gray circles above the line of y=20). Therefore,
we define these earthquakes, i.e. earthquakes with 20 ormore phase detections and azimuthal gap of the stations
less than 240 degree, as well-located earthquakes. Figures
2 and 3 show the well-located earthquakes with larger symbols. They largely represent the distribution of the wholeaftershock catalog.
The seismic network that we use in this study covers thefault rupture area with about 20 km spacing. With thisstation spacing, it may be difficult to resolve the shallow
depth earthquakes. Waveforms with both P-wave and Swave phase detections contribute for a good depth control.
In our dataset, about 60% of the phase detected waveforms
have both P-wave and S-wave observations. Figure S3 shows
the theoretical arrival times for different source depth. The
arrival times at the closer distance shows a larger difference,and the apparent velocity is faster for the deeper earthquake. Suppose we have a station at the epicenter distance
10 km, the difference of the S-P arrival time for the depth
5 km and 15 km is 0.8 s. We think this difference can be
resolved, even if we have an uncertainty in the automatic
phase detection.
Change of Earthquake Distribution by the Relocation
We have performed a two-step hypocenter relocation using the VELEST program. Figure 9 shows the histogramsof the rms error before and after the relocation. By adding
the station correction factor and optimizing the 1D velocity
structure, the average rms error of the well-located earthquakes decreased from 0.66 to 0.46 s. The average rms errorof the whole catalog also was reduced from 0.44 to 0.31 s.
Figure 10 shows the earthquake locations before and after the relocation. The horizontal differences are very small
inside the seismic network, whereas the earthquakes outside of the network move away from the network by a few
tens of km after the relocation. This is because the inverted
velocity model was faster than the initial model.
Figure 11 shows the depth of the earthquakes before andafter the relocation. In order to clarify the effect of the velocity structure and station correction, the catalog with theinitial velocity model and station correction, and the catalog
with the optimal velocity model and no station correction,
are also shown in the Figure.
The optimal velocity structure is faster than the initial
velocity structure. Therefore, to make the apparent velocity slower (i.e., make the theoretical arrival times later), thedepths of the earthquakes tend to be shallow (see Figures
11a and 11c).
The station corrections are positive at the stations in
Terai, where the epicenter distance is large in general.
Therefore, if we consider the station correction, the theoretical arrival times at the distant station become later,
and the apparent velocity becomes slower. This is similar
to change the velocity model slower. Therefore, the depth ofthe earthquakes becomes deeper (see Figures 11a and 11b).
If we consider both the optimal velocity model and station correction, the depth of the well-located earthquakes
does not change so much (see Figures 11a and 11d). However, the distribution of other earthquakes, shown as graycircles in Figure 11, is scattered in the shallow depth before
the relocation, and became more confined after.

Discussion
Comparison with the previous catalog
There is a similar pattern for the distribution of aftershocks in the previous aftershock studies [e.g. Adhikari
et al., 2015; Baillard et al., 2017], but our study has greater
spatial resolution due to the data from the dense seismic
network made available to us. Our catalog has smaller magnitude of completeness than that in the previous study andmore aftershocks than those of the same period [e.g. Baillard et al., 2017], as shown in Figure 2a. Our catalog andBaillard et al. [2017] have an overlapped period for about
one month. Figure S4 shows the 669 earthquakes with thedifference of the origin time less than 3 s and the differenceof the horizontal location less than 20 km. The horizontal
differences are very small inside the seismic network, andlarger at the border between China and Nepal. In the section profile, our catalog shows more tightly clustered depth
estimates than Baillard et al. [2017]. The depth change inthe dip direction of the MHT was not clearly observed in
the previous study.
The spatial resolution was improved owing to the joint
hypocenter determination. The average rms error decreased
from 0.66 s to 0.46 s, and the estimated station corrections
were consistent with the subsurface geological structure. Although our catalog lacks the first two months of aftershocks,
it is a well-resolved catalog resulting from the many nearsource seismic stations in the fault rupture area.
Aftershocks and 3D structure model
In order to examine the relationship between the aftershocks and the 3D fault structure, a 3D MHT structure
model [Hubbard et al., 2016] has been included in Figures
2b and 3. This model was developed from the geologicalcross section of central Nepal and extended laterally using
the MFT, MBT, and GPA surface traces as a proxy [Hubbard et al., 2016]. There are two ramp structures parallelto the MBT and GPA near Kathmandu, and the middle
décollement is also bounded by two ramps.
Comparison of the cross sections and the aftershock
data shows that most aftershocks are confined between the
ground surface and the MHT shear zone. The aftershock
depths not only match the slip surface, but are also scattered throughout the entire hanging wall. The bottom of theaftershock distribution is consistent with the MHT model.
Figure 12 shows the aftershock distribution and mainshock source model [Kobayashi et al., 2016]. According to
Hubbard et al. [2016], the location of the large slip is consistent with the location of the middle décollement between
the two ramp structures. Our aftershock distribution image
shows that there are fewer events in the large slip area andmore events on the two ramp structures surrounding thelarge slip area. This is consistent with the observation ofother large earthquakes [e.g. Das and Henry, 2003; Shearer
et al., 2003; Yukutake and Iio, 2017]. The western aftershocks were much more densely distributed than the easternaftershocks, and the very dense area seems to be parallel to
the GPA. Interestingly, the large clusters of earthquakes at
the northern side of the GPA is located above the deeper
ramp in the model, and the small cluster south of the GPA
is located above the shallower ramp.
Comparison with the previous seismicity before theGorkha earthquake shows that the large clusters at thenorthern side of the GPA exist before the Gorkha earthquake [Ader et al., 2012; Stevens and Avouac, 2015]. The
locations are consistent with areas with a large gradient ofcoupling (the downdip edge of the locked MHT fault zone)
where stress buildup is at a maximum [Stevens and Avouac,
2015; Avouac et al., 2015]. There is almost no seismicity
south of the GPA prior to the Gorkha earthquake [Aderet al., 2012; Stevens and Avouac, 2015]. The Gorkha earthquake ruptured the lower edge (northern end) of the locked
MHT, which, we assume, has activated the southern clusters. We propose that these smaller clusters south of the
GPA are generated by the stress heterogeneities resulting
from the shallower ramp structure.

Conclusions
In this paper, we created a well-resolved aftershock catalog covering the time period up to 11 months after the2015 Gorkha earthquake by processing data from the NAMASTE temporary seismic network that was deployed in a

YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE
near-source region. We identified about 15000 events, withthe magnitude of completeness of M2.
Our catalog shows two major clusters north of the GPA
and smaller clusters in the south. The southern clusters areshallower than those of the northern clusters, which is consistent with the dip direction of the MHT. Most aftershocks
are confined between the ground surface and the MHT shear
zone, and their distribution may reflect structural complexity along the MHT.
Our aftershock distribution shows that there are fewer
events in the large slip area, and the clusters occurred on
the two ramp structures in the 3D MHT model surrounding the large slip area. This may suggest a larger stressaccumulation on the ramp structures. Compared with theprevious seismicity before the Gorkha earthquake, there isalmost no seismicity south of GPA. The Gorkha earthquake
ruptures the lower edge of the locked MHT, which may haveactivated the southern clusters.
We obtain a 1D velocity structure and station correction factors using a joint hypocenter determination. The
P-wave velocity at depths of <20 km was almost constantat 6.0 km/s, which is similar to Vp of the shallow continental crust.

Data and Resources
We use the seismic waveform data obtained by the groupof the Department of Mining and Geology in Nepal, the University of California at Riverside, the University of Texas at
El Paso, Stanford University, and Oregon State University.
The seismic data is available on the IRIS Data Management Center website at https://www.fdsn.org/networks/
detail/XQ_2015/. (last accessed July 2019). Some plots
were made using the Generic Mapping Tools version 4.5.7[Wessel and Smith, 1991]. The supplemental material contains four figures (Figures S1 to S4) and the aftershock catalog (Table S1).
Acknowledgments. We thank the entire NAMASTE teamincluding students, staffs, and scientists from the Departmentof Mining and Geology in Nepal, the University of California
at Riverside, the University of Texas at El Paso, Stanford University, and Oregon State University for operating the temporary
seismic network in Nepal. We appreciate valuable comments fromDr. Jim Mori in Kyoto University.

References
Ader, T., Avouac, J., Liu-Zeng, J., Lyon-Caen, H., Bollinger,
L., Galetzka, J., Genrich, J., Thomas, M., Chanard, K., andSapkota, S. N. (2012). Convergence rate across the Nepal
Himalaya and interseismic coupling on the Main Himalayan
Thrust: Implications for seismic hazard. Journal of Geophysical Research: Solid Earth, 117(B4).
Adhikari, L., Gautam, U., Koirala, B., Bhattarai, M., Kandel,T., Gupta, R., Timsina, C., Maharjan, N., Maharjan, K., andDahal, T. (2015). The aftershock sequence of the 2015 April
25 Gorkha-Nepal earthquake. Geophysical Journal International, 203(3):2119–2124.
Avouac, J.-P., Meng, L., Wei, S., Wang, T., and Ampuero, J.-P.
(2015). Lower edge of locked Main Himalayan Thrust unzipped by the 2015 Gorkha earthquake. Nature Geoscience,
8(9):708.
Bai, L., Liu, H., Ritsema, J., Mori, J., Zhang, T., Ishikawa, Y.,
and Li, G. (2016). Faulting structure above the Main Himalayan Thrust as shown by relocated aftershocks of the 2015
Mw7.8 Gorkha, Nepal, earthquake. Geophysical Research Letters, 43(2):637–642.
Bai, L., Klemperer, S., Mori, J., Karplus, M., Ding, L., Liu,
H., Li, G., Song, B., and Dhakal, S. Lateral variation of
the Main Himalayan Thrust controls the rupture length of
the 2015 Gorkha earthquake in Nepal. Science Advances,5(6):eaav0723.
Baillard, C., Lyon-Caen, H., Bollinger, L., Rietbrock, A., Letort,J., and Adhikari, L. B. (2017). Automatic analysis of the
Gorkha earthquake aftershock sequence: evidences of structurally segmented seismicity. Geophysical Journal International, 209(2):1111–1125.

X-5

Christensen, N. I. and Mooney, W. D. (1995). Seismic velocity structure and composition of the continental crust: A
global view. Journal of Geophysical Research: Solid Earth,
100(B6):9761–9788.
Dahal, R. (2006). Geology for technical students. Bhrikuti Academic Publication, Kathmandu, 756.
Das, S. and Henry, C. (2003). Spatial relation between mainearthquake slip and its aftershock distribution. Reviews ofgeophysics, 41(3).
Galetzka, J., Melgar, D., Genrich, J. F., Geng, J., Owen, S., Lindsey, E. O., Xu, X., Bock, Y., Avouac, J.-P., and Adhikari,
L. B. (2015). Slip pulse and resonance of the Kathmandu
basin during the 2015 Gorkha earthquake, Nepal. Science,
349(6252):1091–1095.
Ghosh, A., Mendoza, M., LI, B., Karplus, M. S., Nabelek, J.,
Sapkota, S. N., Adhikari, L. B., Klemperer, S. L., and Velasco,
A. A. (2017). Fault structure in the Nepal Himalaya as illuminated by aftershocks of the 2015 Mw 7.8 Gorkha earthquake
recorded by the local Namaste network. AGU Fall Meeting
Abstracts.
Hildyard, M. W., Nippress, S. E., and Rietbrock, A. (2008).
Event detection and phase picking using a time-domain estimate of predominate period T pd . Bulletin of the SeismologicalSociety of America, 98(6):3025–3032.
Hildyard, M. W. and Rietbrock, A. (2010). T pd , a damped predominant period function with improvements for magnitude
estimation. Bulletin of the Seismological Society of America,
100(2):684–698.
Hubbard, J., Almeida, R., Foster, A., Sapkota, S. N., Bürgi,
P., and Tapponnier, P. (2016). Structural segmentation controlled the 2015 Mw 7.8 Gorkha earthquake rupture in Nepal.
Geology, 44(8):639–642.
Karplus, M., Sapkota, S., Nabelek, J., Pant, M., Velasco, A.,
Patlan, E., Klemperer, S., Braunmiller, J., Ghosh, A., andSapkota, B. (2015). Aftershocks of the M7.8 Gorkha (Nepal)
Earthquake: Early Results from Project NAMASTE. AGU
Fall Meeting Abstracts.
Karplus, M. S., Pant, M., Velasco, A. A., Nabelek, J., Kuna,
V. M., Sapkota, S. N., Ghosh, A., Mendoza, M., Adhikari,
L. B., and Klemperer, S. L. (2017). Structure and tectonics ofthe Main Himalayan Thrust and associated faults from recent
earthquake and seismic imaging studies using the NAMASTE
array. AGU Fall Meeting Abstracts.
Kissling, E., Ellsworth, W., Eberhart-Phillips, D., and Kradolfer,
U. (1994). Initial reference models in local earthquake tomography. Journal of Geophysical Research: Solid Earth,
99(B10):19635–19646.
Kissling, E., Kradolfer, U., and Maurer, H. (1995). Program
VELEST user’s guide - short introduction. Institute of Geophysics, ETH Zurich.
Kobayashi, H., Koketsu, K., Miyake, H., Takai, N., Shigefuji,M., Bhattarai, M., and Sapkota, S. N. (2016). Joint inversion of teleseismic, geodetic, and near-field waveform datasetsfor rupture process of the 2015 Gorkha, Nepal, earthquake.
Earth, Planets and Space, 68(1):66.
Letort, J., Bollinger, L., Lyon-Caen, H., Guilhem, A., Cano, Y.,
Baillard, C., and Adhikari, L. B. (2016). Teleseismic depth estimation of the 2015 Gorkha-Nepal aftershocks. GeophysicalJournal International, 207(3):1584–1595.
McNamara, D. E., Yeck, W., Barnhart, W. D., Schulte-Pelkum,
V., Bergman, E., Adhikari, L., Dixit, A., Hough, S., Benz,
H. M., and Earle, P. (2017). Source modeling of the 2015
Mw 7.8 Nepal (Gorkha) earthquake sequence: Implications
for geodynamics and earthquake hazards. Tectonophysics,
714:21–30.
Mendoza, M., Ghosh, A., Karplus, M., Nabelek, J., Sapkota, S.,
Adhikari, L., Klemperer, S., and Velasco, A. (2016). Alongstrike variations in the Himalayas illuminated by the aftershock sequence of the 2015 Mw 7.8 Gorkha earthquake using
the Namaste local seismic network. AGU Fall Meeting Abstracts.
Mendoza, M., Ghosh, A., Karplus, M., Klemperer, S., Sapkota, S., Adhikari, L., and Velasco, A. (2016). Duplex in
the Main Himalayan Thrust illuminated by aftershocks ofthe 2015 Mw 7.8 Gorkha earthquake. Nature Geoscience,
doi:10.1038/s41561-019-0474-8 (in press).
Monsalve, G., Sheehan, A., Schulte-Pelkum, V., Rajaure, S.,
Pandey, M., and Wu, F. (2006).
Seismicity and onedimensional velocity structure of the Himalayan collision zone:
Earthquakes in the crust and upper mantle. Journal of Geophysical Research: Solid Earth, 111(B10).

X - 6 YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKENábĕlek, J., Hetényi, G., Vergne, J., Sapkota, S., Kafle, B., Jiang,
M., Su, H., Chen, J., and Huang, B.-S. (2009). Underplating in
the Himalaya-Tibet collision zone revealed by the Hi-CLIMB
experiment. Science, 325(5946):1371–1374.
Negi, S. S., Paul, A., Cesca, S., Kriegerowski, M., Mahesh, P.,
and Gupta, S. (2017). Crustal velocity structure and earthquake processes of Garhwal-Kumaun Himalaya: Constraints
from regional waveform inversion and array beam modeling.
Tectonophysics, 712:45–63.
Pandey, M., Tandukar, R., Avouac, J., Lave, J., and Massot, J. (1995). Interseismic strain accumulation on the Himalayan crustal ramp (Nepal). Geophysical Research Letters,
22(7):751–754.
Pandey, M., Tandukar, R., Avouac, J., Vergne, J., and Heritier, T. (1999). Seismotectonics of the Nepal Himalaya froma local seismic network. Journal of Asian Earth Sciences,17(5-6):703–712.
Pant, M., Velasco, A., Karplus, M., Patlan, E., Ghosh, A.,
Nabelek, J., Kuna, V., Sapkota, S., Adhikari, L., and Klemperer, S. (2016). Aftershock stress analysis of the April 2015
Mw 7.8 Gorkha earthquake from the NAMASTE project.
AGU Fall Meeting Abstracts.
Shearer, P. M., Hardebeck, J. L., Astiz, L., and Richards-Dinger,
K. B. (2003). Analysis of similar event clusters in aftershocksof the 1994 Northridge, California, earthquake. Journal ofGeophysical Research: Solid Earth, 108(B1).
Stevens, V. and Avouac, J. (2015). Interseismic coupling on
the Main Himalayan Thrust. Geophysical Research Letters,
42(14):5828–5837.
Takanami, T. and Kitagawa, G. (1988). A new efficient procedure
for the estimation of onset times of seismic waves. Journal of
Physics of the Earth, 36(6):267–290.
Takanami, T. and Kitagawa, G. (1991). Estimation of the arrival times of seismic waves by multivariate time series model.

Annals of the Institute of Statistical mathematics, 43(3):407–
433.
Tamaribuchi, K. (2018). Evaluation of automatic hypocenter determination in the JMA unified catalog. Earth, Planets andSpace, 70:141.
Tilmann, F., Ni, J., and INDEPTH III, S. T. (2003). Seismic
imaging of the downwelling Indian lithosphere beneath central Tibet. Science, 300(5624):1424–1427.
Ueno, H., Hatakeyama, S., Aketagawa, T., Funasaki, J., andHamada, N. (2002). Improvement of hypocenter determination procedures in the Japan meteorological agency. Quarterly
Journal of Seismology, 65:123–134.
Wang, X., Wei, S., and Wu, W. (2017). Double-ramp on the Main
Himalayan Thrust revealed by broadband waveform modelingof the 2015 Gorkha earthquake sequence. Earth and Planetary
Science Letters, 473:83–93.
Wessel, P. and Smith, W. (1991). Free software helps map anddisplay data. Eos, 72(441):445–446.
Wu, S., Yamada, M., Tamaribuchi, K., and Beck, J. (2014).
Multi-events earthquake early warning algorithm using a
Bayesian approach.
Geophysical Journal International,
200(2):791–808.
Yamada, M. 2017). Prototype of P-wave Detection Algorithm for
Earthquake Early Warning. Seismological Society of Japanannual meeting Abstracts.
Yukutake, Y. and Iio, Y. (2017). Why do aftershocks occur?
Relationship between mainshock rupture and aftershock sequence based on highly resolved hypocenter and focal mechanism distributions. Earth, Planets and Space, 69(1):68.
M. Yamada, Kyoto University, Uji, Gokasho, 611-0011, Japan(masumi@eqh.dpri.kyoto-u.ac.jp)

Figures and Tables
Table 1. Initial and optimal models for the P-wave and S-wave velocity structure in km/s.
Depth (km) Vp(ini) Vp(opt) Vs(ini) Vs(opt)
0-2
5.56
6.01
3.18
3.56
2-4
5.56
6.01
3.18
3.56
4-6
5.56
6.01
3.18
3.56
6-8
5.56
6.01
3.18
3.56
8-10
5.56
6.01
3.18
3.56
10-15
5.56
6.01
3.18
3.56
15-20
5.56
6.01
3.18
3.56
20-23
5.56
6.01
3.18
3.56
23-25
6.70
6.69
3.71
3.79
25-30
6.70
6.75
3.71
3.82
30-35
6.70
6.80
3.71
3.83
35-40
6.70
6.92
3.71
3.96
40-45
6.70
6.93
3.71
4.00
45-50
6.70
7.1
3.71
4.12
50-55
6.70
7.22
3.71
4.15
55-60
8.10
8.23
4.63
4.67
608.10
8.23
4.63
4.69

YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE

(a)
80°

84°

88°
Terai

30°

Sub-Himalaya

GPA

Lesser-Himalaya
Higher-Himalaya
Tibetan-Tethys
Paleozoic Granite
Tertiary Leucogranite

28°

MFT: Main Frontal Thrust
MBT: Main Boundary Thrust
MCT: Main Central Fault
GPA: Gorkha-Pokhara Anticlinorium
KTM: Kathmandu
100km
26°

(b)
Terai

S

Sub-Himalaya

Lesser-Himalaya

Higher-Himalaya

Tibetan-Tethys
TibetanPlate

Indian ShieldN
Figure 1. Geology in Nepal modified after Dahal [2006].
(a) Geological map and aftershock area shown in shaded
region. Large and small stars show the location of themainshock and the largest aftershock on May 12, respectively. (b) Schematic NS cross section.

X-7

X - 8 YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE

Figure 2.

YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE

X-9

Figure 2: Spatiotemporal variation of our catalog. Larger symbol represents the well-located earthquakes. Largeand small stars show the location of the mainshock and the largest aftershock on May 12, respectively. (a) Timesequence of the aftershocks. The numbers on the vertical axis represent the year (15 or 16) and the month. The
shaded period lacks some data (less than 25 stations). The catalog of Baillard et al. [2017] is used before June25, 2015. (b) Map view of the aftershock distribution. Major fault systems and the 3D MHT model are added[Hubbard et al., 2016]. The solid triangles show the locations of the temporary seismic stations. Thick graylines show the boundary of the countries. The top right figure shows the map of Nepal and the location of themainshock.

Figure 3. Profiles of the aftershocks. (a) Dashed lines
show the locations of cross sections. Other symbols arein the same format as Figure 2b. (b)-(d) NS and EW sections of the aftershocks. Circles with an errorbar show
the well-located earthquakes and lighter circles show all
the other earthquakes. The solid and open inverted triangles show the locations of MBT and MCT, respectively.
The 3D MHT model [Hubbard et al., 2016] is represented
by the dashed line.

Figure 4. Histogram of the (a) magnitude, (b) number of picks used for location determination, and (c) azimuthal gap of the stations.

X - 10 YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE

Figure 5. Results of joint hypocenter determination.
(a) Initial and optimal 1D velocity models for Vp and Vs.
Gray bars show the depth distribution of all earthquakes
after the relocation. (b) Station correction factors at theaftershock stations. The circles and triangles show thepositive and negative delay time in seconds, respectively.
The background color shows the altitude. The stars markthe locations of the mainshock and the largest aftershockon May 12.

(a)

80

(b)

6

rms=0.14s

rms=0.77s

60

CountsCounts150

abs(Tobs -Test) (s)

200

100
50

40
20

0
-1

0

1

0
-1

Tpauto-Tpmanual (s)

0

1

(c)
Tp rms=0.38s
Ts rms=1.03s

4
2
0

0

100

Tsauto-Tsmanual (s)200Distance (km)Figure 6. Histograms of the residuals for (a) automatic and manual P-wave arrival times, (b) automaticand manual S-wave arrival times, and (c) observed andestimated arrival times for P- and S-phases as a functionof the epicenter distance.

Counts

1500
1000
500
0

0

0.5

RMS (s)

1

4000

4000

3000

3000

2000

2000

1000

1000

0

0

5

ErrorX (km)

10

0

1500

all data
filtered

1000
500

0

5

10

ErrorY (km)Figure 7. Histograms of the rms, the hypocenter location error in the X (longitude), Y (latitude), and Z(depth) directions before the relocation by the VELEST.
Dark color shows the well-located earthquakes and graycolor shows all earthquakes.

0

0

5

ErrorZ (km)10

YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE X - 11

80

80

80

60

60

60

60

40

40

40

40

20

20

20

20

0

0

0.5

1

RMS (s)005ErrorX (km)

10

0

0

5

10

0

all data
filtered0ErrorY (km)5

ErrorZ (km)Figure 8. Relationship between the number of P- andS-phase detections vs the rms, the hypocenter location
error in the X (longitude), Y (latitude), and Z (depth)
directions before the relocation by the VELEST. Blackcircles show the well-located earthquakes and gray circles
show all the other earthquakes.

2000
before JHD
filtered
after JHD
filtered

1500

CountsNo. of picks

80

1000

500

0

0

0.5

1

RMS (s)Figure 9. Histograms of the rms error before (gray) andafter (black) the joint hypocenter determination. The
solid and dashed lines show the histograms for the whole
and well-located catalogs, respectively. Cross symbols
show the average rms errors for the well-located catalog.

1.5

10

X - 12 YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE

Figure 10. Comparison of the well-located earthquake
locations before and after the relocation shown in small
circles. The bars on the map frame show the range (20km) of earthquakes included in the cross sections.The
solid triangles show the locations of the temporary seismic stations. Large and small stars show the locationof the mainshock and the largest aftershock on May 12,
respectively.

YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE X - 13

Figure 11. Comparison of the earthquake locations with(a) initial velocity model and no station correction (before the relocation) (b) initial velocity model and station correction (c) optimal velocity model and no station
correction, and (d) optimal velocity model and station
correction. The format is the same as the Figure 3d.

X - 14 YAMADA ET AL.: A REFINED AFTERSHOCK CATALOG FOR THE 2015 GORKHA EARTHQUAKE

Figure 12. Aftershock locations and slip model of themainshock [Kobayashi et al., 2016]. Gray contour shows
the slip in meters, and the small circles show the frequency distribution of the aftershocks. It is the numberof hypocenters within a radius of 3 km of each earthquake. Major fault systems and the 3D MHT model (ingray contours) are added. Large and small stars show
the location of the mainshock and the largest aftershockon May 12, respectively.