Probing time-dependent afterslip and viscoelastic relaxation following the 2015 Mw7.8 Gorkha earthquake based on the 3-D finite-element model

We analyzed daily displacement time series from 34 continuous GPS stations in Nepal and 5 continuous GPS stations in South Tibet, China, and extracted the first 4.8 years of postseismic motion after the 2015 Mw7.8 Gorkha earthquake. With the longer duration GPS observations, we find that postseismic displacements mainly exhibit southward and uplift motion. To study the postseismic afterslip and viscoelastic relaxation, we built a 3-D spherical finite-element model (FEM) with heterogeneous material properties and surface topography across the Himalayan range, accounting for the strong variations in material properties and surface elevation along the central Himalayan arc. On the basis of the FEM, we reveal that the predicted viscoelastic relaxation of cm level moves southward to the north of the Gorkha earthquake rupture, but in an opposite direction to the observed postseismic deformation in the south; the postseismic deformation excluding viscoelastic relaxation is well explained by afterslip downdip of the coseismic rupture. The afterslip is dominant during 4.8 years after the 2015 Mw7.8 Gorkha earthquake; the contribution by the viscoelastic relaxation gradually increases slightly. The lack of slip on a shallow portion and western segment of the MHT during and after the 2015 Gorkha earthquake implies continued seismic hazard in the future.


Introduction
On 25 April 2015, the Mw7.8 Gorkha earthquake struck, just 77 km northwest of Kathmandu, Nepal, resulting in large economic losses and a death toll of up to 8000 (Earthquake Relief Portal 2015). The Himalayan region, accommodating a high convergence rate of 40-50 mm/ year between the overriding Eurasian Plate and the underthrusting Indian Plate (Bilham et al. 1997;Lavé and Avouac 2000), is prone to earthquakes, such as the 1255and 1934Mw8 + Nepal-Bihar earthquakes, 1956 Mw8.6 Assam-Tibet earthquake, 2005 Mw7.6 Kashmir earthquake, etc. (Sapkota et al.2013;Bilham 2004;Wang and Fialko 2014). The Gorkha earthquake occurred near the E-W trending Main Frontal Thrust (MFT) (USGS 2015), the most active branch of the Main Himalayan Thrust (MHT), which absorbs about half of MHT convergence (Lavé and Avouac 2000;Kumar et al. 2001). Studies show that the 2015 Gorkha rupture was concentrated in the deeper part of the seismogenic zone, between depths of 10 and 15 km, leaving the unbroken shallow part at high seismic risk (Lindsey et al. 2015;Wang and Fialko 2015;Mencin et al. 2016;Sreejith et al. 2015;Galetzka et al. 2015).
When an earthquake occurs, sudden coseismic stress changes in the surrounding rocks may lead to viscoelastic relaxation in the lower crust and upper mantle, and Su et al. Earth, Planets and Space (2020) 72:168 aseismic slip updip and/or downdip of the coseismic rupture (Pollitz et al. 2003;Freed et al. 2004Freed et al. , 2007Freed et al. , 2010Rousset et al. 2012;Perfettini et al. 2004;Barbot 2010;Savage and Svarc 1992). The postseismic deformation due to the Gorkha earthquake provides valuable constraints to investigate the lithosphere strength and constitutive properties beneath the Himalayan region and Tibetan Plateau. Therefore, the exploration of postseismic transients and what mechanisms govern the postseismic process is of a broad interest. A few of studies focused on the afterslip modeling to explain the rapid transient postseismic deformation (Gualandi et al. 2016;Sreejith et al. 2015;Jiang et al. 2019). With the longer span of data available, some studies investigated various mechanical processes contributed by afterslip and viscoelastic relaxation with and without poroelastic rebound (Mencin et al. 2016;Zhao et al. 2017;Wang and Fialko 2018;Jiang et al. 2018;Jouanne et al.2019;Tian et al. 2020;Liu et al. 2020). The previous studies mostly employed the homogeneous elastic half-space or layered elastic model to explore afterslip (Sreejith et al. 2016;Zhao et al. 2017; Wang and Fialko 2018;Jiang et al. 2018;Jiang et al. 2019), and the layered model or lateral heterogeneous model, which incorporates a low viscous layer beneath the Tibetan Plateau to investigate viscoelastic relaxation (Zhao et al. 2017;Wang and Fialko 2018;Jouanne et al. 2019;Tian et al. 2020).
Material heterogeneity affects surface deformation due to fault slip (Masterlark et al. 2001(Masterlark et al. 2003Hughes et al. 2010;Kyriakopoulos et al. 2013;Wiseman et al. 2015;Tung and Masterlark 2016;Hu et al. 2016;Klein et al. 2016;Hsu et al. 2014;Hines et al. 2016;Pratama et al. 2017;Suito 2017;Freed et al. 2017). Therefore, the homogenous or layered earth model assumption for the Himalayan range, which has strong variations in material properties (Cattin et al. 2001;Monsalve et al. 2006;Schul.te-Pelkum et al. 2005), inevitably introduces bias into the modeling (Kyriakopoulos et al. 2013;Tung and Masterlark 2016;Suito 2017;). In this paper, we build a 3-D spherical finite-element model (FEM) considering heterogeneous material properties and surface topography across the Himalayan range to investigate the aseismic afterslip and the viscoelastic relaxation following the 2015 Gorkha earthquake, using extracted postseismic displacement from continuous GPS (cGPS) data in Nepal and southern Tibetan Plateau over 4.8 years after the event.

GPS processing
We collected raw RINEX data from 39 cGPS sites within an area between 78-92° E and 23-34° N (Fig. 1). GPS sites in Nepal include 25 sites operated prior to the Gorkha earthquake and 9 sites deployed immediately Fig. 1 Tectonic setting. The black beach ball shows the mainshock focal mechanism from U.S. Geological Survey. The blue circles are continuous GPS sites used in this study. The color represents the topographical relief, and the color of overlapped fault plane represents the coseismic slip distribution from Tung and Masterlark (2016). The black lines mark the depth of faulting. The red lines represent active faults from Taylor and Yin (2009). The Figure was generated by GMT software (Wessel et al. 2019) after the earthquake to capture postseismic deformation. The raw RINEX files of these sites are openly available at UNAVCO (ftp://data-out.unavc o.org). We additionally collected data from 5 cGPS in South Tibet from the Crustal Movement Observation Network of China (CMONOC) (https ://www.cgps.ac.cn/), which began recording in 2011.
We processed the GPS data using the GAMIT/ GLOBK processing software (Herring et al. 2015). Initially, we estimated loosely constrained daily solutions with GAMIT, together with orbits and Earth orientation, tropospheric delay parameters and full covariance matrices. We used data from the 39 cGPS sites from Nepal and South Tibet together with about 80 global distributed core stations of the International GNSS Service (Johnston et al. 2017). Then, we used GLOBK, a smoothing Kalman filter, to transform the loosely constrained solutions to the ITRF2014 reference frame (Altamimi et al. 2016), which is realized by the 80 core stations. The detailed strategy and applied models followed Su et al. (2019).
We estimated daily displacement time series for the 5 GPS sites in Tibet from 2011 until March 5, 2020, and 34 sites in Nepal from 2011 until January 26, 2019. We compared our time series with solutions from the Nevada Geodetic Laboratory (NGL) (Blewitt et al. 2018) and found the two data sets to be in good agreement (Additional file 1: Fig. S1). Since the GPS coordinate time series for sites in Nepal routinely provided by NGL have longer duration, some of those begin in the 1990s, we used in the time series modeling the NGL solutions for these sites. Considering the possible slight deviation due to the different reference frame realizations of the two solutions, we estimated seven transformation parameters between the NGL and our GAMT/GLOBK solutions to transfer the NGL solutions from IGS 08 (Rebischung et al. 2012) to ITRF2014.

Time series analysis
The obtained daily time series reflects the surface motion due to the velocities, coseismic displacements, postseismic parameters, seasonal terms (annual and semiannual), and non-tectonic offsets primarily due to any changes in antennas and receivers. The component (∆N, ∆E, ∆U) time series at discrete time epochs t i can be modeled independently according to (Nikolaidis 2002): (1) The coefficient a is the value at the initial epocht 0 ,t i denotes the time elapsed (in years) fromt 0 , and the linear rate (slope) b represents the interseismic secular tectonic motion in mm/year. The coefficients c , d , e , and f denote annual and semi-annual variations. The seasonal variations in Himalayan range caused by the hydrological loading are up to 30-50 mm (Fu and Freymueller 2012), and can bias the transient displacements, therefore the periodic terms need to be well dealt with. The magnitudes, g , of n g jumps (offsets, steps, and discontinuities) are due to coseismic deformation and/or non-coseismic changes at epochs T g . H denotes the discrete Heaviside function. With respect to the postseismic deformation term, it is often parameterized by exponential or logarithmic model, or combinations of both. The former is often associated with viscoelastic relaxation in the lower crust and the upper mantle (e.g., Shen et al. 1994;Nikolaidis 2002), while the latter is associated with afterslip within the extended rupture area (e.g., Savage and Svarc 1997;Freed et al. 2010). In this paper, we use the logarithmic parametrization, since related research showed that the afterslip mechanism dominated the Gorkha earthquake (Gualandi et al. 2016;Mencin et al. 2016;Sreejith et al. 2016;Zhao et al. 2017;Wang and Fialko 2018;Jiang et al. 2019). The postseismic term includes the amplitudes k j , postseismic events starting at epochs T k and a decaying constant τ j , which characterizes how fast the postseismic deformation decays with time. The noise term ε i is considered as the sum of errors and unmodeled residuals.
Fitting the time series, we can clearly extract the various motions, in particular, the postseismic deformation. For the sake of reliability and obtaining the postseismic parameters with a high signal-to-noise ratio, two more steps were taken: firstly, we estimated secular velocities of sites which operated for more than 2 years before the earthquake, and combined those with GPS velocities provided by Ader et al. (2012) and Bettinelli et al. (2006), through seven parameter transformations to obtain the integrated velocity field in the ITRF2014 reference frame. We then interpolated the velocities using the VISR software (Shen et al. 2015), and re-estimated the time series parameter (Eq. 1), constraining the prior secular velocities. Secondly, motivated by the work of Savage and Svarc (2009) and Barbot et al. (2009), we employed the principal component analysis (PCA) to extract a common postseismic signal for all the displacement time series and their principal components (Additional file 1: Fig.  S2). The first principal component clearly represented the postseismic process. On account of the non-linearity of the postseismic deformation model parameter τ , we used trial-and-error over a range of decay values to minimize the sum of the squared residuals between the observations and model. We then fit each coordinate time series (Eq. 1) with the optimal τ , and finally extracted the postseismic parameter for each site (Additional file 1: Fig. S3).

Postseismic deformation
The 4.8-year cumulative postseismic displacements ( Fig. 2) associated with the 2015 Gorkha earthquake exhibit similar patterns as the coseismic offsets (Elliott et al. 2016;Gualandi et al. 2016;Mencin et al. 2016). Most GPS sites continue to move southward after the earthquake, showing a dominant thrust feature with minor dextral strike-slip motion along southeast Nepal. As the earthquake ruptured toward the southeast, significant postseismic displacements are primarily distributed across the High Himalaya region. The largest postseismic displacement at site CHLM reaches up to -93.84 mm, -31.87 mm and 36.64 mm for the N, E, U directions, respectively. Displacements in the north of the rupture area are much larger compared to displacements south of the rupture, which has been confirmed by many studies (Zhao et al. 2017;Wang and Fialko 2018;Jiang et al. 2018;Jiang et al. 2019;Jouanne et al. 2019;Tian et al. 2020), implying the possible existence of afterslip downdip of coseismic rupture. In the vertical direction, most sites experienced ongoing uplift after the Gorkha earthquake, while some sites to the south of the rupture area subsided. The extremely large uplift of site NAST is not reliable, since it suffered from rapid subsidence prior to the mainshock (Zhao et al. 2017). Therefore, we exclude the vertical component of site NAST in modeling the fault slip.

Finite-element model
The superiority of finite-element models (FEMs) in simulating the elastic dislocation and the viscoelastic relaxation with heterogeneous material properties has been proven by previous studies (Masterlark et al. 2012;Hughes et al. 2010;Hines et al. 2016;Hsu et al. 2013;Tung and Masterlark 2016;Pratama et al. 2017). In the view of the strong topographical fluctuations and heterogeneity across the Himalayan region, the spherical FEM with topography and heterogeneity was constructed using the CUBIT software (Blacker et al. 2016).
We defined a spherical shell from the mantle to the surface shaped by the irregular topographical geometry in the target region ( Fig. 3) with a resolution of about 1 km obtained from the SRTM30 software (Becker et al. 2009). The 3D FEM extended over 78-92° E and 23-34° N, from the ground surface to a depth of 300 km, incorporating the upper crust, the middle crust, The lower crust and the upper mantle, and the layer interfaces were interpolated based on CRUST 1.0 (Laske 2013). The heterogeneous elastic material properties were characterized by spatially variable P and S wave velocities, and the density based on CRUST 1.0 (Laske 2013). A spatially variable fully anisotropic viscous model was assigned (Additional file 1: Fig.  S4), based on Sun et al. (2013). Furthermore, we assumed zero displacements as side and bottom boundary conditions to take into account their very large extent, along with the stress-free top surface.
The 3D FEM spanned about 1400 km, 1200 km and 300 km in length, width and depth. The blocks were meshed and discretized by 186,079 tetrahedral elements. The elements in the upper crust were of a length of about 11 km, while the element size increased up to a factor of 5 at the deeper depth in the mantle for the sake of computational time saving; the elements near the fault were dense with a length of about 5 km to capture strong variable deformation. Because the FEM configuration used by Tung and Masterlark (2016) basically matched ours, we adopted the same refined coseismic rupture to drive viscoelastic relaxation. On the basis of that, we assigned the same fault geometry with a planar patch centered at (27.95° N, 85.27° E) with a dip and strike of 10° and We performed the calculation with the finite element software PYLITH (Aagaard et al. 2013), which was widely applied to postseismic deformation modeling (Hines et al. 2016;Pratama et al. 2017) and to generate the Green's function (Diao et al. 2014;Hsu et al. 2013;Hines et al. 2016).

Viscoelastic relaxation
Viscoelastic relaxation due to the coseismic stress changes in the lower crust and upper mantle can explain the postseismic deformation at extensive spatial scales (Freed et al. 2004(Freed et al. , 2007Huang et al. 2014;Rousset et al. 2012). The previous studies calculated the viscoelastic relaxation based on the lateral heterogeneous or layered model, and evaluated viscoelastic relaxation in the postseismic process after the Gorkha earthquake (Zhao et al. 2017;Wang and Fialko 2015;Jiang et al. 2018;Jouanne et al. 2019;Tian et al. 2020;Liu et al. 2020). Differing from these models, we designated spatial variable viscosities to account for the heterogeneous rheological properties across the Himalayan region. Furthermore, we simulated the viscoelastic relaxation with the refined rupture model of Tung and Masterlark (2016), which was inverted by a 3D heterogeneous FEM with input from GPS and INSAR data consistent with other published models (Sreejith et al. 2016;Wang and Fialko 2018),

Afterslip
The sudden coseismic deformation triggers the aseismic slip on the extended coseismic rupture area (Marone 1991;Perfettini et al. 2004;Hsu et al. 2006;Barbot 2010). The afterslip plays an important role in the postseismic process from several days to several months (Barbot et al. 2008;Segall 2010;Huang et al. 2014). In order to investigate the distribution of afterslip following the Gorkha earthquake, we also performed a 3D FEM inversion using the cGPS data as input. The fault, 286 km long and 200 km wide, was discretized by a grid of 28 × 20 into 560 subfault patches with a dimension of 10.21 × 10 km, and with Green's functions for each subfault patch due to a unit slip along dip and strike.
With the generated Green's functions, we assumed that the afterslip can be described by a dislocation model of distributed slip. The objective function is: The first term represents the misfit between modeled and observed displacements, where W is the weight matrix inferred from observation uncertainties, G is the Green's functions, m is the estimated slip distribution, and d is the observations. The second term is Laplacian smoothing to avoid abrupt slip variation, where β is smoothness and ∇ 2 is the Laplacian operator.
An algorithm was then developed to invert the smoothed afterslip distribution achieved by a constrained least-squares optimization based on the steepest descent method, motived by Wang et al. (2013). The optimal smoothness was determined by the trade-off between the roughness and data misfit.

Verification of 3D finite-element model and methodology
In order to validate our FEM calculations, we tested a flat layered model with elastic upper and middle crust overlaying a viscoelastic lower crust and mantle, and assigned viscosity values of 1.6 × 10 19 Pa·s for the lower crust and 10 20 Pa·s for the mantle according to Jiang et al. (2018). Comparing the afterslip and viscoelastic relaxation predicted by our FEM with afterslip by the SDM software (Wang 2013) and viscoelastic relaxation by the PSGRN/ PSCMP software (Wang et al. 2006), respectively (Additional file 1: Figs. S5 and S6), we find coherent patterns and equivalent magnitudes of afterslip and viscoelastic relaxation.

Viscoelastic relaxation
We characterized the viscoelastic material as the bi-viscoelastic Burgers model, incorporating a transient relaxing Kelvin element in series with a steady-state Maxwell element. For simplicity, we set a constant ratio of 0.1 between Kelvin and Maxwell viscosities (Hu et al. 2016;Hines et al. 2016;Zhao et al. 2017). The viscoelastic relaxation by the forward model exhibits two quadrants of horizontal displacements, namely, pronounced northward motion in the south of the rupture and southward motion in the north. The directions of the predicted viscoelastic relaxation in the south of the rupture area are opposite to the observations. In the vertical direction, there is obvious subsidence in the near field and uplift in the surrounding areas (Fig. 4). The relatively small horizontal and vertical displacements over 4.8 years after Gorkha earthquake imply that viscoelastic relaxation is not indicative of the postseismic process.
We additionally characterized the rheological property as a Maxwell body, and obtained even less significant viscoelastic relaxation (Additional file 1: Fig. S7). Lacking a Kelvin element, Maxwell rheology usually produces smaller viscoelastic relaxation than a Burgers rheology. In the horizontal component, the Maxwell rheology indicates a similar pattern to that of the Burgers rheology. In the vertical direction, there is subsidence in the rupture area for both models, but insignificant uplift to the south of rupture area.

Afterslip
Next, we modeled the postseismic deformation by afterslip on the fault interface. The knee of the tradeoff curve between the model roughness and data misfit was selected as the optimal smoothing factor (Fig. 5); a smoothing factor of β = 0.06 was chosen in our final solution. With viscoelastic relaxation deducted, the observed horizontal surface displacements are well explained by the downdip afterslip model (Fig. 6). The postseismic displacements due to afterslip are much larger than those associated with viscoelastic relaxation, suggesting that the afterslip plays the dominant role in the postseismic process during the first 4.8 years after the event. The afterslip is mainly distributed downdip of the coseismic rupture, which would be responsible for the much larger postseismic deformation north of the rupture than south of the rupture. The high-slip patch is located at the depth of about 25-30 km with a peak of 0.28 m. Afterslip occurs at the periphery of the coseismic rupture, in the areas characterized by low aftershock activity (Barbot et al. 2009). The moment release by postseismic afterslip is 7.32 × 10 19 N m, equivalent to a moment magnitude of Mw ~ 7.18, approximately 6.7% of the mainshock moment release of 1.09 × 10 21 Nm of the coseismic moment release (Tung and Masterlark 2016).
To examine whether the afterslip is well constrained by the GPS data, we conducted a checkerboard test for the resolution (Additional file 1: Fig. S8). The slip model merged 4 × 4 subpatches into a single pixel to exhibit 7 × 5 checkered asperities (each pixel of 40.85 × 40 km) of 0 and 1.4 m slip. We then compared the checkerboard pattern slip with the slip model inverted by GPS displacements, and found that slip distribution to be well recovered in spite of minimal distortion and discernible artifacts along the down edge of the slipping patches. The far downdip of the fault was not well

Temporal evolution of the deformation and its mechanisms
We divided the postseismic deformation over 4.8 years after the Gorkha earthquake into four equal time periods of 1.2 years (Fig. 7). The overall patterns of observed postseismic deformation, including both afterslip and postseismic relaxation, during the four time periods are  06 m), the third to about 13.9% (0.04 m) and about 10.2% (0.03 m) in the last 1.2 years. Snapshots of the spatial distribution due to afterslip alone indicate that most of it occurs downdip of the coseismic rupture. The viscoelastic relaxation component shows a similar pattern, but smaller slip magnitudes. The afterslip downdip of the coseismic rupture is usually associated with velocity-strengthening behavior, results in enduring aseismic slip following the earthquake (Lienkaemper and McFarland 2017;Tian et al. 2020). Figure 8 shows how the postseismic deformation evolves with time and the individual contributions of afterslip and viscoelastic relaxation. At the site CHLM, the predicted viscoelastic relaxation is opposite to the observation for the N component. For the E component, the contribution rates of viscoelastic relaxation to postseismic deformation at site CHLM during four different time periods are 19.2%, 22.7%, 24.8%, and 26.2%, respectively. At the site XZZF in the Tibet, the contribution rates of viscoelastic relaxation to postseismic deformation during four different time periods are 10.5%, 12.5%, 13.7%, and 14.5% for the N component, and 4.4%, 6.0%, 7.1% and 8.0% for the E component. The viscoelastic relaxation plays insignificant role in the postseismic process over the 4.8 years after Gorkha earthquake, but its contribution to the postseismic deformation gradually increases slightly. Longer spans of GPS displacements after the Gorkha earthquake show that postseismic slip is still ongoing.

Effect of heterogeneity under Tibet
In the context of crustal rheology, evidence from seismic receiver functions, geodetic inversion, resistivity and temperature profiles across the whole Tibetan Plateau show lateral heterogeneous properties (Cattin et al. 2001;Hetényi et al. 2006;Bai et al. 2010;Huang et al. 2014;Avouac et al. 2015;Sun et al.2013). The assumption of lateral uniform viscosity hardly represents the spatially variable viscous properties from the Tibetan plate to the Indian plate, and may inaccurately estimate the extent and magnitude of viscoelastic relaxation in the Himalayan region. Compared with viscoelastic relaxation of our model with that of the flat layered model (Jiang et al. 2018) (Additional file 1: Fig. S5), we find that the softer viscous lower crust over the entire region inevitably brings about an excessive estimation of viscous relaxation in Nepal. The simple heterogeneous rheological models (Zhao et al. 2017;Wang and Fialko 2018;Tian et al. 2020) incorporates the viscous lower crust layer under Tibetan plate. To investigate the difference between a fully heterogeneous model and the simple heterogeneous model, we calculated the viscoelastic relaxation with the software VISCO2.5D (Pollitz 2014), following the simple heterogeneous rheological structure and effective viscosities inferred by Zhao et al. (2017). The viscoelastic relaxation of the simple heterogeneous model exhibits a two-quadrant pattern, northward motion in the south of Fig. 8 The postseismic deformation daily displacement time series of the north and east components at GPS sites CHLM and XZZF. Black dots indicate the observed displacements. The red line indicates the data fitting of the postseismic deformation (Eq. 1). The blue line indicates the predicted viscoelastic relaxation. The green triangles indicate the sum of predicted viscoelastic relaxation and the postseismic deformation due to afterslip the rupture, southward motion in the north of the rupture, and uplift in the north and subsidence in the south (Additional file 1: Fig. S9). The opposite sign of the predicted displacements with the observations in the south of the rupture area is consistent with our 3D FEM model, consistent with other studies (Wang and Fialko et al. 2018;Jiang et al. 2018). However, it is different from the southward motion predicted by the Tian et al. (2020). Their lower viscosity of the lower crust and mantle under the Tibetan plate with the simple heterogeneous model, results in the larger postseismic displacements in Tibet than our model. Additionally, the distance of the transition zone from the MFT is chosen as 167 km, further than the model of Sun et al. (2013), leading to the smaller magnitude of the postseismic deformation in the south of the rupture area. It is difficult to explain the postseismic deformation, occurring both in Nepal and Tibet, by viscous relaxation.
The lower viscosity of the lower crust and mantle under the Tibetan plateau could account for the postseismic displacements in Tibet, but it leads to a larger discrepancy in the south of the rupture area. The opposite signs between observations and model predictions in the south of the rupture area could be reduced by assuming that a transition from strong to weak lower crust occurs farther to the north (Wang and Fialko 2018). However, this would be inconsistent with the assumption that the topographic slopes are controlled by the viscosity of material in the underlying lower crust (Clark and Royden 2000;Royden et al. 1997).
The postseismic viscoelastic relaxation obtained from geodetic measurements sheds some light on the rheology structure of Himalayan region. Wang and Fialko (2018) and Jiang et al. (2018) inferred > 10 18 Pa s and 1.6 × 10 19 Pa s for the lower crust beneath the Tibet. Zhao et al. (2017) and Tian et al. (2020), respectively, inferred a transient viscosity of 8 × 10 18 and 5 × 10 17 Pa s, a steady viscosity of 8 × 10 19 and 5 × 10 18 Pa • s for this layer. Liu et al. (2020) deduced that the effective viscosity would decrease northward from 10 18 to 10 19 Pa s around the rupture zone to ~ 3 × 10 16 -10 18 Pa s ~ 150 km north. Similar studies were conducted for several earthquakes around Tibetan Plateau. The effective viscosities on the order of magnitude of 10 18 -10 19 Pa s were predicted in studies of the 2001 Mw7.8 Kokoxili earthquake (Ryder et al. 2011;Wen et al. 2012 Kashmir earthquake (Wang and Fialko 2014) and 2008 Mw7.9 Wenchuan earthquake (Huang et al. 2014). The above studies, constraining effective viscosities by geodetic postseismic deformation, ignored the dynamic propagating stress changes from viscoelastic relaxation and those induced by any afterslip, and may somewhat estimate lower effective postseismic viscosities (Liu et al. 2020). Strain rates were highest just after the earthquakes, and the effective viscosity were equivalently lower (Liu et al. 2020). The time span is another factor which influences the inferred value of the viscosity. The rheological structure used in this paper, inferred from seismic velocities and GPS velocities, represents the long-term and steady-state viscosity, which is larger than the studies above.

Comparison with the previous afterslip models
The afterslip associated with the Gorkha earthquake is broadly investigated by several studies based on the GPS and/or InSAR data (Gualandi et al. 2016;Mencin et al. 2016;Sreejith et al. 2016;Zhao et al. 2017;Wang and Fialko 2018;Jiang et al. 2018;Jiang et al. 2019;Tian et al. 2020). These studies have found similar afterslip patterns as our study; afterslip mainly occurs downdip of the rupture area. The high afterslip patch constrained by the joint inversion of GPS and InSAR (Wang and Fialko 2015;Sreejith et al. 2016) data are more easterly than the high afterslip patch constrained by GPS alone, which might be due to the lack of GPS sites in the north of the rupture area. The moment released by afterslip is about 7.32 × 10 19 Nm over first 4.8 years, which is echoed by 5.5 × 10 19 Nm over first 1 year ( Zhao et al. 2017), 1.2 × 10 20 Nm over first 1.6 years (Jiang et al. 2018), 6.0 × 10 19 Nm and 1.2 × 10 20 Nm over first 2 years (Wang and Fialko 2018;Jiang et al. 2019), and a moment between 1.21 × 10 20 Nm and 2.37 × 10 20 Nm over 1.16 years (Liu et al. 2020). The low moment release due to afterslip, approximating 6.7% of the mainshock moment release, is also found for several earthquakes where afterslip occurred downdip of the rupture (Hsu et al. 2006;Ryder et al. 2011;Huang et al. 2014;Lienkaemper and McFarland 2017).
Different from the previous studies, we inverted the afterslip based on the heterogeneous elastic model including the topographic relief and a spherical shell. We compared the afterslip inverted by 3D FEM and flat layered models, and found the WRMS misfit between the afterslip and observations to be 3.18 mm for the FEM and 3.19 mm for the flat layered mode. The afterslip distribution of the flat layered model is somewhat smaller from that of the 3D FEM. The high value patch of afterslip inverted by flat layered model is more concentrated and maximum afterslip is up to 0.31 m (Additional file 1: Fig. S10), larger than the 0.28 m of the FEM. However, the moment released by afterslip of the flat layered model is equivalent to Mw7.15, lower than Mw7.18 of the 3D FEM.

The potential poroelastic rebound
When an earthquake occurs, sudden pore fluid pressure changes in the ambient rocks accompany the coseismic pressure changes, and then gradually evolve towards an equilibrium condition as the flow of fluid is restored, leading to time-dependent surface deformation (Peltzer et al. 1998;Wang and Kümpel 2003). It is difficult to solve the problem analytically because the deformation and pore pressure fields are coupled through the equilibrium and diffusion equation (Segall 2010). The characteristics of the anisotropic rocks further complicate the calculation. One common way is to calculate postseismic poroelastic rebound through differencing coseismic deformation models under the undrained and drained conditions, represented by variable Poisson's ratios (Fialko 2004;Barbot et al. 2010;Wang and Fialko 2018;Gonzalez-Ortega et al. 2014;Hughes et al. 2010;Hu et al. 2014). We approximately evaluated the potential poroelastic rebound without considering the time-dependent process and any heterogeneity. The poroelastic layer was assumed to be within the top 20 km of upper crust, while the undrained and drained Poisson's ratios were set to 0.28 and 0.25, respectively, and the rupture model of Tung and Masterlark (2016) was used to drive the poroelastic rebound. Our results show that the poroelastic rebound (Additional file 1: Fig. S11) spreads out from the rupture area; large poroelastic rebound concentrates the rupture area. However, the maximum displacement is on the order of millimeter, far less than the effects of afterslip and postseismic relaxation, which echoes previous studies (Zhao et al. 2017;Wang and Fialko 2018).

Insight into the seismic risk
Due to the rapid plate convergence of Indian-Eurasian collision and high strain accumulation, the Himalayan collision zone is prone to frequent devastating earthquakes. During the interseismic period, the MHT is locked from the surface to a distance of approximately 100 km downdip, corresponding to a depth of 15 to 20 km based on long time geodetic measurements (Ader et al. 2012;Stevens and Avouac 2015;Liu et al. 2016). The background seismicity along the Himalayan arc is clustered along a relatively narrow zone, which approximately coincides with the downdip end of the locked fault zone (Cattin and Avouac 2000;Bollinger et al. 2004) and the zone of greatest shear stress accumulation (Ader et al.2012). More detailed inversions for distributed interseismic coupling find that the MHT appears nearly fully locked to the south of the front of the Higher Himalaya and fully creeping to the north of it. The transition from unstable to stable slip behavior can be related to the temperature at depth (Ader et al. 2012;Stevens and Avouac 2015).
The Gorkha earthquake ruptured the deep part of the fully locked segment of the MHT Galetzka et al. 2015;Gualandi et al.2016;Elliott et al. 2016;Qiu et al. 2016), only releasing a small amount of seismic moment deficit (Feng et al. 2015;Liu et al. 2015); the inferred afterslip was concentrated at the deeper extent of the coseismic rupture, followed 2 weeks later by the Mw7.3 aftershock, which unzipped the eastern edge of the mainshock Wang and Fialko 2018), leaving the shallow portion and the west of the mainshock rupture areas still locked. Considering the interseismic fault coupling and the small amount of energy released during the coseismic and postseismic periods, the unzipped shallow portion and western segment of the MHT are still at the high seismic risk (Mencin et al. 2016;Elliott et al. 2016;Zhao et al. 2017;Wang and Fialko 2018;Tian et al. 2020;Liu et al. 2020).

Conclusion
We processed GPS displacement time series of 39 sites in Nepal and South Tibet of China, applying PCA and prior velocity constraints to extract postseismic displacements over 4.8 years after the 2015 Gorkha earthquake. The postseismic displacements show the southward and upward motion around the epicentral area. We then constructed the 3D FEM to explicitly account for the surface topography, earth curvature and heterogeneous material properties to refine the postseismic deformation after the earthquake. With the more realistic FEM, the Burgers rheology generates the viscoelastic relaxation on the order of a centimeter. The viscoelastic relaxation shows the southward motion in the north of the rupture, but nearly exhibits the opposite direction to the observations in the south of the rupture. However, the observed postseismic motions for this earthquake are dominated by afterslip. The inverted afterslip is mainly distributed downdip of the coseismic rupture with a peak of ~ 28 cm, in the area characterized by low aftershock activities. The moment release by afterslip is 7.32 × 10 19 N m, equivalent to a moment magnitude Mw ~ 7.18, approximately 6.7% of the mainshock moment release. Considering the accumulated moment deficit and the small amount energy release during and after the 2015 Gorkha earthquake, the lack of slip on a shallow portion and western segment of the MHT implies continued seismic hazard into the future.