Supplementary information for A thin crème brûlée rheological structure for the Eastern California Shear Zone

Supplemental details on data analysis, modeling methods, and modeling results, as well as supplemental figures and table.

Several groups have focused on and processed different subsets of the data (Shen et al., 2011). Here, we use GPS time series from the Geodesy Group of the Southern California Earthquake Center (CMM4, Crustal Motion Map Version 4) (Shen et al., 2011) and from the Nevada Geodetic Lab (NGL, geodesy.unr.edu/gps_timeseries/tenv/NA12/, access date: 27 March, 2020). The CMM4 GPS data are mainly derived from 1992-2004 campaign GPS measurements in Southern California ( Figures DR1-DR2). The NGL GPS dataset includes additional data from stations in the Yucca Mountain region and is mainly derived from post-Hector Mine continuous GPS observations (Figure DR2).

s1.1.1 Subsampling NGL GPS time series
We transform the original NGL daily position time series after the first three months into multi-day combination time series by averaging several days' positions. Here, we uniformly set the 1-sigma uncertainties to 2 mm for the horizontal component, which is comparable to those of the CMM4 solution in which time correlated errors have been thoroughly explored. The 1-sigma uncertainties of the vertical component are uniformly set to 4 mm, consistent with the minimum value of the corresponding 1-sigma uncertainties in the CMM4 solution.

s1.1.2 Correcting residual postseismic transients
The difference between the observed position time series and the steady motion at the background rates represents the postseismic relaxation signals. The accuracy of postseismic relaxation time series depends on how well the nominally interseismic GPS velocities represent the background rates. However, the published interseismic velocity fields rely on how well the apparent postseismic transients are removed from post-earthquake GPS observations. As shown in a previous study that considers preseismic triangulation and trilateration data (Liu et al., 2015), residual postseismic transients are not negligible in the nominally interseismic CMM4 and NGL GPS velocities in the southern Mojave. Thus, we have to correct the interseismic GPS velocities, which would otherwise give us biased postseismic relaxation time series. We use the velocities predicted from a block model constrained by pre-Landers triangulation and trilateration data as well as far-field GPS velocities (Liu et al., 2015) instead of the original CMM4 and NGL GPS velocities in the southern Mojave, and retain the CMM4 and NGL GPS velocities outside the southern Mojave.
GPS velocities in Southern California may have been disturbed by other large historical (e.g., the 1857 Mw 7.9 Fort Tejon and 1872 Mw 7.4 Owens Valley earthquakes) and pre-historical earthquakes (Johnson et al., 2007;Pollitz et al., 2008). The pre-Landers rates may include effectively steady motion associated with the long-term viscoelastic relaxation following these events. As both spatial and temporal deformation gradients from these events are likely very modest across the region, we do not attempt to correct for them.

s1.1.3 Combining GPS time series
The CMM4 and NGL GPS time series at co-located GPS sites can be combined, because (1) the reference frame misalignment between CMM4 and NGL solutions is minor, and (2) postseismic relaxation time series from both solutions are consistent during the common observation period from 1999 to 2004. More specifically, post-Hector Mine CMM4 time series can be used to extend the time span of NGL time series at a few sites (e.g., SDHL) which lack data for the early rapid postseismic transients; for the other co-located cGPS sites,

s1.1.4 New GPS time series
The CMM4 and NGL postseismic position time series, after the interseismic velocity correction (Liu et al., 2015), feature enduring late-stage postseismic transients after 2006 in the epicentral area of the two Mojave earthquakes ( Figure DR2). The major patterns are consistent with a solution which is based on a suite of postseismic deformation models (Pollitz, 2015), but ours are better because of the use of pre-Landers data. The GPS data used in most previous studies, however, suggest that the late-stage postseismic transients were possibly complete by 2006 (Herbert et al., 2014).

s2. Supplementary Methods
Coseismic faulting alters the stress field in the crust and upper mantle. The earthquake-induced stress dissipates with time as several physical processes, such as afterslip, poroelastic rebound, and viscoelastic relaxation, occur in response to the stress redistribution, resulting in time-dependent and observable surface deformation. Both the near-field and farfield GPS observations are influenced by viscoelastic relaxation in the upper mantle (Freed et al., 2007;Pollitz, 2015). The near-field GPS signals include additional contributions from shallower deformation processes, including afterslip (Shen et al., 1994;Savage and Svarc, 1997;Pollitz, 2015), poroelastic rebound (Peltzer et al., 1998;Fialko, 2004a), and viscoelastic relaxation in the lower crust (Deng et al., 1998). We therefore have to explore the relative importance of afterslip, poroelastic rebound, and viscoelastic relaxation when interpreting the post-earthquake GPS observations.
We use layered half-space coseismic slip models of the Landers (Fialko, 2004b) and Hector Mine (Simons et al., 2002)  by ~20%, and is therefore not adopted.
To ensure consistency, we use the seismic velocity structure model of the Mojave Desert of Jones and Helmberger (1998)-the same as that used to derive the coseismic slip models (Simons et al., 2002;Fialko, 2004b)-throughout our modeling of kinematic afterslip, poroelastic rebound, and distributed viscoelastic relaxation. To explore how the derived viscosities of the lower crust and upper mantle vary, we also consider other combinations of coseismic slip models and earth structures. For the 2019 Ridgecrest earthquake sequence we use a coseismic slip model constrained by GPS and InSAR (Wang et al., 2020) to drive the viscoelastic relaxation. The seismic velocity model is the same as the one for the Mojave case.

s2.1 Poroelastic rebound
We use the software package PEGRN/PECMP (Wang and Kümpel, 2003) (Peltzer et al., 1998), and the rate of decay is dictated by the hydraulic diffusivity. In the model implementation, we need to consider the trade-off between the contribution from changes in Poisson's ratios and the depth range of fluid flow (controlled by the depth-dependent hydraulic diffusivity). Typically, the permeability (and hence the hydraulic diffusivity) decreases with depth following a power law and perhaps only the very shallow portion hosts active porous flow at the time scale of our observations, with high permeabilities on the order of 10 -13 -10 -14 m 2 (Townend and Zoback, 2000;Jónsson et al., 2003). Another end-member case is that the entire upper crust is abundant in pores and inter-connected fluids and hence hosts poroelastic rebound after the coseismic stress loading (Fialko, 2004a).

s2.2 Viscoelastic relaxation
We use the software package PSGRN/PSCMP (Wang et al., 2006) to model distributed viscoelastic relaxation. Despite geodetic and geophysical evidence of heterogeneous lithosphere beneath Southern California (Pollitz, 2015), for simplicity, we only consider laterally homogeneous models. The layered viscoelastic structures tested in this study have a 16-km-thick elastic layer overlying a 12-km-thick viscoelastic lower crust underlain by a viscoelastic half space ( Figure DR4).
We first consider rheological models with Burgers rheology for both the lower crust and upper mantle, because (1) Burgers and Maxwell rheologies macroscopically characterize creep deformation due to dislocation and diffusion creep, respectively (Bürgmann and Dresen, 2008), (2) dislocation creep may prevail in the entire viscoelastic substrate, (3) Burgers rheology yields non-linear deformation patterns similar to the observed postseismic deformation characterized by rapid exponential decay in the first year and more gradual change afterwards (Pollitz, 2003), and (4) the observed transient deformation in the Yucca Mountain region (250 km away from the epicentral regions) indicates that biviscous (e.g., Burgers) or stress-dependent rheology must be considered for the upper mantle (Freed et al., 2012). However, linear diffusion creep may also occur, and therefore we also consider the case with Maxwell rheology for the lower crust and Burgers rheology for the upper mantle.
For comparison, we also consider the simplest case with Maxwell rheology for both layers.
For the Burgers rheology in our models, the steady-state shear modulus μ MB and bulk modulus k are determined from the previously established layered P-wave velocity, S-wave velocity, and density ( Figure DR4). The yet-to-be-determined parameters are the transient shear modulus μ KB , transient viscosity η KB , and steady-state viscosities η MB of the Burgers body. To reduce the number of unknown parameters (considering the resolving power of GPS data), we fix some of the parameters and only explore the steady-state viscosities of the Burgers body (η MB ) using a grid search. The ratio of the transient viscosity to the steady-state viscosity (η KB /η MB ) is fixed at 0.1, consistent with previous studies (Freed et al., 2012;Pollitz, 2015). We also test other values of η KB /η MB less than 0.1 (e.g., 0.02, 0.05, and 0.075). We find that the estimates of the steady-state viscosities (η MB ) are not sensitive to η KB /η MB if η KB /η MB is less than or equal to 0.1. We find that η KB /η MB should be 0.1 for the upper mantle, otherwise the models fit the far-field data worse. We also find that the preferred η KB /η MB is smaller than 0.1 but greater than 0.02 (e.g., 0.05) when only viscoelastic relaxation is considered for modeling the near-field data. To estimate the effects of transient rheology, we also consider the models constrained by GPS time series in different time periods (the first 1 month, 1 year, and the full 18 years for the Mojave case; the first 1 month and 1 year for the Ridgecrest case.). The ratio of the fully relaxed shear modulus (μ MB μ KB /(μ MB +μ KB )) to unrelaxed shear modulus (μ MB ), parameterized by a in PSGRN, also influences the amount of early rapid transient deformation. We examine a between 0.1 and 0.9 in 0.1 increment (the optimal a is ~0.7).
Because the chosen viscoelastic rheology is linear, the modeled viscoelastic relaxation for both earthquakes can be summed up. However, if power-law rheology applies, the effective viscosity is stress-dependent, and the post-Landers GPS-derived viscosities would differ from those from post-Hector Mine GPS data. Our modeling results show that the rheological structures constrained separately by the post-Landers and post-Hector Mine data are slightly different but consistent, which suggests that the earthquake-induced stress changes are not significantly larger than the background stress. This is an open question that should be further addressed in follow-up studies.

s2.3 Other coseismic slip models and seismic velocity structures
How the estimates of the effective viscosities of the lower crust and upper mantle vary is mainly dictated by the patterns and magnitudes of deviatoric stress change induced by the coseismic rupture, which is determined by the choice of coseismic slip model and seismic velocity model. For that reason, we also consider the coseismic slip models in homogeneous half-space earth model (Wald and Heaton, 1994;Jónsson et al., 2002) and other seismic velocity models such as the one used in Pollitz's previous studies (Pollitz, 2003(Pollitz, , 2015 and the 1D seismic velocity model ( Figure DR4) (Pollitz, 2003(Pollitz, , 2015, except for a difference less than 15% in the upper crust (<16 km depth, Figure

s2.4 Afterslip s2.4.1 Stress-driven afterslip
We consider stress-driven deep afterslip at 16-34 km depth that is simulated using the software package Unicycle (Barbot et al., 2009;Barbot et al., 2017). Rate-strengthening friction applies and the initial slip velocity, V, is determined by in which V 0 is a reference velocity, Δτ is coseismic shear stress change, σ is effective normal stress, and a is the friction parameter. We test V 0 in the range 10 -4 -10 0 m/yr and a in the range 10 -5 -10 0 . The data constraints are the residual time series from the viscoelastic relaxation modeling, with steady-state viscosities allowed to range from 1×10 18 to 1×10 22 Pa s.
Stress-driven afterslip is dictated by the coseismic stress change and the two free model parameters V 0 and aσ, thus the number of free parameters is usually much less than that of kinematic afterslip models (see Table DR1). Therefore, stress-driven afterslip models generally fit the data worse than the kinematic afterslip model. However, kinematic afterslip models are often physically implausible because the distribution of afterslip is not related to the available stress change. In this study, our preferred combined model is the joint inversion of stress-driven afterslip and viscoelastic relaxation.

s2.4.2 Kinematic afterslip
We also estimate the distribution of afterslip on and below the coseismic rupture plane(s) following a kinematic slip inversion approach. The kinematic slip model serves as independent references to our stress-driven afterslip model. In the kinematic afterslip inversion, the observed data d are related to afterslip s by Green's functions G for rectangular dislocations in a layered half-space . To regularize the inversion problem, the second order Laplacian operator L is applied to smooth the afterslip distribution. These can be represented by the following equation where ε denotes data errors that are assumed to follow a Gaussian distribution, e characterizes the uncertainties applied to the smoothing constraints (equivalent to pseudoobservations) that are uniformly set by λ. The trade-off between data misfit and roughness of the afterslip should be explored in this optimization, which is determined by choosing an optimal λ from the data misfit versus model roughness curve.
We subtract viscoelastic response (as described in the Section s3. We impose positivity constraint to get right-lateral slip in afterslip inversion. When the viscosities of the viscoelastic layers are small such that the model predicted fault-parallel motion from viscoelastic relaxation exceed the observed magnitudes, the residual time series (observed minus modeled) require left-lateral slip, which does not make sense. In that case, afterslip is not allowed to occur. For the other cases, however, right-lateral afterslip is required to fit the data. We choose the extent of smoothing constraint applied in the inversion by examining the trade-off curves between slip roughness and data misfit. According to the trade-off curve, a unified value of the smoothing parameter can be used in all the joint inversions. When the smoothing parameter varies around the optimized value, the inferred viscosities almost do not change, but the magnitude of peak afterslip may change a bit (e.g., when λ doubles, the estimated afterslip is rougher, with peak slip increasing by less than 100%, but the estimated viscosities remains at approximately the same values).

s2.5 Fitting criterion
We choose the best-fit models that minimize weighted residual sum of squares (WRSS) of the GPS position time series: where the summation is over Ns stations (j=1,Ns) whose east-west (i=1) and north-south (i=2) components have Nt epochs of data (k=1,Nt, each station has a different Nt), and are the observed and modeled displacements, respectively, ijk is the 1-sigma uncertainty, and is the offset to align with that is estimated from the misalignments ( -) by least-squares. As an alternative to estimating from -, we estimate from , which compulsively align with in the early period; the resultant overall data misfit is worse (increasing WRSS by a factor of ~4) and the estimates of viscosities of the lower crust decrease by a factor of 1.0-1.3.

s3.1 Poroelastic rebound
Are the signals induced by poroelastic rebound discernable by post-Landers campaign GPS and post-Hector Mine cGPS measurements? To address this issue, we first obtain the best-fit viscoelastic relaxation and stress-driven afterslip models and then test whether the addition of poroelastic rebound helps to reduce the data misfit (both horizontal and vertical components). We test two scenarios of poroelastic rebound ( Figure DR5A First, we test the model whether the inclusion of poroelastic rebound reduce the misfit to horizontal component of data. When poroelastic rebound and the viscoelastic relaxation model with Burgers rheology are used to fit the data, we find that the scenario (1) (blue circles in Figure DR5B) fit the data better than the scenario (2) (red circles in Figure DR5B).
We also find that the scenario (1) may help to reduce the data misfit of the models with only viscoelastic relaxation (data misfit is shown by the green line in Figure DR5B), if the hydraulic diffusivity is greater than 0.5 m 2 /s. For the scenario (2), the predicted deformation is distributed more than 20 km from the coseismic rupture ( Figure DR5A), and although the predicted fault-perpendicular motion can be used to remedy model misfit (due to viscoelastic relaxation) to the observed fault-perpendicular motion, the predicted fault-parallel motion will deteriorate the model fit to the data. When the viscoelastic relaxation model with Maxwell rheology is first used to fit the data, the addition of poroelastic rebound does reduce the data misfit (dashed green line versus red and blue squares in Figure DR5B-5C).
The addition of shallow poroelastic rebound also helps to explain the observed vertical deformation in the first two years following the 1992 ( Figure DR6) and 1999 ( Figure   3 in the main text) events. The optimal stress-driven afterslip and viscoelastic relaxation model cannot give a perfect match to the observed vertical deformation within ~10 km from the fault ( Figure 3D and Figure DR6D). The spatial patterns of the mismatch are opposite to model predictions of shallow poroelastic rebound ( Figure 4A and Figure DR6A). When combined with additional contribution from shallow poroelastic rebound, the combined mode fits the post-Landers CMM4 campaign GPS and post-Hector Mine cGPS data better at the 99.9% and 89.4% confidence level. Therefore, consideration of shallow poroelastic rebound, stress-driven afterslip, and viscoelastic relaxation can explain the spatial patterns of the observed vertical deformation ( Figure 4E and Figure DR6E). Note that the major patterns of the modeled shallow poroelastic rebound do not vary ( Figure DR7), even though coseismic slip models with better resolution on the distribution of shallow slip are considered (Xu et al., 2016).
How the poroelastic rebound-induced deformation decays with time is mainly dictated by hydraulic diffusivity. We examine the case with poroelastic rebound and the viscoelastic relaxation model with Burgers rheology. We test several candidate values of hydraulic diffusivity: 0.1, 0.5, 1.0, 1.5, 2.0, and 2.5 m 2 /s. For the cases with poroelastic rebound in the shallow 2.5 km, the goodness of fit is best when hydraulic diffusivity is ~2.5 m 2 /s, but the data misfit is essentially equivalent when the hydraulic diffusivity varies from 1.0 m 2 /s to 2.5 m 2 /s (blue squares in Figure DR5C).
In summary, GPS measurements, although sparsely distributed, support shallow poroelastic rebound, rather than the upper crustal scale poroelastic rebound (Fialko, 2004a), for the Mojave case. InSAR data with better spatial coverage along the 1992 and 1999 ruptures (in particular the stepovers) is required for a better understanding of how poroelastic rebound developed after the two earthquakes. The following models with viscoelastic relaxation include a contribution from poroelastic rebound in the upper 2.5 km with a hydraulic diffusivity of 1.0 m 2 /s. However, the viscosity estimates are not sensitive to whether the shallow poroelastic rebound is considered.

s3.2 Univiscous lower crust and biviscous upper mantle
We consider the cases with Maxwell rheology for the lower crust and Burgers rheology for the upper mantle (the ratio of the transient viscosity to the steady-state viscosity is fixed at 0.1), because (1) the microscale physical processes responsible for ductile deformation may be diffusion creep and dislocation creep for the lower crust and upper mantle, respectively, and (2) the observed transient deformation in the Yucca Mountain region (250 km away from the epicentral regions) support biviscous (e.g., Burgers) or stressdependent rheology for the upper mantle (Freed et al., 2012).
The steady-state viscosity of the lower crust can be ~1.0×10 20 Pa s, less than three times that of the upper mantle, in a joint inversion of both kinematic/stress-driven afterslip and viscoelastic relaxation ( Figure DR8 and Table DR1). The estimated stress-driven afterslip lasts for more than 4 years and does not exhibit rapid early transients within the first year ( Figure DR9). The best-fit models fit all the 18-years GPS data worse than the corresponding best-fit models with Burgers rheology for the lower crust (see Section s3.3) at the 99% confidence level.
We also test the effects of the ratio of the transient viscosity to the steady-state viscosity (Burgers rheology for the upper mantle). When that ratio is less than 0.1 (e.g., 0.02, 0.05, and 0.075), the best-fit models with only viscoelastic relaxation fit the data (both the near-field and far-field data, in particular, the far-field data) worse than the best-fit model with that ratio being 0.1. The data do not allow us to make a distinction between the models with that ratio being 0.1 and 0.075, because data misfits are statistically comparable at the 55.55% confidence level. However, the models with that ratio ≤0.05 fit the data worse than the model with that ratio being 0.1 at >95.00% confidence level. Note that the estimates of the steady-state viscosity are not sensitive to the tested ratios of the transient viscosity to the steady-state viscosity (similar conclusions apply to the Section S3.3).

s3.3 Biviscous lower crust and upper mantle
We consider the cases with Burgers rheology (the ratio of the transient viscosity to the steady-state viscosity is fixed at 0.1) for both the lower crust and upper mantle, because dislocation creep might occur in both layers. The rheological parameters are presented in the Table DR1. Shown in Figure DR10 is misfit to GPS data in the time period 1992-1999, 1999-2010, and 1992-2010. Shown in Figure DR11 are the distribution of afterslip of the 1992 and 1999 events. Shown in Figure DR6  When additional contribution from stress-driven afterslip is considered ( Figure DR10-DR11), the best-fitting steady-state viscosity of the lower crust is ~2×10 20 Pa s (Table DR1), ~5 times of that of the upper mantle. The stress-driven afterslip induces early rapid deformation and ends within the first 3 years (Figures DR14-DR15). Viscoelastic relaxation in the lower crust and upper mantle induce the slowly decaying but long-lasting displacements in the near field.
When additional contribution from kinematic afterslip is considered ( Figure DR10-DR11), the steady-state viscosity of the lower crust can be ~9×10 20 Pa s (Table DR1), an order of magnitude larger than that of the upper mantle. That does not indicate that the relaxation of the lower crust has not come into play in the first two decades. We believe that the role of the derived kinematic afterslip is exaggerated and the effects of viscoelastic relaxation in the lower crust are underestimated. That is because: (1) The derived kinematic afterslip lasts for at least one decade (in particular for the afterslip associated with the 1992 event, see Figure DR11) and induces steady-deformation that should be attributed to bulk viscoelastic relaxation.
(2) The distribution of afterslip is not dictated by the coseismic stress loading (especially for the kinematic afterslip of the 1999 Hector Mine event, see Figure   DR11G).
We also test the effects of the ratio of the transient viscosity to the steady-state viscosity. When that ratio is less than 0.1 (e.g., 0.05 and 0.075) but greater than 0.02, overall, the best-fit models with only viscoelastic relaxation fit the data better than the best-fit model with that ratio being 0.1, at >94.25% confidence level (The data misfit decreases as that ratio decreases until when the ratio is 0.02. When the ratio is 0.05, the model fit the data better at the 99.47% confidence level.). The relevant models with smaller values of that ratio fit the near-field data better, but fit the far-field data worse. However, F-test shows that the far-field data misfits are statistically comparable at the <71% confidence level, when that ratios are in the range from 0.05 to 0.1. The transient viscosity of the lower crust constrained by all the 1992-2010 data is estimated to be ~1×10 19 Pa s, less than 4 times that of the upper mantle.
The transient viscosity of the lower crust is estimated to be ~1×10 18 Pa s and ~3×10 18 Pa s from the GPS data in the first month and in the first year following the two earthquakes, respectively. These modeling tests, together with other tests in this section, reveal that transient deformation in the lower crust is not negligible after the two Mojave earthquakes.

s3.4 Transient rheology inferred from post-Ridgecrest GPS measurements
To years (Pollitz et al., 2001;Pollitz, 2003). Because of insufficient far-field data constraints, the transient viscosity of the upper mantle is not well constrained, but is estimated to be on the order of 10 17 Pa s. The data used in our test include more far-field GPS data, which provide robust constraint on the upper mantle rheology. We therefore believe that our estimates of the short-term upper mantle rheology are reliable.
To sum up, based on the data in the first 1 year, the transient rheology of the upper mantle beneath the Ridgecrest region is consistent with that of the upper mantle beneath the southern Mojave. However, because that (1) 1992-1999CMM4, 1999-2010NGL, 1999-2010CMM4+NGL, and 1992-2010 CMM4+NGL GPS time series.

s4. Supplementary Table
b Green and yellow shaded rows indicate the models with the same bulk rheology but constrained by different datasets. Maxwell and Burgers rheologies have 1 and 3 unknown parameters, respectively. Stress-driven afterslip inversion has 2 unknowns. Kinematic afterslip inversion for each event has ~15 well-resolved parameters. For the Burgers rheology, the ratio of the transient viscosity to the steady-state viscosity is 0.1.
c Viscoelastic relaxation and shallow poroelastic rebound. The same applies to the first entry of the three entries for each modeling test. So do the second and third entries in d and e.
d Kinematic afterslip and viscoelastic relaxation (+shallow poroelastic rebound). e Stress-driven afterslip and viscoelastic relaxation (+shallow poroelastic rebound).    (Xu et al., 2016) and a modified slip model in the same layered Earth model (Fialko, 2004b), respectively. (B) Blue and green vectors are the modeled displacements in a layered Earth model, based on coseismic slip model in homogeneous half-space (Xu et al., 2016) and a modified slip model in the same layered Earth model (Simons et al., 2002), respectively.     and lower crust (η LC ), respectively. The colored lines show the statistical 95% confidence regions of the models with respect to the optimal models marked by the color-coded diamonds with only viscoelastic relaxation (Model 1, red contours), with additional contribution from stress-driven afterslip (Model 2, blue contours), and with additional contribution from kinematic afterslip (Model 3, green contours). The yellow diamond and contour are for Model 2 only constrained by vertical GPS data. The black diamond marks the optimal viscosities without afterslip and without considering the interseismic velocity correction (Liu et al., 2015).  Horizontal and vertical axes show the steady-state viscosity of the upper mantle (η UM ) and lower crust (η LC ), respectively. The colored lines show the statistical 95% confidence regions of the models with respect to the optimal models marked by the color-coded diamonds with only viscoelastic relaxation (Model 1, red contours), with additional contribution from stress-driven afterslip (Model 2, blue contours), and with additional contribution from kinematic afterslip (Model 3, green contours). The yellow diamond and contour are for Model 2 only constrained by vertical GPS data. The black diamond marks the optimal viscosities without afterslip and without considering the interseismic velocity correction (Liu et al., 2015).