Coseismic and Postseismic Slip of the 2005 Mw 8.6 Nias‐Simeulue Earthquake: Spatial Overlap and Localized Viscoelastic Flow

We present coseismic slip and afterslip inversion models based on the same fault geometry for the Mw 8.6 2005 Nias‐Simeulue earthquake at the Sumatran subduction zone. We estimate the coseismic slip using near‐field static GPS offsets, and vertical displacements based on satellite and coral data, while we estimate the afterslip simultaneously with viscoelastic flow using approximately nine years of GPS data following the event. With the current spatial resolution of our GPS network it is difficult to accurately resolve contributions from different postseismic mechanisms, that is, afterslip and viscoelastic relaxation from oceanic or continental mantle. We thus run many synthetic tests and models with various setups to find features that consistently appear in all our models, which we consider as robust. We find that the estimated afterslip is located primarily updip and downdip of the coseismic rupture patch and partially overlaps the updip region of the coseismic slip. We also find that the viscoelastic flow in the mantle wedge following this event was likely localized beneath the downdip region of the coseismic slip, rather than uniformly layered across the area as assumed by forward models. This localized viscoelastic flow coincides with a low‐velocity zone below Toba volcano, as imaged by tomography studies; it is possible that the viscoelastic flow beneath the volcano accelerated following this event.


Introduction
The 28 March 2005 M w 8.6 Nias-Simeulue earthquake was generated by rupture of a~400-km section of the Sunda megathrust offshore Sumatra, Indonesia (Figure 1). A portion of the megathrust in this section since then has continued slipping, as postseismic afterslip, generating surface deformation that has been recorded by the Sumatran GPS Array (SuGAr).
Afterslip often occurs in areas surrounding the coseismic rupture patch (e.g., Iinuma et al., 2016;Johanson et al., 2006;Miyazaki et al., 2004;Wang et al., 2012;Yagi & Kikuchi, 2003). While afterslip releases coseismic stress increases in some areas, it introduces time-varying stress increases to other areas. These additional stress increases can trigger aftershocks and potentially induce shallow tsunami earthquakes through loading of the shallow portion of the megathrust (e.g., Tsang et al., 2016). Understanding long-term patterns of afterslip can help us gain a comprehensive insight into the slip budget for future seismic and tsunami hazard assessment.
The coseismic slip of the Nias-Simeulue earthquake (e.g., Briggs et al., 2006;Konca et al., 2007) and its ensuing postseismic deformation have been previously investigated by a number of studies (e.g., Hsu et al., 2006;Kreemer et al., 2006). Most published postseismic studies covered less than a year of the postseismic period, with viscoelastic flow being neglected in the modeling (Hsu et al., 2006;Kreemer et al., 2006). One exception is Sun et al. (2018), who used the same length of postseismic data as we do and included viscoelastic effects from both oceanic and continental mantle. However, because their focus was the migration of the dividing boundary of seaward and landward motions, their estimates of afterslip and viscoelastic flow were obtained via trial-and-error visual inspection instead of inversion.
Globally, seafloor measurements (Watanabe et al., 2014) and postseismic modeling from many recent great earthquakes (e.g., Sun et al., 2014;Sun et al., 2018;Sun & Wang, 2015) have shown that viscoelastic flow plays an important role in postseismic deformation. Neglecting viscoelastic flow affects afterslip estimates on the megathrust, and hence interpretations of the physical properties of the fault and slip budget assessment. In addition, previous studies on the Nias-Simeulue event used different inversion techniques and fault geometries, making it difficult to accurately compare the coseismic and postseismic slip distributions. ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
Consequently, self-consistent estimates of the coseismic and postseismic slip together with viscoelastic flow are needed to draw a clear and complete picture in terms of the slip budget over the seismic cycle.
In this paper, we remodel the coseismic and postseismic deformation following the 2005 M w 8.6 Nias-Simeulue earthquake with a consistent Slab1.0-based fault geometry that is discretized into triangular patches. We incorporate viscoelastic flow estimation simultaneously into the afterslip inversion by inverting strains within finite volume cuboids Moore et al., 2017) distributed in the continental mantle wedge, using a new technique as described in Qiu et al. (2018). We use self-consistently estimated GPS coseismic offsets and postseismic time series spanning approximately nine years .
Our results show that the coseismic slip and afterslip are generally spatially complementary, with partial overlap updip of the estimated coseismic slip patch. With viscoelastic flow incorporated into our models, the results suggest that our limited SuGAr stations at the shallow depth cannot resolve the contributions from shallow afterslip and oceanic mantle, but are sufficient to resolve a certain degree of contributions from deep afterslip and the continental mantle wedge. Our results also show that simple layered viscoelastic forward models predict uplift at GPS station SAMP, in contrast to the subsidence that was recorded. The subsidence can be explained by localized, deep viscoelastic flow below the downdip area of coseismic slip in the continental mantle wedge.

Geodetic Data and Inversion Methods
In this section, we detail our data and methods for the coseismic and postseismic slip inversion, including fault geometry, the joint inversion of both afterslip and viscous flow, and the discretization strategy of the Nias-Simeulue earthquake. The brown beach ball represents the gCMT focal mechanism for the main shock. Vectors indicate the coseismic displacements, and postseismic displacements over different time periods (note that different vectors have different scales). White triangles and pink circles indicate the locations of SuGAr GPS stations, and the locations of coseismic vertical displacements estimated from ASTER images and coral data, respectively. Dashed brown box outlines the region for seismicity shown in Figure S10 (International Seismological Centre, 2015;NCEDC, 2016). Green dots show the background seismicity from the ANSS catalog that spans from 1 January 2000 to the day before the M w 9.2 2004 Sumatra-Andaman earthquake. The GPS data error ellipses show 95% confidence level. The gray thick arrow shows the plate convergence rate between the Indo-Australian and Sunda plates (DeMets et al., 2010). oceanic mantle and continental mantle wedge using 3-D cuboid volumes. We also describe synthetic and resolution tests in which we investigated the resolvability of oceanic and continental mantle, and potential trade-offs in the modeling of afterslip and viscoelastic flow.
2.1. Coseismic Modeling 2.1.1. Coseismic Deformation We used GPS, satellite, and coral data to constrain the coseismic deformation and invert for the coseismic slip distribution. The coseismic offsets estimated from SuGAr daily position time series are taken from Feng et al. (2015) (Table S1 and Figure 1). The coseismic vertical displacements estimated from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images (Table S2) and coral microatolls (Table S3) are compiled in Briggs et al. (2006).

Coseismic Slip Inversion
Compared with previous studies that used either a planar fault geometry (Briggs et al., 2006;Konca et al., 2007), two planes with different dip angles (Hsu et al., 2006), or a curved geometry based on slab contours (Kreemer et al., 2006;Linette et al., 2010), we created a curved fault mesh based on Slab1.0 (Hayes et al., 2012) using triangular elements. The Slab1.0 geometry agrees well with the relocated seismicity recorded by a temporary seismic array ( Figure S1) (Tilmann et al., 2010). We chose a maximum depth of~80 km, slightly deeper than 30-60 km which corresponds to the 350-450°C isotherm based on thermal modeling along a profile at Simeulue (Hippchen & Hyndman, 2008). Temperatures between 350°C and 450°C have been suggested to indicate the downdip limit of coseismic rupture, where velocity-weakening (seismic rupture) transitions into velocity strengthening (aseismic creep; e.g., Wang, 2007;Wang et al., 2003). Our final fault mesh is divided into 400 equilateral triangular fault patches, with the length of all three sides being 30 km.
We used a linear least squares method to estimate coseismic slip. We constrained the rake to 118+/−45°b ased on the rake given in the global centroid moment tensor catalog (Ekström et al., 2012), and applied a Laplacian constraint to spatially smooth the slip. We calculated Green's functions for triangular elements in a half-space using the analytic solution from Nikkhoo and Walter (2015). Green's functions for triangular elements in a layered elastic crust have not been developed yet. Although previous studies show that a layered Earth structure affects the magnitude of slip and moment estimates, layered structures produce a spatial pattern of slip estimates similar to homogeneous half-space Earth models (e.g., Briggs et al., 2006;Hsu et al., 2006;Konca et al., 2007;Kreemer et al., 2006). Using the half-space triangular Green's functions is thus most suitable for our goal of investigating the spatial extent of coseismic and postseismic slip estimates on a curved fault. We tested different relative weightings between GPS offsets and the coseismic vertical displacements from ASTER and coral. We find weighting GPS offsets 3.5 times more produced the best model results that fit both data sets reasonably well within their uncertainties.
We used a Bootstrap method (e.g., Feng et al., 2016;Parsons et al., 2006;Wright et al., 2003) to quantify uncertainties in our estimated coseismic slip model. This technique uses a Monte Carlo (MC)-based resampling scheme (see Text S1; Feng et al., , 2014. We penalized the fault boundary triangular elements and conducted two sets of experiments using the MC method. For the first set, we simulated the uncertainties of geodetic data (GPS offsets, and ASTER and coral vertical displacements) and varied the spatial smoothing, while for the second set we only simulated the uncertainties of geodetic data. We performed a total of 100 iterations for each set . The mean and uncertainties of the coseismic slip estimates from the MC simulations are presented in Figure 2 for the first set and Figure S2 for the second set. Comparing results from the two sets, we find that the estimated uncertainty of slip is dominated by the geodetic data rather than the spatial smoothing.

Postseismic Modeling 2.2.1. Postseismic GPS Time Series
We obtained approximately nine-year-long postseismic GPS time series following the 2005 M w 8.6 Nias-Simeulue earthquake, for 13 stations, from Feng et al. (2015). The 13 GPS stations comprise 12 SuGAr stations and 1 station (SAMP) operated by the Indonesian Geospatial Information Agency (BIG, formerly known as BAKOSURTANAL; Figure 1).
In addition to the 2005 event, the stations used in this study were affected by the postseismic deformation of five other earthquakes including the 2004 M w 9.2 Sumatra-Andaman earthquake (e.g., Meltzner et al., 2006), the 2008 M w 7.4 Simeulue earthquake (Morgan et al., 2017), the 2010 M w 7.8 Banyak Islands earthquake (Morgan et al., 2015), the 2010 M w 7.2 Simeulue earthquake (Wong et al., 2017), and the 2012 M w 8.6 Wharton Basin sequence (e.g., Hill et al., 2015). In order to obtain the deformation signal resulting only from the 2005 Nias-Simeulue event, we fit the original time series with a nonlinear functional model that contains offsets and logarithmic functions for all recorded earthquakes, along with long-term rates, annual and semiannual cycles. In other words, the coseismic offsets and postseismic decays for all events were estimated simultaneously in a self-consistent manner. Full details of the GPS processing and analysis techniques are given in Feng et al. (2015). We filled data gaps with these model fits, and used the gap-filled time series for our joint inversion, which is required for achieving stable inversion results when we have a sparse GPS network and many model parameters. We removed from the original time series the postseismic deformation associated with four of the five events excluding the 2004 event. Unfortunately the 2004 event occurred too close to the 2005 event in time (approximately three months apart) to allow the separation of the 2004 and 2005 postseismic signals . To give some idea of the impact on our results of assuming that all the deformation was from the 2005 event, we calculated the cumulative displacements over the same data period from a simple postseismic model of the 2004 event. The resulting trench-parallel motions in the 2005 rupture region are <10 cm as compared with >50-cm trench-normal motions observed within the same time period ( Figure S3; Banerjee et al., 2007;Pollitz, 1997;Pollitz et al., 2006Pollitz et al., , 2008. We also compared the cumulative postseismic displacements with and without correcting the 2004 event, and found that the overall pattern of the 2005 postseismic deformation does not change much ( Figure S3b). In order not to introduce errors due to mismodeling of the 2004 event, we chose not to remove the postseismic effect of the 2004 event using a model.

Simultaneous Inversion of Afterslip and Viscous Strain
Possible postseismic mechanisms include afterslip, poroelastic rebound, and viscoelastic relaxation. Poroelastic rebound is usually confined to the areas adjacent to the rupture and lasts for less than a few months (Hughes et al., 2010): in most cases, this process contributes significantly less surface deformation over long time periods and across large spatial scales compared to the other two processes, and thus, we ignored it in this study.
Traditionally, afterslip has been modeled with GPS time series that are "corrected" for viscoelastic deformation predicted from forward viscoelastic models. The forward predictions rely on having an accurate coseismic source model and precise knowledge of the rheological structure of the mantle. Given that these are not usually available, forward predictions might introduce large uncertainties into afterslip estimates. As an alternative to this traditional method, Qiu et al. (2018) proposed an inversion technique that simultaneously inverts for afterslip and viscous strain in the mantle, without a priori assumptions about the mantle rheology. They demonstrated that this inversion approach is capable of estimating temporal and spatial heterogeneity in the viscosity structure of the mantle wedge. The novelty of this technique is that it can intrinsically take into account the mechanical coupling of afterslip and viscoelastic relaxation, and most importantly without a priori assumptions about the rheology of the mantle wedge or coseismic source, lessening the possible uncertainties and trade-offs introduced by forward models. We applied this joint inversion approach to model the postseismic deformation of the 2005 event.
The joint inversion technique implements a modified version of the Extended Network Inversion Filter (Hsu et al., 2006;McGuire & Segall, 2003;Miyazaki et al., 2004;Qiu et al., 2018;Segall & Matthews, 1997). We used the Extended Network Inversion Filter to model the time evolution of both afterslip on the megathrust (using the triangular Green's functions) and viscous strains in the mantle wedge (using the Green's functions for strains within finite volume cuboids) based on the approximately nine-year postseismic GPS time series. We ran the filter with a 10-day sampling rate to reduce the computational burden. Our selected sampling rate might have some effect on the estimation of transient slip rates or strain rates in the early time stages, but it should not significantly change the estimates of the final cumulative afterslip on the megathrust and viscous strain in the mantle wedge. We allowed all fault patches to slip, and applied spatial smoothing regularization for both slip rates and strain rates to avoid overfitting the data. We also penalized isotropic strain in the inversion, and constrained strain directions to be perpendicular to the induced coseismic stress changes at each strain component Qiu et al., 2018;Tsang et al., 2016).
We discretized the lower crust and upper mantle beneath the coseismic rupture patch into three layers of cuboids that extend from the bottom of the fault to a depth of~180 km ( Figure S4). Each layer consists of three cuboids that are 570 km long along strike, 60 km wide along dip, and 40 km tall along depth ( Figure S4). We used elongated cuboids that extend the whole length of the Nias-Simeulue rupture area instead of shorter cuboids, as we have only two stations PBLI and SAMP directly above the cuboids (Figures 1 and S4). While each triangular patch has two components of slip estimates, each cuboid has six independent strain components; thus, we do not have enough resolution to use smaller cuboids. We calculated the resolution matrix for this model configuration, which shows the thrust motion-related components e11, e13, and e33, are well resolved ( Figure S4). This configuration of cuboids is our preferred model setup for the mantle; however, it does not include the oceanic mantle. In the next section, we show our tests with oceanic mantle added to the model and explain why the model setup without oceanic mantle is chosen.

Model Setup and Resolution Tests for Postseismic Modeling
In comparison to the large number of data points we have for coseismic modeling, we have only 13 GPS stations for postseismic modeling. So we face a challenge of resolving afterslip and viscoelastic deformation using a limited number of data points. In order to understand what we can or cannot resolve using the available data, we conducted three sets of tests. The first set including nine different model setups for the mantle with only synthetic viscoelastic inputs focuses on testing whether we can resolve contributions from both the oceanic mantle and continental mantle wedge (section 2.3.1), the second set using the same nine model setups but with real data is used to determine our preferred model setup and identify which features of afterslip are robust regardless of mantle discretization (section 2.3.2), and the last set tests the trade-offs between afterslip and viscous strain using our preferred model setup for the mantle (section 2.3.3).

Can We Resolve Viscoelastic Deformation From Both the Oceanic Mantle and Continental
Mantle Wedge? While afterslip always produces seaward motions after great subduction earthquakes, viscoelastic deformation following great subduction earthquakes is characterized by landward motions in the fore-arc region and seaward motions in the back-arc region, with the dividing boundary migrating landward with time (Sun et al., 2014;Watanabe et al., 2014). Landward motions are contributed by the oceanic mantle, whereas seaward motions are contributed by the continental mantle wedge (Barbot, 2018). Therefore, in order to model 10.1029/2018JB017263 Journal of Geophysical Research: Solid Earth viscoelastic relaxation correctly, ideally we should include both the oceanic mantle and continental mantle wedge in our models.
However, the question is whether the SuGAr network can resolve simultaneously viscoelastic deformations of the landward motions from the oceanic mantle and seaward motions from the continental mantle wedge while resolving afterslip. To address this question, we tested nine model setups that use various different cuboid sizes and include the continental mantle only, the oceanic mantle only, and both (see Text S2 for more information). We discretized the oceanic mantle from~60-to~220-km depth (similar to previous studies, e.g., Han et al., 2008;Pollitz et al., 2006Pollitz et al., , 2008 into elongated cuboids with a length of 570 km and varying widths and thicknesses.
We first created synthetic time series for the SuGAr stations that represent the pure viscoelastic contribution using Relax (Barbot & Fialko, 2010) based on our preferred coseismic slip model, with a Maxwell rheology (viscosity 10 19 Pa s) under an elastic lid of~60 km (e.g., Han et al., 2008;Pollitz et al., 2006Pollitz et al., , 2008. We then inverted the viscoelastically "induced" time series for both afterslip on the fault and viscous strain in the cuboids to investigate the possible trade-offs between afterslip and viscous strain. Details about these tests are given in Text S2. The viscoelastic forward predictions and their corresponding inversion results are shown in Figure S5 and Table S4. Regardless of the model setup for the mantle, the results always include a significant contribution from afterslip (gray vectors in Figure S5) even though the synthetic data are generated by viscous flow only. This means that the viscoelastic deformation cannot be fully recovered because much of it is mapped to afterslip on the fault, leading to significant trade-offs between afterslip and viscous flow. Due to lack of geodetic measurements in the fore-arc region, particularly near the trench, it is difficult to disentangle contributions from each part of the mantle and afterslip on the fault.

Preferred Model Setup for the Mantle
In addition, we used the same nine model setups mentioned above to invert the GPS postseismic displacement time series ( Figure S6). The results show that none of these model setups can fully predict both the landward and seaward motions. Although landward motions from the oceanic mantle are not resolved, the seaward motion and subsidence at SAMP (green vectors) are correctly predicted by almost all the model setups except model setup OMCM-1, suggesting that the SuGAr network is sufficient to resolve the contribution of the continental mantle wedge to a certain degree ( Figure S6). Considering the degree of the trade-offs and model fits with the data, we chose to use model setup CM-1 ( Figure S6) that contains only the continental mantle wedge as our preferred model setup (Text S2). We note that ignoring the oceanic mantle would underestimate shallow afterslip. We use model setup CM-1 for further checkerboard and resolution tests, and our final postseismic modeling.

Synthetic Tests Using the Preferred Model Setup
We conducted 10 synthetic tests ( Figure S7) to understand our model resolution and investigate trade-offs between afterslip and viscous flow in the continental mantle wedge. These tests include trade-offs between afterslip and viscous strain in the continental mantle wedge (Tests 1-3), checkerboard (Tests 4 and 5), resolution of the shallow and deep afterslip (Tests 6-8), sensitivity to input strain magnitude (Tests 9 and 10). We generated synthetic displacements that were calculated from either zero or unit slip on the triangular subfault patches. Following Qiu et al. (2018), we used the spatial distribution of the coseismic stress changes of the Nias-Simeulue earthquake as the synthetic strains, with their magnitudes scaled by the amplitude of the second invariant of coseismic stress. We used our coseismic slip model to calculate the stress changes (Nikkhoo & Walter, 2015) on these cuboids. To optimize the afterslip and viscous strain solution, we obtained the spatial and temporal smoothing parameters for the Kalman filter in an iterative way until they converged (Qiu et al., 2018). We implemented these optimized smoothing parameters in all the experimental tests, and also the postseismic modeling of the 2005 event.
Our resolution tests show that the SuGAr network can resolve the synthetic slip location at places where GPS stations are available, but slip amplitudes are underestimated by~20% to 50% depending on the location and the dimension of the synthetic slip patches (Tests 4-8 in Figure S7). When afterslip is underestimated, the result is overestimation of viscous strain estimates (Tests 2-10 shown in Figure S7). We observe the overestimation of viscous strain in all tests. These trade-offs affect the magnitude estimates, but do not affect the spatial pattern of either slip or viscous strain; thus, the spatial pattern of existence of both the shallow and 10.1029/2018JB017263 Journal of Geophysical Research: Solid Earth deep slip seem to be a robust feature of our models (Tests 6-10 in Figure S7). Details of all these tests are given in Text S2.

Results
Our estimated coseismic slip model shows two areas of high slip on the megathrust, one underneath Nias and the other under southern Simeulue, with an area of low slip between them (Figure 2a). The estimated afterslip distribution of the Nias event, in general, spatially complements the coseismic rupture, located primarily updip and downdip of the coseismic slip, with updip afterslip partially overlapping the coseismic slip. Our model suggests that localized rather than broad viscoelastic flow downdip of the coseismic rupture area is required to explain the observed subsidence at mainland station SAMP.

Coseismic Slip
We determine the final coseismic slip estimates (Figure 2a) from the mean of the MC simulations. The coseismic slip distribution has two high-slip regions and a low-slip region in-between, forming a saddleshaped rupture area (Figure 2a). The slip maxima for the two high-slip regions are located directly under the northern tip of Nias (~14 m) and downdip of the southern tip of Simeulue (~10 m; Figure 2a). The low-slip region has reduced slip (<5 m).
The uncertainties of the slip estimates are controlled by the spatial configuration of the SuGAr network, and ASTER and coral points (Figures 2b and S2). The slip uncertainties for the main rupture area are <50 cm, suggesting that the overall feature of the main rupture is relatively well constrained. For areas north of the Banyak Islands, and east and south of Nias, the slip uncertainties are >1 m, suggesting that the slip in these areas is not well constrained, possibly due to the large uncertainties of the coral measurements and the paucity of GPS stations. However, as the areas that have large slip uncertainties are not in the main slip area of slip >2 m, they do not affect the robustness of our main slip estimate.
In order to test whether the saddle shape is a robust feature of the coseismic rupture, we conducted two checkerboard tests: one with uniform slip extending from northern Simeulue to southern Nias and the other with two high-slip regions underneath Simeulue and Nias, respectively ( Figure S8). The results of the two tests show that the spatial coverage of our data points is dense enough to differentiate between the two high-slip regions and one uniform-slip region. This further suggests that the spatial pattern of the coseismic rupture is a robust feature. Overall, the spatial pattern of our coseismic slip model is consistent with published slip models (e.g., Briggs et al., 2006;Hsu et al., 2006;Konca et al., 2007;Kreemer et al., 2006). Our coseismic slip model suggests a coseismic moment release of M w~8 .6, assuming a shear modulus of 30 GPa.

Afterslip
The model setup and resolution tests in section 2.3 indicate that the SuGAr network can resolve only the role of deep afterslip and the continental mantle wedge, with motions at SAMP correctly predicted. Shallow afterslip is largely underestimated and deep afterslip is slightly underestimated due to the missing landward motions that would be expected from the oceanic mantle ( Figure S6, CM-1). However, the time evolutions of the updip and downdip afterslip is persistent regardless of the mantle configuration ( Figure S9), so we conclude that the spatiotemporal pattern of the afterslip is a robust feature. We show the cumulative slip pattern during four time periods in Figure 3. Within the first six months (Figure 3a), the afterslip occurred mainly offshore of Nias and Simeulue Islands, with some occurring southeast of Nias Island. In six months to the first year, the afterslip continued to occur offshore of Nias and Simeulue. In the first year following the event (Figures 3a and 3b), the downdip afterslip was significantly lower than the shallow afterslip, suggesting that the slip rate of the shallow megathrust was faster than that of the deep portion. This slip rate difference would be even larger if we include the oceanic mantle ( Figure S6). The majority of the deep afterslip downdip of the coseismic rupture occurred after a year following the earthquake (Figures 3c and 3d). Our estimated spatial pattern of the afterslip is consistent with previous studies (e.g., Hsu et al., 2006;Kreemer et al., 2006;Linette et al., 2010); however, our obtained temporal evolution of the afterslip has not been revealed by any previous studies.
The estimated afterslip spatial pattern correlates well with the aftershock distribution in the same period ( Figure 3). We also show the time series of the depth distribution of the seismicity from ANSS, USGS, and 10.1029/2018JB017263 Journal of Geophysical Research: Solid Earth ISC catalogs in Figure S10. In the year following the earthquake, we observe a seismically quiet period at depths of~55 to 100 km where the deep afterslip occurred ( Figure S10). Within the same time period and at the same depth range, we estimate low levels of afterslip values (<25 cm), suggesting that spatiotemporal changes of the afterslip may drive the onset of seismicity in the fault zone. We would need better resolution to confirm whether there is a genuine downdip migration of the afterslip, or whether this is a trade-off with the viscoelastic flow. Overall, however, these observations suggest that the associated aftershocks were primarily driven by the afterslip process (e.g., Barbot et al., 2009;Hsu et al., 2006;Perfettini & Avouac, 2004;Savage & Yu, 2007).
The model-predicted postseismic decays fit the time series fairly well ( Figure S11). Especially for station SAMP (Figure S11), the vertical component is poorly explained by the afterslip-only model and better explained by adding viscoelastic flow ( Figure S12, and this is further discussed in section 4).

Afterslip and Its Relationship With Coseismic Slip
Overall, the afterslip spatially complements the coseismic slip (Figure 4a). However, a spatial overlap between coseismic slip and afterslip may exist offshore of the two high-slip regions close to BSIM and LHWA (Figure 4a). Resolution tests ( Figure S7), and also error analysis (Figure 2b), suggest that the slip amplitude might be underdetermined, but the spatial locations of coseismic slip and afterslip are

10.1029/2018JB017263
Journal of Geophysical Research: Solid Earth reasonably well determined at these locations. Coseismic uplifts at BSIM and LHWA indicate that the slip must have occurred beneath Nias and Simeulue on the megathrust, but fail to pinpoint how much further it had extended offshore due to the lack of seafloor measurements. These sudden uplifts were reversed to subsidence for the postseismic period following the event. The vertical reversals at these stations indicate that the afterslip occurred on the megathrust offshore of Nias and Simeulue. We note that the estimated uncertainty of the coseismic slip in the main rupture area is <50 cm. If we use the 2-m coseismic slip contour to delineate the main rupture area (the dashed purple line in Figure 4a), the area of overlap offshore Simeulue is reduced, but the area of overlap offshore Nias remains large (>30 km in length along dip; Figure 4a). The location of the coseismic slip can be affected by some factors, such as the geometry of the fault. The dip of the shallow megathrust is poorly determined due to less accurate or missing seismicity for determining the location of the plate interface (Hayes et al., 2012). However, such a large spatial overlap is unlikely to be introduced by uncertainties in the fault geometry: in the case of Nias, to make the overlap disappear, the dip angle must vary by >50°, assuming a shallow 10-degree dipping fault. With these lines of evidence, we conclude that the spatial overlap offshore Nias is a robust feature.

Inferred Localized Viscoelastic Flow
The estimated viscous strains for the cuboids, inferred primarily across the boxes of the first layer (Figure 4b), decrease both horizontally away from the megathrust and with depth. From the top layer to the third layer, the model resolution for viscous cuboids decreases with increasing distance to SAMP, resulting in less reliable model estimates for the cuboids in both the second and third layers ( Figures S4 and S13). Except cuboids 1, 4, and 7 that are relatively far from SAMP ( Figure S13), the motions of the rest of the cuboids are well resolved and have similar orientation to the coseismic stress changes. Based on the resolution matrix of the top layer of cuboids ( Figure S4, thrust motion-related strain components e11, e22, e13, and e33), and also the resolution tests ( Figure S7), the viscous strain magnitude estimates suffer from trade-offs with the afterslip estimates. Nevertheless, the spatial location of the strain is robust.
The afterslip can generate <20% of the observed subsidence at SAMP (Figures S11 and S12), and it cannot be deeper than our estimated location as deeper afterslip would predict uplift at SAMP. In addition, the chemical and temperature environment at this depth favors ductile flow. Consequently, we conclude that the vertical subsidence at SAMP is caused predominantly by viscoelastic flow, but has to be localized as shown in our inversion results (Figure 4b).

Physics-Based Forward Models Predict Opposite Vertical Motions Recorded at SAMP
To understand the general spatial pattern of the viscoelastic deformation following this event, with focusing on the vertical component of SAMP, we additionally created a suite of physics-based coseismic stressdriven forward viscoelastic models. As most existing programs for viscoelastic modeling (VISCO1D (Pollitz, 1997), PSGRN (Wang et al., 2006), Relax (Barbot & Fialko, 2010)) do not use triangular slip patches, we thus combined the slip estimates of neighboring triangular patches to form rectangular patches that maintain the same seismic moment. We used the combined rectangular slip model to predict surface displacements from viscoelastic relaxation. We ran a suite of forward models using Relax (Barbot The open and red triangles represent the relative locations of GPS stations PBLI and SAMP, and Toba volcano, respectively. The pink shaded-transparent area represents the approximate location of an inferred lower velocity zone from tomography results of . The isolated single cuboid shows the thrust motion related shear strain component e13. The error ellipses for the data are shown with 95% confidence. As and Vf stand for afterslip contribution and viscoelastic contribution to the surface displacements, respectively. & Fialko, 2010), with two end-member Maxwell rheology models based on the minimum and maximum viscosities estimated by various studies of the neighboring M w 9.2 2004 Sumatra-Andaman earthquake (e.g., Broerse et al., 2015;Han et al., 2008;Hoechner et al., 2011;Hu & Wang, 2012;Pollitz et al., 2006;Qiu et al., 2018). We also tested different elastic thicknesses based on previous studies (e.g., Broerse et al., 2015;Han et al., 2008;Hu et al., 2016;Pollitz et al., 2006;Wiseman et al., 2015).
For all the forward models we obtain postseismic uplift for mainland Sumatra, regardless of the elastic thickness and viscosity ( Figure 5). A more complex FEM postseismic model for the Nias earthquake (Wiseman et al., 2015), which included a subducting slab and Maxwell rheology of the upper mantle, also shows uplift motion at SAMP. On the other hand, the GPS vertical time series for SAMP show continuous subsidence following this event, over the entire approximately nine years (Figures 4a and S11). Therefore, forward models that use a simple layered rheology over a broad spatial scale cannot explain the recorded subsidence at SAMP, suggesting that another mechanism (e.g., afterslip or local heterogeneity in the rheological structure) may control this subsidence. We will discuss this in more detail in section 4.2. Spatial complementarity between coseismic slip and afterslip has been observed globally at many subduction zones (e.g., Iinuma et al., 2016;Johanson et al., 2006;Miura et al., 2006;Miyazaki et al., 2004;Wang et al., 2012), suggesting different physical properties for coseismic patches that favor stick-slip and afterslip regions that prefer aseismic creep. However, a number of recent studies have suggested spatial overlap between coseismic slip and afterslip (e.g., Barnhart et al., 2016;Bedford et al., 2013;Bürgmann et al., 2002;Johnson et al., 2012;Pritchard & Simons, 2006;Salman et al., 2017;Tsang et al., 2016). Our results suggest the coseismic slip and afterslip also overlap for the Nias earthquake. Factors that control the rupture characteristics of an earthquake include the frictional properties on the fault, chemical and thermal boundaries (Lay et al., 2012), and subducting barriers (Meltzner et al., 2012). In addition, dynamic modeling (Helmstetter & Shaw, 2009) suggests that frictional differences are not necessary to fit postseismic deformation; instead, a large range of rate and state model parameters can explain the postseismic observations.
The frictional properties that control the shallow part of the megathrust are poorly understood (e.g., Lay & Schwartz, 2004;Wang & Timothy, 2004). The shallow megathrust has been described as velocity strengthening, arresting fast seismic rupture, due to the presence of unconsolidated or semiconsolidated sediments (e.g., Byrne et al., 1988;Marone & Scholz, 1988;Wu et al., 1975). According to this theory, the shallow segment would be persistently aseismic. The transition between deep velocity-weakening and shallow velocitystrengthening regions is likely marked by the seismic front (Byrne et al., 1988) that is located along the edge of the estimated updip afterslip patches (Figures 3 and 4). Seismicity at the Nias-Simeulue segment is concentrated around depths of~20 km, at the downdip edge of our estimated shallow afterslip patch ( Figure 3). However, the shallow megathrust can also be seismogenic (Hubbard et al., 2015), as is evidenced by tsunami earthquakes such as the 2010 M w 7.8 Mentawai earthquake (e.g., Hill et al., 2012) and the 1994 M w 7.6 and 2006 M w 7.8 Java events (e.g., Abercrombie et al., 2001;Ammon et al., 2006). Here at the Nias segment, the shallow part also ruptured during the 1907 M w 7.8 tsunami earthquake (Kanamori et al., 2010) (Figure 6). Clearly, both fast and slow slip (e.g., Saffer & Wallace, 2015;Wallace et al., 2016Wallace et al., , 2017 can rupture the shallow portion of the fault, and essentially, these rupture styles cannot be explained by fixed frictional properties. Possible complex ambient conditions of the fault zones such as fluid pressure (Moreno et al., 2014) and material complexity from fault gouge of mature fault zones (Hyndman et al., 1997) allow dynamic weakening mechanisms (Noda & Lapusta, 2013) to promote coseismic rupture extension to the shallow weak segment without reaching to the trench. This may apply to the coseismic slip and afterslip overlap of the Nias event and possibly to other large and great earthquakes (e.g., Barnhart et al., 2016;Bedford et al., 2013;Bürgmann et al., 2002;Johnson et al., 2012;Pritchard & Simons, 2006;Salman et al., 2017;Tsang et al., 2016). The stress loading of the shallow segment from the Nias rupture together with residual stress left over from previous events within seismic cycles may trigger the release of the cumulative strain (Rafael et al., 2018) as afterslip occurring after the 2005 event ( Figure 6), or occasionally, the cumulative strain could be released as tsunami earthquakes, as evidenced by the 1907 M w 7.8 event ( Figure 6).

Localized Viscoelastic Flow Accelerated by the Earthquake
Different tectonic settings correspond to different physical conditions of mantle rocks and therefore different behaviors in their underlying mantle. Physical conditions also vary along the length of a subduction zone. Figure 7b shows the vertical displacement time series for station UMLH, which was affected primarily by the 2004 M w 9.2 Sumatra-Andaman earthquake in the north, and for LNNG, MKMK, and JMBI, which were affected by the 2007 M w 8.4 Bengkulu event (Figure 7a). These stations show vertical uplift, which is successfully explained by both inversion that incorporates viscoelastic relaxation (Qiu et al., 2018) and layered forward models (Wiseman et al., 2015). In comparison, station SAMP (Figure 7b),~400 km away from UMLH (Figure 7a), shows vertical subsidence, which cannot be explained by the broad viscoelastic flow from any forward models ( Figure 5). These spatially opposite vertical motions may suggest varying viscosity of the mantle wedge along the length of the Sumatran subduction zone. Close to SAMP, a region with low S wave velocities and high V p /V s ratios was imaged by seismic tomography beneath Toba volcano (e.g., Figure 4b), suggesting the presence of a localized high temperature zone due to the presence of a magma reservoir. A high-temperature anomaly at 50-km depth was also suggested beneath SAMP ( (Chlieh et al., 2007), 2005 M w 8.6 Nias-Simeulue earthquake (Konca et al., 2007), 2007 M w 8.4 Bengkulu earthquakes (Konca et al., 2008), and the 2010 M w 7.8 Mentawai tsunami earthquake (Hill et al., 2012). Shaded areas represent the published coseismic rupture areas. Triangles represent the SuGAr network. Moment tensors show the main shock of each event from gCMT. Red triangle shows the locations of the example GPS sites. (b) Vertical time series of selected GPS stations, with the same color code as the coseismic rupture patches. Vertical dotted lines show the times of the Nias and Bengkulu earthquakes.

10.1029/2018JB017263
Journal of Geophysical Research: Solid Earth dehydrated water from the subducting slab migerates upward to the the magma reservoir beneath Toba. High temperature together with excess fluid may lower the viscosity and produce localized viscoelastic flow; this may have caused the subsidence of SAMP. An earthquake-triggered viscoelastic flow beneath volcano chains was also observed following the 2010 M w 8.8 Maule earthquake at Chile subduction zone (Li et al., 2017).

Conclusions
We remodel, in a self-consistent manner, the coseismic slip and long-term postseismic deformation following the 2005 M w 8.6 Nias-Simeulue earthquake. Our results suggest that, on a large spatial scale, the estimated coseismic slip distribution and afterslip are spatially complementary; however, locally, they partially overlap offshore of the Simeulue and Nias Islands. The deep afterslip has limited slip within the first year following the earthquake, which may explain missing seismicity at similar depths. The correlation between the afterslip and seismicity suggests that the associated aftershocks are primarily driven by spatiotemporal changes in the afterslip process. Localized instead of broad viscoelastic flow in the mantle wedge is needed to explain subsidence observed in mainland Sumatra, which cannot be explained by layered forward models.