*5.2. Aftershocks*

We consider both aftershock sequences as composite events and we obtain elastic dislocation models for aftershock sequence 1 and 2 so that their contributions can be removed from the post-seismic displacement fields. We do not consider the 25 August 2018 aftershock located east of the epicenter because there is no clear signal associated with it in ascending and descending data (Figure S1) and no geodetic focal mechanism can be estimated. For both ascending and descending orbits, we use interfergrams reconstructed from the raw time series. They are less noisy than the observed interferograms because the network inversion of interferograms filters out temporal decorrelation noise.

For aftershock sequence 1, the best-fitting dislocation has a length and width of 12.9 km and 5.4 km and 0.4 m slip. For aftershock sequence 2, the best-fitting dislocation has a length and width of 12.7 km and 8.8 km and 1.8 m slip (Table 2). The observations, modelled displacement and difference between observations and model are shown in Figure 3a,b. The red lines are nearby faults. The geodetic focal mechanisms (shown in Figure 3) are very similar to the Centroid-Moment-Tensor (CMT) solutions (shown in Figure 1).

#### *5.3. Afterslip for Time Period 1*

To test whether the surface displacements are consistent with afterslip, we use the ascending and descending data until June 2018 (Figure 4a) after removal of the contribution from the January 2018 aftershock sequence. We estimate the best-fitting uniform dislocation using a homogeneous elastic half space model. The best-fitting solution (Table 2, black solid rectangle in Figure 4a) shows that the observations are well explained by 0.8 m of uniform slip along one dislocation with a length and width of 45.5 and 6.4 km located on the near-horizontal, western, up-dip continuation of the coseismic fault.

**Figure 3.** Uniform slip inversion result of aftershock sequence 1, aftershock sequence 2 for ascending data and descending data. (**a**) Aftershock sequence 1, ascending data from 10 January 2018 to 16 January 2018, descending data from 5 January 2018 to 17 January 2018; (**b**) aftershock sequence 2, ascending data from 18 November 2018 to 30 November 2018, descending data from 25 November 2018 to 7 December 2018. Black solid rectangle: best-fitting dislocation. Red lines: nearby faults. Black dot: upper edge of fault. The beach balls for aftershock sequence 1 and 2 are obtained based on the inverted fault parameters. Blue points for aftershock sequence 1 and 2: reference point.


**Table 2.** Range, optimal and uncertainties of InSAR-inversion fault slip parameters for the aftershock sequence that occurred on 11 January 2018, and the aftershock sequence that occurred on 25 November 2018 and afterslip for time period 1.

a maximum posteriori probability solutions of 2.5% and 97.5% from the posterior probability density function of fault parameters. b The positive value means the hanging wall showed relative uplift; c the negative value means that the hanging wall direction of the motion was opposite to the strike; d longitude and latitude of the corner of upper edge; e calculated based on the optimal inversion parameters.

**Figure 4.** Best-fitting modelling results for (**a**) afterslip period (until June 2018) and (**b**) viscoelastic relaxation period (from June 2018 to February 2020). Black thin dashed rectangle: coseismic fault. Black solid rectangle: best-fitting afterslip dislocation. Black dot: upper edges of fault. Red lines: nearby faults. The reference point is shown in Figure 2.

#### *5.4. Viscoelastic Relaxation for Time Period 2*

We assume that the displacements in time period 2 are caused by viscoelastic relaxation. In order to constrain the rheological structure we use a 25 km-thick elastic upper crust and a 20 km-thick viscoelastic Maxwell lower crust over a viscoelastic Maxwell upper mantle (Figure 5a; the average depth of the regional Moho is 45 km according to the CRUST 1.0 model) [35,36]. In the grid search the viscosities of the lower crust and upper mantle range from 1 × 10<sup>17</sup> Pas to 1 × 10<sup>21</sup> Pas; we use a step size for the logarithm of viscosity of 0.25.

Our result (Figure 5b) yields the best estimates of the parameters η*lc* = 1+0.8 −0.4 × 10<sup>19</sup> Pas, where the confidence interval is calculated using F–Test [37–39]. Our data resolve the viscosity of the lower crust but not the viscosity of the uppermost mantle, except that it is larger than 0.6 × 10<sup>19</sup> Pas.

The comparison between the observed and simulated displacement (Figure 4b) shows that the model explains the descending data well but there is a discrepancy for the ascending data. This discrepancy can be attributed to the larger tropospheric noise for the ascending data which are acquired during the day when the tropospheric variability is high (local acquisition time of 15:00 and 3:00 for ascending and descending data, respectively). What is more, the color discontinuity in the ascending data (around center-bottom of Figure 4b, top) is caused by the subtraction of coseismic displacement of aftershock sequence 2.

Furthermore, the viscoelastic relaxation signal during this 1.75 year period is relatively small (up to 4 cm). The small signal is a consequence of the fault geometry. A horizontal fault causes largely horizontal flow.

We also explored afterslip models and concluded that they can be ruled out. The fault plane of the best-fitting model locates a few kilometers above the coseismic rupture in the phanerozoic cover (Figures S2 and S3), which is geologically not plausible. Besides, the best-fitting viscoelastic relaxation model is characterized by a significantly better rms (*rms* = *n i*=1 *residual*<sup>2</sup> *i n* ) than the best-fitting afterslip model (0.015 m versus 0.025 m). We can also rule out afterslip on the coseismic rupture or its downdip extension because of even higher rms. We did not consider combined models with both afterslip and viscoelastic relaxation because the data lack the resolution to constrain all model parameters.

**Figure 5.** (**a**) Regional layered model. The upper crust depth is shown at 4:1 scale for visual clarity; the lower crust and upper mantle depth are shown at 1:1 scale; (**b**) misfit contour map in η*lc* and η*um* space. The yellow star is model 1 (best-fitting model) in Figure 6; the yellow square is model 2 in Figure 6; the yellow triangle is model 3 in Figure 6. The black dash line in (**a**) is the coseismic fault. The black solid line in (**a**) is the afterslip fault. The gray-shaded area in (**b**) denotes 95% confidence interval.

#### *5.5. Summary of Post-Seismic Models*

To visualize the fits of the post-seismic models we combine the ascending and descending displacement fields for the two time periods to the quasi-vertical and quasi-horizontal displacements in east–west direction [40] (Figure 6a,b). The north–south displacement cannot be calculated because

only right-looking InSAR data are available and because polar-orbiting sensors are not sensitive to north–south displacements. We generate simulations for models with lower and higher viscosities in both the lower crust and uppermost mantle. We use η*lc* = η*um* =10<sup>18</sup> Pas for the low viscosity model (model 2) and η*lc* = η*um* =10<sup>20</sup> Pas for the high viscosity model (model 3), respectively.

For time period 1, the low viscosity model produces vertical displacements comparable to the observations, but the east–west displacements are different than observed (model 2, Figure 6a). This reinforces that afterslip was the dominating process during this period. For time period 2, the low and high viscosity models predict significantly more and less displacements than observed, respectively (models 2, 3 in Figure 6f,g), reinforcing that the viscosities are of the order of 10<sup>19</sup> Pas (model 1, Figure 6d).

**Figure 6.** InSAR-derived horizontal east–west and vertical displacements with aftershocks' contribution removed and modelling results. (**a**) and (**b**) InSAR observation; (**c**) afterslip model; (**d**) viscoelastic relaxation model 1; (**e**) and (**f**) viscoelastic relaxation model 2; (**g**) viscoelastic relaxation model 3. Black dot: upper edge of the fault. Red line: nearby faults. Thin dashed rectangle: best-fitting afterslip dislocation and coseismic fault, respectively.
