Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk

TECTO-127121; No of Pages 14 Tectonophysics xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

Tectonophysics journal homepage: www.elsevier.com/locate/tecto

Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk Chengli Liu a, Yong Zheng a,⁎, Rongjiang Wang b, Bin Shan a, Zujun Xie a, Xiong Xiong a, Can Ge c a b c

State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China GFZ German Research Centre for Geosciences, D-14473 Potsdam, Germany School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China

a r t i c l e

i n f o

Article history: Received 23 February 2016 Received in revised form 26 April 2016 Accepted 19 May 2016 Available online xxxx Keywords: Gorkha earthquake Gedetic datasets Rupture process Seismic risk

a b s t r a c t The rupture processes of the 2015 April 25 Gorkha earthquake and its strongest aftershock occurred on May 12 in Nepal are investigated by joint inversion of seismological and geodetic data. Synthetic test shows that the sedimentary layers in the source region play an important role in the rupture process inversion. Our optimized model of the mainshock shows that the rupture has a unilateral propagation pattern. The dominant mechanism is pure thrust with maximum slip of 5.8 m, the rupture scale extends ~60 km along dip and ~150 km along strike, and the largest static stress change is ~7.6 MPa. The total seismic moment is 7.87 × 1020 N m, equivalent to Mw 7.9. Most seismic moment was released within 80 s and the majority seismic moment was released at the first 40 s. The rupture propagated in main slip asperity with a velocity of ~3.0 km/s. The strong aftershock magnitude is about Mw 7.3, and the peak slip is about 5.0 m, close to the peak slip of the mainshock. Moreover, the slips of the mainshock and the aftershocks are in good complementary, suggesting a triggering relationship between them. Considering the strain accumulation, the Gorkha earthquake ruptured only part of the seismic gap alone, thus still poses high earthquake risk, especially in the west side of the mainshock rupture zone. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The Indian Plate is subducting northward beneath the Eurasia Plate with a convergence rate of ~ 3.6 cm/yr (Altamimi, 2009). Due to this rapid convergence, the Himalayan arc is prone to strong earthquakes frequently (Bilham, 2004; Lavé et al., 2005; Kumar et al., 2006; Bettinelli et al., 2006). The central Nepal is located at the central part of the Himalayan collision zone, an active tectonic zone. Based on comparative studies of historical seismic records and the convergent rate, the central Nepal has accumulated significant tectonic strain; and it is highly possible that large earthquakes with estimated seismic magnitude larger than Mw 7.7 will occur there in the future (Feldl and Bilham, 2006; Bilham et al., 2001; Mukhopadhyay et al., 2011). On 25 April 2015, a devastating earthquake struck the Gorkha region in the central Nepal (hereinafter as the Gorkha earthquake), resulting in numerous casualties and properties loss. The epicenter (28.147°N, 84.708°E) reported by USGS (http://earthquake.usgs.gov/earthquakes/ eventpage/us20002926#general_summary) was ~77 km northwest of the Nepalese capital of Kathmandu, with a hypocenter depth of 15 km. The global Centroid-Moment Tensor (GCMT) solution of this event (http://www.globalcmt.org/CMTsearch.html) indicates that it is ⁎ Corresponding author. E-mail address: [email protected] (Y. Zheng).

a nearly pure thrust event with fault geometry of strike 293°, dip 7°, and rake 108°. Its centroid depth is ~12 km and the centroid epicenter is to the southeast of the hypocenter (27.7°N, 85.37°E), with a seismic moment of 7.76 × 1020 N m, corresponding to a moment magnitude of Mw 7.9. The Gorkha earthquake was followed by a substantial aftershock sequence, among which, the largest aftershock occurred on May 12, 2015 and had a magnitude of Mw 7.3 (07:05:19 UTC, 27.819°N, 86.080°E, and the hypocenter depth 15 km) to ~ 150 km east of the mainshock (Fig. 1). An accurate rupture model is of fundamental importance for further studies on source properties and seismic hazards, such as computing Coulomb stress transfer, characterizing the seismic cycle and hazard assessment related to the Himalayan thrust fault system. Moreover, following the mainshock, many aftershocks occurred; some of the aftershocks were quite big and also caused severe damages to the source region. Considering the active tectonic environment and the seismicity in the Nepal region, studying on the source parameters and the relationship between the mainshock and the strong aftershocks may help understand the seismogenic environment and evaluate the potential seismic hazard in the Himalayan fault zone. Lots of rupture models have been derived for this event by using teleseismic data (Zhang et al., 2015; Wang et al., 2015; Fan and Shearer, 2015; Yagi and Okuwaki, 2015), geodetic data (Wang and Fialko, 2015; Lindsey et al., 2015; Diao et al., 2015; Feng et al., 2015),

http://dx.doi.org/10.1016/j.tecto.2016.05.034 0040-1951/© 2016 Elsevier B.V. All rights reserved.

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

2

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

78˚E 33˚N

84˚E

81˚E

87˚E

90˚E

93˚E 33˚N

Tibet Plateau 2015/04/25_Mw7.9 30˚N

Nepal 15

05

India

30˚N Un rup tur ed

2015/05/12_Mw7.3

1

Kathmandu 833 19

27˚N

34

27˚N

3.6 cm/yr

km 24˚N

24˚N 0

78˚E

81˚E

84˚E

87˚E

90˚E

100 200 93˚E

Fig. 1. Tectonic setting of the 2015 Mw 7.9 Gorkha earthquake and its aftershocks distribution. The red barbed line indicates the Himalayan Thrust fault. The black arrow shows Indian plate motion relative to Eurasia computed using the rotation poles of Eurasian plate in ITRF 2005 from Altamimi (2009). Black ellipses show estimated locations and possible rupture areas of historical earthquakes (Ader et al., 2012). Red stars show the relocated epicenter of the mainshock and its largest aftershock, and its aftershocks from April 25, 2015 to June 8, 2015 indicated by yellow circles (Adhikari et al., 2015), which are shown scaled with symbol size proportional to seismic magnitude. The focal mechanisms are GCMT best double-couple solutions for the mainshock and the May 12, 2015 Mw 7.3 aftershock. The green square is the location of the Nepalese capital of Kathmandu. The black rectangle is the fault plane used in this study. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

or joint inversion with multiple datasets (Galetzka et al., 2015; Zhang et al., 2015; Grandin et al., 2015; Avouac et al., 2015). Most of models indicate that the Gorkha earthquake is a unilateral rupture event. Finite fault rupture models determined by teleseismic indicate that the largest rupture slips are concentrated in the area to the east of hypocenter with magnitude of 5–7.5 m (http://www.earthobservatory.sg/news/april25-2015-nepal-earthquake; Wang et al., 2015; Yagi and Okuwaki, 2015). However, due to the well-known propagation effects, teleseismic body waves are insufficient for resolving rupture processes in detail. On the other hand, geodetic data (e.g., InSAR data and static GPS coseismic offsets) can provide near field constraints, and is more sensitive to the detailed slip pattern of the rupture model. The inversions of geodetic data show that the large slips of the Gorkha earthquake tend to concentrate in the central region of the rupture zone, and the largest rupture is also located east to the hypocenter (Lindsey et al., 2015; Diao et al., 2015; Feng et al., 2015; Wang et al., 2015). The distributions of along strike slip in these models are generally similar, with unilateral rupture expanding ~150 km along strike. There are remarkable differences in the total slip and the distribution of along dip slip among the geodesy-based models. Most of these models do not have significant surface ruptures; however, Galetzka et al. (2015) and Lindsey et al. (2015) found 1–2 m slip near the surface at the updip part of the southern slip patch. Avouac et al. (2015) investigated the coseismic slip pattern by inverting teleseismic and InSAR datasets; their rupture model is consistent with most of the geodetic based models, with no significant slip near the surface. Grandin et al. (2015) developed a coseismic rupture model using teleseismic waves, strong motion data, high-rate GPS, static GPS, and synthetic aperture radar (SAR) data. The rupture pattern appears simple in their slip model, largely as a result of low-frequency waveform data used in their studies or the smoothing factor used in their inversion. These differences may be mainly due to nonuniform spatial sampling, different data types and difference between InSAR data from different satellites. Moreover, an accurate location of the hypocenter is important for the inversion of rupture process. For most large earthquakes, near-field or regional recordings are not available timely. So the current proposed rupture models almost adopt the location of U.S. Geological Survey National Earthquake Information Center (NEIC),

and its error usually is as large as 25 km (Engdahl et al., 1998), which might be another reason for the model differences. Here, we jointly used teleseismic waves, strong motion data, high-rate GPS, static GPS displacement, and Interferometric Synthetic Aperture Radar (InSAR) data, and take into account the influence of the thick sediments in Kathmandu basin to seismic waveforms, to obtain a refined model of the spatiotemporal history of slip for the Gorkha earthquake and its strongest aftershock. Based on our rupture models and the distribution of the aftershocks, we then discuss the potential seismic risk in the Himalayan Thrust fault after the occurrence of the Gorkha earthquake. 2. Data and methods 2.1. Geodetic data Near-field observations, especially near-field static displacements, are particularly useful for constraining the slip distribution of large complex ruptures. The Gorkha earthquake is the first occurrence of a large continental thrust earthquake which is well recorded by a group of high-rate GPS stations closely encompassing the rupture area (Fig. 2a). In addition, the rupture area of the Gorkha earthquake is also well covered by InSAR line-of-sight displacement data from ALOS-2 observations, with both the ascending and descending paths (Fig. 2b and c), which can provide well-constrained slip distribution. The combination of these measurements provides an extraordinary opportunity to study the kinematics of the source process. In this work, both GPS and InSAR datasets are collected and inverted for the rupture processes of the mainshock together with seismic data. The geodetic dataset consists of 19 coseismic GPS displacements and 5 high-rate GPS from the Tectonics Observatory network in Nepal (Galetzka et al., 2015), and displacement in the Line of Sight direction obtained from ALOS-2 resampled at 4180 points (Lindsey et al., 2015). 2.2. Seismic data We downloaded teleseismic data from the Incorporated Research Institutions for Seismology (IRIS) network within epicentral distances

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

78˚E 30˚N

80˚E

82˚E

(a)

84˚E

3

86˚E

88˚E

90˚E

SMKT

2015/04/25_Mw7.9

JMSM

DNGD NPGJ

80 cm

PYUT

28˚N

2015/05/12_Mw7.3

DNSGGHER BESI DMAU BELT

CHLM SYBC KKN4 KIRT NAST DAMA KATNP SNDL

Strong Motion

TPLJ

RMTE

Static GPS

km

High−rate GPS

26˚N

0

(b)

Path 48 2015/02/22−2015/05/03

(c)

75

150

Path 157 2015/02/21−2015/05/02

29˚N

2015/04/25_Mw7.9

2015/04/25_Mw7.9 2015/05/12_Mw7.3

28˚N

2015/05/12_Mw7.3

km 27˚N

cm −100 −50 84˚E

0

0

40

km 80

50 100 85˚E

cm −100 −50

86˚E

87˚E

84˚E

0

0

40

80

50 100 85˚E

86˚E

87˚E

Fig. 2. (a) Map view of the Mw 7.9 Gorkha earthquake and some datasets used in this study. Red stars and black beachballs denote the relocated epicenter (Adhikari et al., 2015) and GCMT Solution of the mainshock and aftershock. The green triangle and red triangles represent the strong-motion stations and high-rate GPS stations, respectively. The blue squares indicate all the GPS coseismic observations. The upper inset shows the location of the teleseismic stations used in this study. (b) and (c) are the line-of-sight displacement in centimeters from ALOS-2 along descending Path and ascending paths, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

between 30o and 90o. After removing instrument responses, we selected 42 P and 29 SH waveforms with high signal-to-noise ratios and relative evenly azimuthal distribution (Fig. 2a, inset). For the azimuths with dense stations, only those records showing coherent and clear phases in most of the neighboring records were retained. We integrated the broadband seismograms into displacement and then band-pass filtered them at 0.0033–1 Hz band to remove high frequency noise and long period disturbances. Besides teleseismic records, one USGS strong motion station is located close to the causative fault, and provides near field acceleration records (Fig. 2a). We collected the accelerograms from this station and corrected the baseline shifts of the accelerograms following the method of Wang et al. (2011). The corrected acceleration records were then converted into displacements. In order to remove high frequency noise, a low-pass filter with a corner frequency of 1 Hz was applied to filter the displacements records. 2.3. Modeling strategy and fault parameterization In order to obtain the source kinematic models, we adopted the method proposed by Ji et al. (2002; 2003), which performs the waveform inversion in the wavelet domain, along with a simulated annealing method to simultaneously invert for slip amplitude, rake angle, rise

time, and average rupture velocity. The advantage of this approach is a reduction in the number of free parameters, and is capable of combining multiple datasets, including geodetic and seismic data, to retrieve the rupture evolution on the fault, which has been witnessed by numerous previous studies (Delouis et al., 2002; Sladen et al., 2010; Konca et al., 2010; Wei et al., 2011; Shao et al., 2011; Fielding et al., 2013; Yue et al., 2014; Liu et al., 2015). Due to the nature of the geodetic data and the seismic datasets, the geodetic datasets from GPS and InSAR usually provide strong constraints on the spatial distribution of the total coseismic slip, while the seismological datasets are more sensitive to temporal evolution of the earthquake ruptures. Moreover, the sensitivities of the near field geodetic data and seismic data to the spatial distribution are also different; the near field geodetic data usually are more sensitive to the shallower slip, while the seismic data are more sensitive to the deeper slip. So these two types of datasets are complementary. Since the determination of a source rupture model for a given fault geometry is an underdetermined problem due to numerous trade-offs among model parameters, if the geodetic observations of coseismic deformation are available and inverted jointly with the seismological observations, these trade-offs can be significantly reduced. In addition, in order to effectively eliminate singular values between neighboring subfaults effectively and more smooth the visual effect, we imposed additional constraints by minimizing the gradient difference of

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

4

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

neighborhood subfaults using the spatial Laplacian smoothing and by assuming the seismic moment in the Global Centroid Moment Tensor (GCMT) catalog (7.659 × 1020 N m) as the reference moment. Based on the GCMT Solution (http://earthquake.usgs.gov/ earthquakes/eventpage/us20002926#scientific_tensor:us_us_ 20002926_mwc_gcmt), and the geological environment in the source region, we chose a slip model that consists of a single rupture plane with a length of 210 km along the strike and 140 km down dip, and with strike and dip angles of 293° and 7°, respectively (Fig. 1). Its top edge has a depth of 7.23 km. The fault plane is further discretized into 294 subfaults, each having a size of 10 km along strike and 10 km along dip. We adopted the epicenter (84.7647°E, 28.2264°N) determined from seismic signals recorded by the Nepal nationwide seismological network (Adhikari et al., 2015). Based on this epicenter, we further conducted a series of preliminary finite fault inversions using waveform datasets to search for the optimal hypocenter depth. The search result shows that the model best fit the waveform at the depth of 15 km. During the inversion process, the slip value for each subfault was allowed to vary from 0 to 8 m, and the search for the rake angle was performed in the range from 78° to 138°, with a search step of 2°. The rise time was allowed to change from 1.6 to 16 s, with a time step of 1.6 s.

3.1. Static inversion We first performed static inversion using the InSAR and static GPS data. The Green's functions for the geodetic data were calculated using the 1D elastic layered earth structure model of Mahesh et al. (2013) (hereafter called the V_Model 1). The fault parameters are described in Section 3.3. The static slip inversions only estimate the total slip on each subfault of the modeled fault, not the temporal evolution of slip, thus reducing two free parameters per each subfault. Hence, it is more rapid and stable to compute them than to model slip history. In our preferred geodetic-only solution (Fig. 3a) (hereafter as Model I), the total length of the rupture reaches about 160 km along the strike, and 60 km along the dip. The slip distribution reaches the maximum slip of 5.2 m, about 70 km east of the epicenter, and is characterized by pure thrust motion. No slip is observed at the shallow part of the fault, which is consistent with previous static studies (Diao et al., 2015; Feng et al., 2015; Wang et al., 2015). In Model I, the slip generally is confined at a depth of 11–19 km, much deeper than the model of Lindsey et al. (2015). However, the static displacements cannot constrain the slip at deep depths quite well, so differences between our static model

(b) Model II

(Static inversion)

Strike= 293˚ Dip= 7˚

(Waveform inversion)

Strike= 293˚ Dip= 7˚

km 0

40

-160

West

East

15 80

Depth (km)

10 40

20 120

-120

-80

-40

0

100

200

(c) Model III

300

400

500

10

20

40

20 80

20

15

80 20

20

40

120

Slip (cm) 0

40

40

-40

60

-80

Distance along dip (km)

-120

40

Distance along dip (km)

-160

km

Depth (km)

(a) Model I

3. Results and discussion

Slip (cm)

600

0

100

200

300

400

500

600

(d) Moment-rate Function

(Joint inversion)

Strike= 293˚ Dip= 7˚

km -40

0

10

20

60

40

40

40

15

20

80 40 60

20 40

120

20

Moment rate (1019N•m/sec)

-80

Depth (km)

-120

80

Distance along dip (km)

-160

2.5 2.0

Model III

1.5

Model II

1.0 0.5 0.0

Slip (cm) 0

100

200

300

400

500

600

0

20

40

60

80

100

Time (s)

Fig. 3. Comparison of slip distributions and moment rate functions of Models II and III. Slip distributions calculated from static datasets (a), seismic datasets (b), and all available datasets (c). The red star indicates the hypocenter, the color bar shows the scale of the slip amplitude, the white arrows represent the slip directions and the contours outline the rupture propagation time in second. (d) Comparison of moment-rate functions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

from seismic signals recorded by the Nepal nationwide seismological network (Adhikari et al., 2015). The Green's functions used to invert seismic data were calculated using V_Model 1. Previous studies found that an inaccurate velocity structure could strongly bias the inversion results which are determined by seismic waveform data (Graves and Wald, 2001; Wald and Graves 2001), especially in the areas with complex crustal structures. Furthermore, local site conditions may act strong effects on ground motion in the basin which has thick sediment. In general case the thick sediments in basin can significantly amplify seismic waveforms, especially in the basin area. In the source region, there is a Kathmandu basin which contains relatively thick sediment (Fig. 4a; Paudyal et al., 2013). For this reason, we built a test to analyze the effects of the sedimentary basin during the inversion. In this test, we applied two velocity models to do the rupture process inversion: V_Model 1 is a 1D crustal model without sedimentary layers (Fig. 5b). Since there are two stations (accelerometer station KATNP, high-rate GPS station NAST) that are located in the Kathmandu basin, it is necessary to quantitatively evaluate the sedimentary effects in the basin. We built a combined model with sedimentary layers (Shengji Wei, personal communication) on top of the basis of model 1 (hereafter called V_Model 2, Fig. 5c). Then we compared the synthetic seismograms which are calculated based on these two velocity models with the observed seismograms. Fig. 5d shows the three-component waveforms at stations KATNP and NAST with those simulated using V_Model 1. Fig. 5e shows the three-component

and Lindsey et al. (2015) may be caused by the insufficient resolution of the static observations. With a few exceptions for GPS points above the largest slip patch, the fit to both the GPS and InSAR data is excellent in both the near and far field (Figs. S1 and S2). During the static inversion, we tested different smoothing values and chose the largest smoothing factor that allowed a good fit to the geodetic data. To evaluate the resolution of our inversions of the geodetic data, we performed a checkerboard test, which use the same parameterization as in the real data inversion. Comparison of the input checkerboard slip model pattern with the slip model obtained from the inversion of the synthetic data indicates that most of the slip can be recovered in both shape and slip amplitude at a scale of larger than 20 × 20 km pattern (Fig. 4). For the scale of the 20 × 20 km pattern, as expected, the resolution for the deepest subfaults of the whole fault is decreased. Overall, this resolution test demonstrates that the slip distribution can be well resolved, as expected from the very large amount of geodetic datasets. The inclusion of waveform datasets, especially the near field waveform data, in our joint inversions should somewhat improve the resolution of deep slip on the faults. 3.2. Waveform inversion We conducted a kinematic slip model inversion for the mainshock using teleseismic data, strong-motion data and High-rate GPS displacement data together. We started with the epicenter location estimated

Static inversion

Joint inversion

-120

-80

-40

0

-160

-120

-80

-40

0

-160

-120

-80

-40

0

-160

-120

-80

-40

0

-160

-120

-80

-40

0

-160

-120

-80

-40

0

-160

-120

-80

-40

0

-160

-120

-80

-40

0

-160

-120

-80

-40

0

20 × 20 km

-160

40

80

30 × 30 km

120

40

80

120

40 × 40 km

Distance along dip (km)

Distance along dip (km)

Distance along dip (km)

Input model

5

40

80

120 Strike= 295˚ Dip= 7˚

km

Slip (cm) 0

100

200

300

400

500

Fig. 4. Checkerboard tests for different slip scales by static datasets and joint all datasets inversion. The subplots in left column are the input models, the subplots in the middle column are the static inversion, and the subplots in the right column are the jointly inversion. From upper to lower the spatial resolution are 20 km × 20 km, 30 km × 30 km, and 40 km × 40 km, respectively.

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

(a)

(b) 85.2˚E

85.4˚E

85.6˚E

(c)

Velocity (km/s)

0 Vs

10 Depth (km)

KATNP Kathmandu NAST

10 30 40

27.6˚N V_Model 1

50

Basin model Vs

Vp

Vs

Vp

20 30 40

V_Model 2

50

60

(d)

0.0 0.4 0.8 1.2

20 Depth (km)

27.8˚N

Vp

Velocity (km/s) 0 1 2 3 4 5 6 7 8 9

0 1 2 3 4 5 6 7 8 9 Depth (km)

6

60

E−W

U−D

N−S

88.64 cm

87.65 cm

NAST

179.93 cm

V_Model 1

154.35 cm

123.92 cm

240.54 cm

KATNP

0

20

(e)

40

60

80

0

20

Time (s)

40

60

80

0

20

Time (s)

40

60

80

Time (s) 88.64 cm

87.65 cm

NAST

179.93 cm

V_Model 2

154.35 cm

123.92 cm

240.54 cm

KATNP

0

20

40

60

80

Time (s)

0

20

40 Time (s)

60

80

0

20

40

60

80

Time (s)

Fig. 5. Velocity model tests. (a) The map view of the Kathmandu basin, the blue triangles show stations located in this basin. (b) The 1D velocity model in the source region (Avouac et al., 2015). (c) The combined 1D velocity model with considering the basin model. (d) and (e) Three-component waveform comparison between the data (black) and the synthetics (red), the synthetics are computed using the V_model 1 and V_model 2, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

waveforms at stations KATNP and NAST with those simulated using V_Model 2. As shown in Fig. 5d, the synthetic seismograms calculated from V_Model 1 cannot well fit the observed waveforms on the horizontal components. In the observed seismograms, we can clearly found some reverberation phases after the first arrival phase, which may be caused by the site effects of the thick sediment in the source region. Since V_Model 1 doesn't involve the sedimentary layers, the observed reverberation phases are not well modeled. On the other hand, the synthetic seismograms calculated with V_Model 2 can fit the observations much better, viewing from Fig. 5e, we can find that most of the reverberation phases are fitted very well. Based on this test we make sure that the thick sediment in the source region plays important role in the near field strong motion, it is necessary to take into account the sedimentary layers effects when using the near field waveform observations in the rupture process inversion. For this reason, we use V_Model 2 as the velocity model to calculate the Green's functions for stations KATNP and NAST to overcome the site effects of the sedimentary layers in the joint inversion. In addition, accurate estimation of earthquake rupture velocity is of importance for slip history and stress drop, because it is an important dynamic parameter of faulting process. To resolve this parameter, we

conducted finite fault inversions with different constant rupture velocities. Rupture velocity versus normalized waveform misfit is shown in Fig. 6. The waveform misfit has a minimum value of 2.8 km/s, which is consistent with the result of Avouac et al. (2015). So the average rupture velocity was allowed to vary from 1.5 to 3.5 km/s during the finite fault inversion, which was also used in the final joint inversion. The inversion result of seismic waveform data is shown in Fig. 3b (hereinafter called Model II). In the inversion, the observed waveforms are fitted very well. Model II reveals that the Gorkha earthquake exhibits a unilateral eastward propagating rupture with pure thrust mechanism. The slip distribution within the asperity is heterogeneous with a larger slip concentrating in its upper left of the hypocenter, which is at odds with the relatively robust geodetic solution. The peak slip value is about 5.6 m. The total seismic moment is 7.82 × 1020 N m. The blue line in Fig. 3d shows the inverted moment rate function. Most of the seismic moment was released within the first 80 s, which is slightly longer than the result of Yagi and Okuwaki (2015). Moreover, the teleseismic data based solutions of Wang et al. (2015), Fan and Shearer (2015) and Yagi and Okuwaki (2015) differ significantly from our waveform solution not only in the slip location but also in the slip pattern. The differences are mainly attributed to a limited resolution

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

1.0

Normalized Fit Error

0.9

0.8

Teleseismic

0.7

High−rate GPS & SM Combined 0.6 0

1

2

3

4

5

Rupture Velocity (km/s) Fig. 6. Normalized inversion residual versus rupture velocity for different dataset.

of teleseismic data and the absence of near field datasets in their rupture model. The waveform source model is consistent with the first order with the geodetic data in terms of the length of the rupture and slip magnitude (Fig. 3). However, because of the numerous trade-offs among the model parameters, when we use the Model II to simulate previous geodetic datasets, the seismic-only models do not match the near field geodetic data well (Fig S3 and S4). This is simply because, although seismic data well constrain the temporal process, the shape of the moment release (relative to the origin time), and the pattern of radiated energy, they do not constrain the spatial distribution of moment release, especially when the distribution of near field data is sparse. Although including regional waveforms from strong motion or high-rate GPS stations can provide strong constraints on slip evolution models (Wei et al., 2012), it also requires detailed information about 3-D velocity structure and fault geometry. 3.3. Preferred rupture model of the mainshock The consistency between Model I and Model II suggests that it is appropriate to constrain the slip model using seismic and geodetic data together. We also have conducted checkerboard tests to analyze the resolution of our joint inversion for different scale slip patterns. The test results are plotted in Fig. 4. As expected, the entire slip patch with different scales can be well retrieved. Compared the geodetic checkerboard tests, even the resolution of deep part for the scales 20 km × 20 km is obviously increased (Fig. 4). This proves that the joint inversion can resolve the rupture slip with spatial scale at least 20 km × 20 km. Our optimized joint inversion model, hereinafter as Model III, shows that the dominant mechanism is pure thrust with a maximum slip of 5.8 m (Fig. 3c). The largest slip patch is located to the left of the hypocenter, extending ~60 km along the dip and ~150 km along the strike. As shown in Fig. 3c, the rupture may not propagate to the ground surface. The total seismic moment is 7.87 × 1020 N m, equals to a magnitude of Mw 7.9. Most of the seismic moment is released within the first 80 s, and the inverted moment rate functions of Models III and II are generally agree with each other (Fig. 3d). Fig. 7 shows the snapshots of Model III with time interval of 10 s, which reveals a complex rupture process. At the beginning period of rupture, especially in the first 10 s, the initiated failure on the west side of the hypocenter is relative slow, with speed of 1–1.5 km/s, while the failure on the east side is faster (~ 2.0–2.5 km/s). During the first 20 s, the rupture propagated in a speed apparently larger than 2.0 km/s, and exceeded the rupture front at

7

20 s. The apparent rupture propagation is significantly speeded up to 3 km/s from 20 s to 40 s, accompanying with larger slip rate, which are consistent with the result of Grandin et al. (2015). The moment rate rapidly increases and reaches 2.0 × 1019 Nm/s at 35 s from the initial rupture after the initial stage, and then rapidly decreased with time (Fig. 3d); and the main rupture exhibit unilateral towards the east along the strike. The reason for the unilateral rupture may be caused by the accumulated strain in the surrounding area of the epicenter. In 1505, a Mw 8.2 earthquake occurred in the west of the epicenter (Fig. 1; Ader et al., 2012), which might partially release the tectonic strain, while in the east side of the epicenter, no big earthquakes occurred, and thus the accumulated stress should be higher. Of course, other factors, such as fault geometry and barriers in the fault plane, may also affect the rupture pattern. Our preferred model can well explain the static displacements (Fig. 8), and fit both strong motion and high-rate GPS observations very well. Larger mismatch between the synthetic seismogram and the corrected displacement in the N-S component of SYBC and the E-W component of RMTE are seen in Fig. 9, which is likely due to the local site effects that we didn't consider in our modeling. A comparison of synthetic and teleseismic datasets are given in Fig. 10, all datasets are in good consistence. As shown in Fig. 11a and d, the InSAR observations provide good spatial coverage in and around the source region, and the line-of-sight deformations are well fitted for both the descending (Fig. 11b) and ascending path (Fig. 11e) with small residuals (Fig. 11c and f). The largest misfits are located at the far field of the main rupture, which are probably caused by contaminations of the InSAR observations. Since the line-of-sight deformation data observed by the ALOS-2 satellite is acquired 8 days after the mainshock, the afterslip and the early aftershocks could contaminate the InSAR data and cause some disturbances for the inversion. To further verify the reliability of our rupture model, we compared our results with previous studies. The geodesy-based inversions have similar along-strike slip patterns as our results, but there are differences in the along-dip slip distributions. In the inversions of Wang et al. (2015); Feng et al. (2015), and Diao et al. (2015), there is no significant slip near the surface, which is consistent with our Model III. However, the rupture model of Lindsey et al. (2015) has ~2 m of slip near the surface in the north, as does the solution of Galetzka et al. (2015). Due to the constraints from InSAR, static GPS and high-rate GPS observations, the slip distribution in our preferred model is more dispersive than the teleseismic data based models by Wang et al. (2015) and Zhang et al. (2015). Grandin et al. (2015) combined teleseismic waves, strong motion data, high-rate GPS, static GPS, and synthetic aperture radar (SAR) data to invert for the slip model of the Gorkha earthquake, and obtained good fits to all datasets, even for the stations located in the Kathmandu basin, which may result in part from the lower data frequency in their joint inversion. 3.4. Rupture process of the Mw 7.3 strong aftershock It is noteworthy that aftershock sequence can help improve the understanding of the seismogenic environment of the mainshock. Moreover, the relationship between the mainshock and the strong aftershock can also provide important information on the seismic risk on May 12, 2015 the strongest aftershock (Mw 7.3) occurred in the east side of the mainshock. We also investigated the rupture model for this aftershock, analyzed the relationship between the mainshock and the strongest aftershock, and evaluate the potential seismic risk in the source region by taking into account all the information about the mainshock, the strongest aftershock and the aftershock sequence. For the Mw 7.3 strong aftershock, we also collected the line-of-sight deformation data from the ALOS-2 satellite (Lindsey et al., 2015), and 1120 re-sampled points were used during the joint inversion. Moreover, we selected 37 P and 27 SH teleseismic waves with even azimuthal coverage and high signal-to-noise ratio to constrain the slip distribution of this strongest aftershock. All datasets are shown in Fig. S5. For P waves,

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

8

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

84˚E

85˚E

86˚E

87˚E

84˚E

85˚E

86˚E

87˚E

84˚E

85˚E

86˚E

87˚E

29˚N 10−20s

20−30s

30−40s

40−50s

50−60s

60−70s

70−80s

80−90s

0−10s

28˚N km 0

27˚N 29˚N

40 80

28˚N

27˚N 29˚N

28˚N

27˚N Slip (cm) 0

60 120 180 240 300 360 420 480 540 600

Fig. 7. Snapshots in a time interval of 10 s. The color shows the fault slip. The white and red dashed contours denote the pseudo-rupture front for a rupture velocity of 2.0 km/s and 3.0 km/s respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

we used ground velocities instead of displacements to improve the spatiotemporal resolution of the inversion (Wald et al., 1996). For SH waves, we converted the observations into displacements. All the selected body waves are band-pass filtered from 0.0033 Hz to 1 Hz. The aftershock fault plane is parameterized with 15 and 13 subfaults along strike (312°) and dip (11°), respectively, and the total area of fault model is 60 km × 52 km. The rupture is supposed to initiate at the NEIC hypocenter (27.819°N, 86.080°E, 15 km). The result shows that the slip is concentrated in a compacted area (Fig. 12a), and the amplitude of peak slip in

80˚E 30˚N

82˚E

84˚E

(a)

86˚E 80 cm Obs Syn

the model is ~5.0 m with a static stress drop of about 4.6 MPa. The peak slip of the aftershock is close to the mainshock (Fig. 3c and Fig. 12a), which could be the main reason why the shaking of the aftershock is so strong and caused relatively heavy damages (Dixit et al., 2015; Martin et al., 2015). The total seismic moment is 8.5 × 1019 N m (Mw 7.3), and most of the seismic moment is released in the first 20 s (Fig. 12b). The average rupture velocity in the first 5 s is about 2.5 km/s, and then slows down rapidly with time. In general, all of the InSAR data and seismic waveforms are well fitted (Fig. S6 and S7).

88˚E 80˚E

82˚E

84˚E

(b)

Mw7.9

86˚E Obs Syn

88˚E 50 cm

Mw7.9

28˚N

km 26˚N

0

75 150

km 0

75 150

Fig. 8. Comparison between the observed and the synthetic displacements. (a) The synthetic (red) and data (black) for the horizontal GPS components. (b) Vertical GPS synthetics (red) and the data (black), the rectangle is the assumed fault plane. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

E−W

9

U−D

N−S

High-rate GPS 11.72 cm

14.37 cm

18.64 cm

37.63 cm

80.04 cm

54.94 cm

37.15 cm

25.69 cm

31.53 cm

87.65 cm

179.93 cm

88.64 cm

94.86 cm

199.62 cm

157.00 cm

123.92 cm

240.54 cm

154.35 cm

SYBC

SNDL

RMTE

NAST

KKN4

Strong-motion data

KATNP

0

20

40

60

80 0

20

Time (s)

40

Time (s)

60

80 0

20

40

60

80

Time (s)

Fig. 9. Comparison of the high-rate GPS and strong-motion records (black line) and synthetic seismograms (red line) derived from our model. Both data and synthetics are aligned by the first P arrivals. The number at the first of each seismogram indicates the station name; the number at the right top is the maximum displacement of the records in cm. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In order to verify whether the aftershocks were triggered by the mainshock or not, we used the program PSGRN/PSCMP (Wang, 1999; Wang et al., 2006) to calculate the coseismic stress changes of the mainshock with an effective friction coefficient of 0.4. Usually, large earthquake occurrences represent the back-ground stress of these areas, so, we use the focal mechanism of the Gorkha earthquake (293°, 7° and 108° for strike, dip and rake) as the optimal mechanism of the given background stress during the calculation. The result shows that Coulomb failure stress (CFS) increment in the area of the strongest aftershock is 0.19 MPa (Fig. 13b profile AA′), which is much larger than the threshold of earthquake triggering, ~ 0.01 MPa (Heidbach and Ben-Avraham, 2007; King et al., 1994). The distribution of other aftershocks which occurred after the mainshock and before the Mw 7.4 aftershock is consistent with the coseismic CFS changes caused by the mainshock. In Fig. 13a, we can find that the aftershocks occurred in the increased CFS areas, especially in the un-rupture gap of the mainshock, and the surrounding regions of the major slip areas. So we speculate that the Mw 7.3 aftershock and most of the early aftershocks were triggered by the mainshock. Along the three profiles shown in Fig. 13a, the distributions of the aftershocks are a little bit different from the pattern along the fault. Along A-A′ profile in the east end of the rupture zone, we can find that most of the aftershocks with magnitudes larger than 6 occurred in the CFS increased regions, except some small part of aftershocks occurred in the shadow region which is close to the fault plane. However, in B-B′ and C-C′ profiles, the patterns are different. Besides the aftershocks occurring in the CFS increased regions, a number of aftershocks with magnitudes smaller than 6 occurred in the shadow zone, especially in the upper side of the fault (Fig. 13b). This phenomenon is consistent with the results given by Chan et al. (2016) and Feng et al. (2015). The seismogenic structure in the source region is very complex, consisting of the Main Boundary Thrust (MBT), Main Central Thrust (MCT) and Main Himalayan Thrust (MHT) (Yin, 2006); and the MBT and MCT is quite close to each near the source region, which is also identified by Chan et al. (2016). So, the earthquakes in

the shadow regions might occur on the other thrusts instead of on the fault of the mainshock. 3.5. Comparison of coseismic slip and aftershocks distribution The 2015 April 25 Nepal earthquake was followed by a strong aftershock activity. Hundreds of aftershocks with magnitude greater than ML 4.0 were recorded by the Nepal seismological network. Here, we discussed the aftershocks in the first 45 days following the mainshock (Adhikari, et al., 2015), which is assumed to be long enough to identify the final stable patterns of an aftershock sequence (Das and Henry, 2003). During this period of time, the catalogs identify 553 events with magnitude ML ranging from 4.0 to 7.3. Comparison of our preferred rupture model with aftershock characteristics from the Nepal seismological network is shown in Fig. 14. The rupture distributions of the mainshock and the strong aftershock exhibit a strong complementary. Aftershock clusters are found at the edge of unbroken areas, and in the regions of rapid transition from high to low slip within the main fault area. Moreover, the lack of large thrust aftershocks within the coseismic large slip areas is consistent with nearly complete stress release. The low-slip region down-dip from the hypocenter has fewer aftershock activity, but the stress is increased caused by the coseismic rupture (Fig. 13), suggesting that it would be likely to rupture in the near future. Another remarkable phenomenon is that few aftershocks occurred in the shallow depth. Since the rupture does not reach the surface, the potential of shallow events in the future is high, or the accumulated strain in the shallow depth could be released in aseismic creep. Another cluster of aftershocks are located in the eastern part of the model with no coseismic rupture, which is just in the area where the largest aftershock occurred, and clearly complementary with the areas of coseismic slip. In 1833, an event occurred in the mainshock area (Ader et al., 2012). Ever since then, before the 2015 April 25 Nepal Gorkha earthquake, there is a seismic gap in this area. So the

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

10

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

(a) 46.2

P 239 LSZ 69

52.4

P 347 KBS 59

171.7

P 91 GUMO 57

65.4

P 222 ABPO 59

83.4

P 343 SFJD 79

108.6

P 85 KWAJ 79

93.0

P 202 CRZF 80

69.8

P 338 KEV 52

164.6

P 76 WAKE 74

95.7

P 200 DGAR 37

162.2

P 324 KONO 58

132.2

P 65 MAJO 45

175.6

P 189 PAF 78

79.2

P 321 ESK 66

111.3

P 61 MIDW 83

93.5

P 186 AIS 65

108.3

P 312 ECH 61

94.6

P 50 YSS 48

226.5

P 162 COCO 41

149.0

P 305 PAB 72

81.9

P 42 PET 58

172.9

P 150 NWAO 68

81.1

P 299 ANTO 44

105.9

P 39 ADK 73

123.8

P 141 MBWA 59

84.5

P 285 TAM 70

51.6

P 35 MA2 54

209.8

P 132 CAN 87

47.2

P 282 KOWA 82

47.2

P 28 YAK 44

274.5

P 129 WRAB 67

74.1

P 275 DBIC 86

25.8

P 27 KDAK 81

113.7

P 121 CTAO 76

59.9

P 271 RAYN 35

31.9

P 20 COLA 78

137.8

P 112 PMG 70

68.1

P 250 MBAR 59

43.0

P 16 TIXI 49

244.7

P 110 DAV 43

115.1

P 242 TSUM 80

39.7

P 105 HNR 81

−20 0

20 40 60 80 100 120

−20 0

20 40 60 80 100 120

−20 0

Time (s)

Time (s)

20 40 60 80 100 120

Time (s)

(b) SH 76 WAKE 74

812.2

SH 65 MAJO 45

965.2

SH 61 MIDW 83

510.3

SH 50 YSS 48

747.2

SH 42 PET 58

296.5

SH 39 ADK 73

434.7

SH 35 MA2 54

SH 242 TSUM 80

277.0

SH 222 ABPO 59

277.1

SH 347 KBS 59

182.1

SH 150 NWAO 68

487.3

SH 338 KEV 52

307.6

SH 141 MBWA 59

786.6

SH 324 KONO 58

362.5

SH 132 CAN 87

583.1

SH 312 ECH 61

299.0

SH 129 WRAB 67

994.8

SH 299 ANTO 44

457.4

553.8

SH 121 CTAO 76

858.5

SH 285 TAM 70

319.8

SH 28 YAK 44

476.6

SH 112 PMG 70

1036.6

SH 275 DBIC 86

513.4

SH 27 KDAK 81

468.7

SH 105 HNR 81

813.1

SH 271 RAYN 35

428.8

SH 20 COLA 78

264.7

SH 250 MBAR 59

309.1

−20 0

20 40 60 80 100 120

Time (s)

876.0

SH 91 GUMO 57

−20 0

20 40 60 80 100 120

Time (s)

−20 0

20 40 60 80 100 120

Time (s)

Fig. 10. Comparison of teleseismic P (a) and SH (b) displacement records in black and synthetic seismograms in red predicted by the slip model. Both data and synthetic seismograms are aligned on the P and SH arrivals. The number at the end of each trace is the peak velocity of the data in micrometers per second. The azimuth and distance in degrees are shown at the beginning of each record with the azimuth on top. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

84˚E

86˚E

88˚E

30˚N

(a)

84˚E

Data

Descending Path 48

86˚E

88˚E Model

(b)

11

84˚E

86˚E

88˚E 30˚N Misfit

(c)

28˚N

28˚N

km

km 0

26˚N

100

km 100

0

−50

0

50

100

30˚N

(d)

−100

Data

Ascending Path 157

26˚N

cm

cm −100

100

0

−50

0

50

100

−10

Model

(e)

cm −5

0

5

10 30˚N

Misfit

(f)

28˚N

28˚N

km

km 0

26˚N

0

100

km 100

0

−50

0

50

100

26˚N

cm

cm −100

100

−100

−50

0

50

100

cm −10

−5

0

5

10

Fig. 11. Comparison between observed sight displacements of the mainshock and those predicted from our best-fitting model. Each data point is represented by a color coded circle. (a) and (d) Line-of-sight deformation data from the ALOS-2 satellite along descending path and ascending path respectively. (b) and (e) Predicted line-of-sight deformation from the preferred model. (c) and (f) Misfit between the InSAR data and modeled result, the black rectangle indicates the surface projection of the fault plane.

distribution of the aftershocks and the mainshock is consistent with the seismogenic environment in this region. Considering the slip distribution of the mainshock and the strongest aftershock, as well as the distribution of the other aftershocks, most of the accumulated strain in the seismic gap should be released, so that the seismic risk in this region should decrease because of the mainshock and the aftershocks. On the contrary, few aftershocks occurred around the epicenter, especially in the west part, distributions of fault locking depth and slip rate deficit inverted from GPS measurements show that the strain accumulation on the western Nepal segment is slow and stable (Li et al.,

2016), it may be another reason that only few aftershocks occurred in the west of the hypocenter. However, the effect of stress change triggered by mainshock show positive stress change in the west part (Fig. 13a), suggesting that it would be likely to rupture in the near future when the stress or strain accumulation is larger enough. The lack of aftershocks at one side of the epicenter, and their concentration above or in the one side of the regions with high coseismic slip is similar to what was also observed for several continental events such as, The 1999 ChiChi earthquake and the 2008 Wenchuan earthquake (Chang et al., 2000; Huang et al., 2008). Nevertheless, the lack of aftershocks in the region

Strike= 312˚ Dip= 11˚

km

(a) -30

-20

-10

(b) 10

20

30

20

30

20

10

20

30 40

10

10

50 20 10

60

Slip (cm) 0

100

200

300

400

Moment rate (1019N•m/sec)

20

10

Distance along dip (km)

10

0

Moment rate function 1.0

0.5

0.0 0

10

20

30

40

Time (s)

500

Fig. 12. (a) The Mw 7.3 aftershock slip distribution calculated from joint inversion. The direction of the fault plane is indicated by the long black arrow on the top of the figure, and the star indicates the hypocenter, the color bar scales the slip amplitude, the white arrows represents the slip directions and contours display the rupture initiation time in second. (b) The moment-rate function of the Mw 7.3 aftershock. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

12

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

(b) Depth (km) Depth (km) Depth (km)

(a) C’

29˚N

B’ A’

28˚N

C 27˚N

B

km 0

40

A

80

84˚E

0 10

MPa

A

A’

B

B’

C

C’

20 30 0 10 20 30 0 10 20 30 50

85˚E

86˚E

10.00 5.00 1.50 0.50 0.10 0.03 0.01 −0.01 −0.03 −0.10 −0.50 −1.50 −5.00 −10.00

100

150

Distance (km)

87˚E

Fig. 13. (a) The coseismic coulomb stress change on the fault plane. The white dots indicate the aftershocks before the May 12, Mw 7.3 aftershock (Adhikari et al., 2015). The red stars show the positions of the mainshock and Mw 7.3 aftershock, the yellow stars indicate the aftershocks with magnitude range of 5.5–6.0, and the green stars indicate the aftershocks with magnitude larger than 6.0. The lines AA′, BB′ and CC′ are the locations of the vertical profiles. (b) The three vertical profiles show the Coulomb failure stress (CFS) caused by the mainshock. The stars indicate the aftershock sequence shown in figure (a). The black line indicate the fault trace used in this study. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

with increased Coulomb failure stress still may indicate the high seismic hazard there, which has been witnessed by the 2013 Lushan earthquake which occurred in the southwest of the 2008 Wenchuan earthquake. 3.6. Seismic hazard assessment The active deformation of the Himalayan Thrust belt has been witnessed by the occurrence of repeated earthquakes. At least two large earthquakes (M N 8.0) might take place in the areas west of Kathmandu in Nepal (Molnar, 1987; Bilham et al., 1995; Bettinelli et al., 2006), and large earthquakes (M N 8.6 with rupture length of ~ 400 km) are also speculated in eastern Nepal (Feldl and Bilham, 2006; Sapkota et al., 2013). Based on the distribution of the slip patterns between the mainshock and the strongest aftershock as well as the

83˚E

84˚E

85˚E

86˚E

distribution of the aftershock sequence, we proposed that most of the strains, accumulated from the epicenter of the mainshock to the east end of the strong aftershock are released by the mainshock and the aftershocks, and thus the risk of strong earthquakes in this area should be low. However, considering the features of the Coulomb failure stress caused by the mainshock and the strongest aftershock, although a large number of aftershocks have occurred after the strongest aftershock (Fig. 15), there are still some areas with high potential of active seismicity in the surroundings of the mainshock and the strongest aftershock, especially in the Regions A, B and C, which would be likely to rupture in the near future (Fig. 15). A convergence motion of about 36 mm/yr (Altamimi, 2009) along the Himalayan arc represents more than 7 m of displacement accumulated since the last big interplate subduction event in this area over

87˚E 30˚N

Depth=15 km µ=0.4

MPa 10.00 1.00 0.50 0.30 0.10 0.03 0.01 −0.01 −0.03 −0.10 −0.30 −0.50 −1.00 −10.00

29˚N Mw7.9

A 28˚N

B

28˚N Mw7.3

C

km 27˚N 0

40

km

80 26˚N

Slip (cm) 0

100 200 300 400 500 600

Fig. 14. The surface projection of the inverted slip distribution. White dots indicate the aftershocks from Adhikari et al. (2015). The red contour shows the slip (N1.0 m) distribution of the Mw 7.3 aftershock occurred on May 12, 2015. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

0 40 80 84˚E

86˚E

88˚E

Fig. 15. The Coulomb failure stress (CFS) changes at 15 km depth caused by the mainshock and the strongest aftershock. The white dots indicate the aftershocks from Adhikari et al. (2015), the yellow stars show the epicenters of the mainshock and the Mw 7.3 aftershock. Regions A, B and C the three main areas with increased Coulomb stress without obvious aftershocks. The black rectangles indicate the surface projection of the mainshock and the strongest aftershock fault planes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

180 years ago (1833 Gorkha earthquake), while our results show that the maximum slip of 2015 Gorkha earthquake is only about 5.8 m. There are two possible reasons for this discrepancy: one reason is that the rest of the strains are probably released by the aftershocks or afterslip, and the other one is probably that some of the accumulated strains have been released during the earthquake of 1934 that occurred in the east of Nepal, which yield ~4 m slip (Bilham et al., 2001). In addition, compared with seismic gap in central Nepal (Bilham et al., 2001), this Gorkha earthquake ruptured only part of the seismic gap. There still remains an earthquake vacant zone with a total length of ~150 km west of the mainshock (Fig. 1). Moreover, from the distribution of historical earthquakes (Fig. 1), we can see that few large earthquakes occurred in this segment in recent 500 years. In addition, an interesting phenomenon is that the aftershocks only migrate towards the east direction after the main shock of the 2015 Gorkha earthquake, but few aftershocks occurred in the west. Although the inter-seismic locking studies show that the stress and strain accumulation is relative slow on this segment, and seismic activity of this segment is still in a stable stage (Chan et al., 2016), based on our CFS analysis, we clearly found that the west region of the mainshock shows positive stress change caused by the mainshock (Fig. 15). Therefore, the 2015 Gorkha earthquake may have improved the seismic hazard in this segment. If the slip distributions were well known for the 1505 Mw 8.2 and 1934 Mw 8.1 events, then a slip versus time analysis could be done, which could provide insights into the interseismic cycle beneath the central Nepal; however, there is no knowledge of the slip distributions for the 1505 and 1934 events, except for their estimated magnitudes and rough locations. Any estimation of slip amount for these two ancient events requires assumptions about stress drop or rupture scale, as well as location of the slip. Anyway, geodetic inversions for recent slip deficit suggest relatively uniform coupling in the central Nepal (Ader et al., 2012; Li et al., 2016). Given the paucity of large earthquakes over the previous decades and the current high slip rate deficit, thus, seismic hazard on the eastern Nepal areas remains, failure in a single event is certainly plausible, but the possibility of several smaller events or even a larger event releasing strain cannot be ruled out. However, because of the complex tectonic environment of the Main Himalayan Thrust (Upreti, 1999; Hodges, 2000), there are a large number of active faults near the source region, and the mechanisms of these faults could be quite different (Shanker et al., 2011). Thus, as receiver faults, their Coulomb failure stress changes should be quite different because of the different mechanisms of the faults. So that the Coulomb failure stress changes in the surrounding receiver faults should be more complex than the Coulomb stress pattern shown in Fig. 15. For this reason, the detailed seismic hazard could be a little different to the present analysis, But in anyway, our work can provide an first order reference for future seismic assessment. 4. Conclusion The rupture process of the April 25, 2015 Gorkha earthquake and its largest aftershock that occurred on May 12, 2015 have been investigated by joint inversion of seismic data and geodetic observations, including teleseismic seismograms, near field strong motion data, GPS data and InSAR data. The mainshock exhibits a unilateral propagation pattern, which ruptured a 150 km × 60 km area and released a seismic moment of 7.87 × 1020 N m, equivalent to a magnitude of Mw 7.9, and most of the moment rates were released within the first 80 s, especially at the first 40 s. The rupture speed is relatively slow at the beginning of the rupture (about 11.5 km/s), and accelerates to 2.0–3.0 km/s to the east end. Its largest aftershock ruptured in concentrate region, with the peak slip of ~5.0 m. The magnitude is Mw 7.3, and most of the seismic moment is released in the first 20 s. The static stress drop of the aftershock is ~4.6 MPa, which is a little larger than that of the mainshock (~4.2 MPa). The analysis of Coulomb failure stress shows that the largest aftershock and the most of the early aftershocks should be triggered by

13

the mainshock. From the rupture models of the mainshock and the largest aftershock, there is a good complementary between the rupture slips of the mainshock and the strongest aftershock. Moreover, the majority of aftershocks are located in low-slip regions of the coseismic rupture of the mainshock. By analyzing the rupture slip distributions of the mainshock and the strongest aftershock, the distribution of aftershocks, and the Coulomb stress field of the mainshock and the strongest aftershock, as well as the historical earthquake activities in the MHT, we proposed that the seismic risk in the east end of the fault could decrease, while in the west side of the epicenter, there is a high potential for seismicity, and special attention should be paid there. Acknowledgements The teleseismic data used in this work are downloaded from IRIS, the strong motion data are provided by USGS, the GPS data processed by the JPL/Caltech ARIA project. The InSAR data were compiled by Lindsey et al. (2015) from the ALOS-2 satellite. We show our respect to Prof. Yingjie Yang in Macquarie University for careful revision of this manuscript and good suggestion. We also show our respect to Prof. Rongjiang Wang in GFZ for providing us tools for correcting baseline shifts of the accelerograms. We also thank Dr. Liming Jiang and his student for providing us some useful InSAR data for this work. This work is supported by NSFC funding (41541034, 41422401, and 41321063), grant from Chinese Earthquake Administration (201308013), China National Special Fund for Earthquake Scientific Research in Public Interest (201308011) and Extinguished Young Scientist funding from Hubei Province (2012FFA026). Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.tecto.2016.05.034. References Ader, T., Avouac, J.P., Liu-Zeng, J., et al., 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: implications for seismic hazard. J. Geophys. Res. Solid Earth 117 (B4) (1978–2012). Adhikari, L.B., Gautam, U.P., Koirala, B.P., et al., 2015. The aftershock sequence of the 2015 April 25 Gorkha–Nepal earthquake. Geophys. J. Int. 203 (3), 2119–2124. Altamimi, Z., 2009. The International Terrestrial Reference Frame (ITRF2005). Geodetic Reference Frames, Int. Assoc. of Geod. Symp., Vol. 134, Part 2. Springer, Dordrecht, Netherlands, pp. 81–82. Avouac, J.P., Meng, L., Wei, S., Wang, T., Ampuero, J.P., 2015. Lower edge of locked Main Himalayan Thrust unzipped by the 2015 Gorkha earthquake. Nat. Geosci. 8, 708–711. Bettinelli, P., Avouac, J.P., Flouzat, M., et al., 2006. Plate motion of India and interseismic strain in the Nepal Himalaya from GPS and DORIS measurements. J. Geod. 80 (8– 11), 567–589. Bilham, R., 2004. Earthquakes in India and the Himalaya: tectonics, geodesy and history. Ann. Geophys. 47, 839–858. Bilham, R., Bodin, P., Jackson, M., 1995. Entertaining a great earthquake in western Nepal: historic activity and geodetic test for the development of strain. J. Nepal Geol. Soc. 11, 73–88. Bilham, R., Gaur, V.K., Molnar, P., 2001. Himalayan seismic hazard. Science 293 (5534), 1442–1444. Chan, C.H., Wang, Y., Almeida, R., et al., 2016. Enhanced stress and changes to regional seismicity due to the 2015 Mw 7.8 Gorkha, Nepal, earthquake on the neighboring segments of the Main Himalayan Thrust. J. Asian Earth Sci. 2016. http://dx.doi.org/ 10.1016/j.jseaes.2016.03.004. Chang, C.H., Wu, Y.M., Shin, T.C., Wang, C.Y., 2000. Relocation of the 1999 Chi-Chi earthquake in Taiwan. Terr. Atmos. Ocean. Sci. 11, 581–590. Das, S., Henry, C., 2003. Spatial relation between main earthquake slip and its aftershock distribution. Rev. Geophys. 41 (3). Diao, F., Walter, T.R., Motagh, M., Prats-Iraola, P., Wang, R., Samsonov, S.V., 2015. The 2015 Gorkha earthquake investigated from radar satellites: slip and stress modeling along the MHT. Front Earth Sci. 3, 65. http://dx.doi.org/10.3389/feart.2015.00065. Dixit, A.M., Ringler, A., Sumy, D., Cochran, E., Hough, S.E., Martin, S.S., Gibbons, S., Luetgert, J., Galetzka, J., Shrestha, S.N., et al., 2015. Strong motion observations of the M 7.8 Gorkha, Nepal, earthquake sequence and development of the N-SHAKE strong-motion network. Seismol. Res. Lett. 86. http://dx.doi.org/10.1785/0220150146 (no. 6). Engdahl, E.R., van der Hilst, R., Buland, R., 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bull. Seismol. Soc. Am. 88 (3), 722–743.

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034

14

C. Liu et al. / Tectonophysics xxx (2016) xxx–xxx

Fan, W., Shearer, P.M., 2015. Detailed rupture imaging of the 25 April 2015 Gorkha earthquake using teleseismic P waves. Geophys. Res. Lett. 42, 5744–5752. Feldl, N., Bilham, R., 2006. Great Himalayan earthquakes and the Tibetan plateau. Nature 444, 165–170. Feng, G., Li, Z., Shan, X., Zhang, L., Zhang, G., Zhu, J., 2015. Geodetic model of the 2015 April 25 Mw 7.8 Gorkha Gorkha earthquake and Mw 7.3 aftershock estimated from InSAR and GPS data. Geophys. J. Int. 203 (2), 896–900. Fielding, E.J., Sladen, A., Li, Z., Avouac, J.-P., Burgmann, R., Ryder, I., 2013. Kinematic fault slip evolution source models of the 2008 M7.9 Wenchuan earthquake in China from SAR interferometry, GPS and teleseismic analysis and implications for Longmen Shan tectonics. Geophys. J. Int. 194 (2), 1138–1166. Galetzka, J., Melgar, D., Genrich, J.F., et al., 2015. Slip pulse and resonance of the Kathmandu basin during the 2015 Gorkha earthquake, Nepal. Science 349 (6252), 1091–1095. Grandin, R., Vallée, M., Satriano, C., Lacassin, R., Klinger, Y., Simoes, M., Bollinger, L., 2015. Rupture process of the Mw = 7.9 2015 Gorkha earthquake (Nepal): insights into Himalayan megathrust segmentation. Geophys. Res. Lett. 42. http://dx.doi.org/10.1002/ 2015GL066044. Graves, R.W., Wald, D.J., 2001. Resolution analysis of finite fault source inversion using 1D and 3D Green’s functions. I. Strong motion. J. Geophys. Res. 106, 8745–8766. Heidbach, O., Ben-Avraham, 2007. Stress evolution and seismic hazard of the Dead Sea Fault System. Earth Planet. Lett. 257, 299–312. Hodges, K.V., 2000. Tectonics of the Himalaya and southern Tibet from two perspectives. Geol. Soc. Am. Bull. 112 (3), 324–350. Huang, Y., Wu, J.P., Zhang, T., Zhang, D., 2008. Re-determination of the large Wenchuan earthquake(Ms 8.0) and its sequence of its aftershocks (in Chinese). Sci. China Ser. D Earth Sci. 38, 1242–1249. Ji, C., Wald, D.J., Helmberger, D.V., 2002. Source description of the 1999 Hector Mine, California, earthquake, part I: wavelet domain inversion theory and resolution analysis. Bull. Seismol. Soc. Am. 92, 1192–1207. King, G.C.P., Stein, R.S., Lin, J., 1994. Static stress changes and the triggering of earthquakes. Bull. Seismol. Soc. Am. 84, 935–953. Konca, A.O., Leprince, S., Avouac, J.-P., Helmberger, D.V., 2010. Rupture process of the 1999 Mw 7.1 Duzce earthquake from joint analysis of SPOT, GPS, InSAR, strong-motion, and teleseismic data: a supershear rupture with variable rupture velocity. Bull. Seismol. Soc. Am. 100, 267–288. Kumar, V.P., Chauhan, N.S., Padh, H., et al., 2006. Search for antibacterial and antifungal agents from selected Indian medicinal plants. Journal Ethnopharmacol. 107 (2), 182–188. Lavé, J., Yule, D., Sapkota, S., et al., 2005. Evidence for a great medieval earthquake (~1100 AD) in the central Himalayas, Nepal. Science 307 (5713), 1302–1305. Li, Y.C., Song, X., Shan, X., et al., 2016. Locking degree and slip rate deficit distribution on MHT fault before 2015 Nepal Mw 7.9 earthquake[J]. J. Asian Earth Sci. 119, 78–86. Lindsey, E.O., Natsuaki, R., Xu, X., et al., 2015. Line-of-sight displacement from ALOS-2 interferometry: Mw7.8 Gorkha earthquake and Mw7.3 aftershock. Geophys. Res. Lett. 42, 6655–6661. http://dx.doi.org/10.1002/2015GL065385. Liu, C.L., Zheng, Y., Xiong, X., et al., 2015. Rupture process of the 23 October 2011 Mw7. 1 van earthquake in eastern Turkey by joint inversion of teleseismic, GPS and strongmotion data. Pure Appl. Geophys. 172, 1383–1396. Mahesh, P., Rai, S.S., Sivaram, K., et al., 2013. One‐dimensional reference velocity model and precise locations of earthquake hypocenters in the Kumaon–Garhwal Himalaya. Bull. Seismol. Soc. Am. 103, 328–339.

Martin, S.S., Hough, S.E., Bilham, R., Hung, C., 2015. Ground motions from the 2015 Mw 7.8 Gorkha, Nepal, earthquake, constrained by a detailed assessment of macroseismic data. Seismol. Res. Lett. 86. http://dx.doi.org/10.1785/0220150138 (no. 6). Molnar, P., 1987. The distribution of intensity associated with the 1905 Kangra earthquake and bounds on the extent of the rupture zone. J. Geol. Soc. India 29, 211–229. Mukhopadhyay, B., Acharyya, A., Dasgupta, S., 2011. Potential source zones for Himalayan earthquakes: constraints from spatial–temporal clusters. Nat. Hazards 57 (2), 369–383. Paudyal, Y.R., Yatabe, R., Bhandary, N.P., Dahal, R.K., 2013. Basement topography of the Kathmandu Basin using microtremor observation. J. Asian Earth Sci. 62, 627–637. Sapkota, S.N., Bollinger, L., Klinger, Y., Tapponnier, P., Gaudemer, Y., Tiwari, D., 2013. Primary surface ruptures of the great Himalayan earthquakes in 1934 and 1255. Nat. Geosci. 6, 71–76. Shanker, D., Paudyal, H., N., S.H., 2011. Discourse on seismotectonics of Nepal Himalaya and vicinity: appraisal to earthquake hazard. Geosciences 1 (1), 1–15 (2011). Sladen, A., Tavera, H., Simons, M., et al., 2010. Source model of the 2007 Mw 8.0 Pisco, Peru earthquake: Implications for seismogenic behavior of subduction megathrusts. J. Geophys. Res. 115, B02405. http://dx.doi.org/10.1029/2009JB006429. Upreti, B.N., 1999. An overview of the stratigraphy and tectonics of the Nepal Himalaya. J. Asian Earth Sci. 17, 577–606. Wald, D.J., Heaton, T.H., Hudnut, K.W., 1996. The slip history of the 1994 Northridge, California, earthquake determined from strong-motion, teleseismic, GPS, and leveling data. Bull. Seismol. Soc. Am. 86, S49–S70. Wald, D.J., Graves, R.W., 2001. Resolution analysis of finite fault source inversion using 1D and 3D Green’s functions. II. Combining seismic and geodetic data. J. Geophys. Res. 106, 8767–8788. Wang, R., 1999. A simple orthonormalization method for the stable and efficient computation of Green's functions. Bull. Seismol. Soc. Am. 89, 733–741. Wang, R., Lorenzo, F., Roth, F., 2006. PSGRN/PSCMP — a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic gravitational dislocation theory. Comput. Geosci. 32 (4), 527–541. Wang, R., Schurr, B., Milkereit, C., Shao, Z., Jin, M., 2011. An improved automatic scheme for empirical baseline correction of digital strong-motion records. Bull. Seismol. Soc. Am. 101 (5), 2029–2044. Wang, W.M., Hao, J.L., He, J.K., Yao, Z.X., 2015. Rupture process of the Mw7.9 Gorkha earthquake April 25, 2015. Sci. China Earth Sci. http://dx.doi.org/10.1007/s11430015-5170-y. Wei, S., et al., 2011. Superficial simplicity of the 2010 El Mayor-Cucapah earthquake of Baja California in Mexico. Nat. Geosci. 4, 615–618. Wei, S., Graves, R., Helmberger, D., Avouac, J.P., Jiang, J., 2012. Sources of shaking and flooding during the Tohoku-Oki earthquake: a mixture of rupture styles. Earth Planet. Sci. Lett. 333, 91–100. http://dx.doi.org/10.1016/j.epsl.2012.04.006. Yagi, Y., Okuwaki, R., 2015. Integrated seismic source model of the 2015 Gorkha, Nepal, earthquake. Geophys. Res. Lett. 42, 6229–6235. Yin, A., 2006. Cenozoic tectonic evolution of the Himalayan orogeny as constrained by along-strike variation of structural geometry, exhumation history, and for land sedimentation. Earth Sci. Rev. 76, 1–131. Zhang, Y., Xu, L.S., Chen, Y.T., 2015. Rupture process of the 2015 Nepal Mw7.9 earthquake: fast inversion and preliminary joint inversion. Chin. J. Geophys. 58 (5), 1804–1811 (in Chinese).

Please cite this article as: Liu, C., et al., Rupture processes of the 2015 Mw 7.9 Gorkha earthquake and its Mw 7.3 aftershock and their implications on the seismic risk, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.05.034