Earth and Planetary Science Letters Crustal rheology of southern Tibet constrained from lake-induced viscoelastic deformation

We probe the rheology of the Tibetan lithosphere using the rebound that accompanied climate-driven lake level variations. At the modern decadal time scale, we used deformation around Siling Tso measured from InSAR. At the millennial time scale, we use Holocene paleoshorelines around Siling Tso and Zhari Nam Tso. We use chronological constraints from the literature and Digital Elevation Models to constrain their ages and geometry. We observe a small post-highstand subsidence of the area near the center of mass of the paleolake-load and a low-amplitude short-wavelength outer bulge. In the context of a model consisting of an elastic lid over a viscous channel with a rigid base, these observations preclude the existence of a thick low viscosity channel and require a thin elastic lid. Based on Monte Carlo inversion, we constrain the range of possible equivalent elastic thickness of the lid ( < 5 km), the viscosity (2 × 10 18 –10 20 Pa.s) and thickness of the crustal channel ( < 10–20 km). By contrast, the modern data requires a stiffer lid with equivalent elastic thickness > 20 km and a > 20 km thick channel with lower crustal viscosity ( < 5 × 10 18 Pa.s). The different rheologies inferred at these different time-scales could be explained by a Burgers body rheology of the middle and lower crust, with a short-term viscosity of 10 18 Pa.s and long-term viscosity of 10 20 Pa.s, or even better by vertical variations of viscosity. To illustrate the latter claim, we show that the observations at the decadal and Holocene time scales can be reconciled by assuming a low viscosity zone (10 18 Pa.s) at mid-crustal depth (between ∼ 10 and 30 km depth) embedded in a higher viscosity crust ( > 10 20 Pa.s). In both cases, the interferences in space of the deformation signals induced by the lakes geometry, and in time through the viscoelastic response to the lake level variations results in limited distortion of the paleo-shorelines. While the elastic lid in the upper crust needs in any case to be thin ( < 10 km), the low amplitude distortion requires signiﬁcant viscoelastic support from the lower crust and upper mantle; this explains the relatively high effective elastic thickness ( > 20 km) inferred in some previous studies of Holocene paleoshorelines. In the longer term, the effective elastic thickness of the lithosphere must drop asymptotically to the value of the elastic lid in the upper crust ( < 10 km); this explains the low effective elastic thickness derived from gravity studies. © 2018 The by Elsevier This is CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The rheological stratification of the continental crust remains a subject of debate (e.g., Bürgmann and Dresen, 2008;Burov et al., 2014;Jackson, 2002). This debate is central to our understanding of continental tectonics, the formation and evolution of mountain ranges and orogenic plateaus in particular (e.g., Beaumont et al., 2001;Copley et al., 2011) and of the seismic cy-argued for "channel flow" invoking a very weak and laterally mobile lower crust. For instance, a large horizontal flux of middle and lower crust squeezed out from beneath the plateau could explain some aspect of the morphology and tectonics of eastern Tibet and the Himalaya (Beaumont et al., 2001;Clark and Royden, 2000). Such channel flow tectonics would imply a viscosity possibly as low as 10 17 Pa.s (e.g., Clark et al., 2005;Royden et al., 1997). However, several geological and geophysical data seem inconsistent with such a weak lower crust. Seismic anisotropy observations suggest coherent deformation throughout the eastern Tibetan lithosphere and hence little decoupling at mid-crustal depth (León Soto et al., 2012). The contrast between the northern part of the plateau, which is dominated by strike-slip faulting, and the southern part, dominated by normal faulting is also an evidence for a strong coupling between the lower crust and upper mantle in southern Tibet (Copley et al., 2011).
Direct constraints on the rheology of the lithosphere may be derived from observing a time-dependent deformation response to a known stress perturbation. At the decadal timescale the viscosity can be estimated from postseismic observations following large earthquakes. Such studies have yielded values between 10 17 Pa.s to more than 10 21 Pa.s (Hilley et al., 2005;Huang et al., 2014;Ryder et al., 2011Ryder et al., , 2014Yamasaki and Houseman, 2012;Zhang et al., 2009). A significant source of uncertainties in such studies is due to the fact that it is difficult to separate the contributions of afterslip and viscous relaxation to postseismic deformation. In addition, we do not know whether viscous relaxation beneath the seismogenic zone is broadly distributed or localized. The rheology inferred from postseismic studies might therefore not be representative of the bulk rheology of the crust.
Crustal rheology may alternatively be probed from the surface deformation induced by lake level variations (Bills et al., 1994;Kaufman and Amelung, 2000). For instance Doin et al. (2015) studied the Siling Tso example (Fig. 1a). The water level has been rising up recently by about 1 meter per year, flexing down the topography around it. According to Doin et al. (2015), the model that best fits the deformation signal measured from InSAR has an elastic upper crust with an equivalent elastic thickness (Te) of ∼30 km and a lower crust with a viscosity of 1-3 × 10 18 Pa.s. It is also possible to estimate the rheology of the crust at the millennial time scale based on surface deformation resulting from lake level fluctuations induced by late Quaternary climate change (e.g., Bills et al., 1994Bills et al., , 2007. This method has also recently been used in Tibet where numerous closed basins bear well-preserved paleoshorelines (England et al., 2013;England and Walker, 2016;Shi et al., 2015). Surprisingly, the paleoshorelines are hardly distorted implying either a large equivalent elastic thickness (≥25 km) or a high viscosity crust (≥10 19 Pa.s).
In this study, we re-analyze and seek to reconcile the rheology inferred from the lake-induced deformation signal at the decadal and millennial time scales. The motivations for the re-analysis are multiple. Firstly, the trade-off between the equivalent elastic thickness of the crust and the thickness and viscosity of the viscous channel was not fully explored in previous studies. Secondly, previous studies have assumed very simple paleo-lake level variations which can now be refined based on recent studies of Tibet paleolakes (Ahlborn et al., 2016;Chen et al., 2013;Hudson et al., 2015;Lee et al., 2009;Rades et al., 2015;Shi et al., 2017). Thirdly, the possible effect of loading and unloading by mountain glaciers was ignored in previous studies; there is however clear geomorphic evidence that the mountain ranges around these lakes underwent glacial advances. Finally, there are now good quality images and better DEMs which can be used to measure the distortion of the paleoshorelines more extensively.

Previous results inferred from lake-induced deformation in Tibet
At the modern decadal time scale (Doin et al., 2015) studied the ground deformation due to Siling Tso level rise of 1.0 m/yr from 2000 to 2006. They measured the ground deformation using interferometric synthetic aperture radar (InSAR) for the 1992-2011 period. Their models considered a crust consisting of an upper elastic lid over a viscoelastic channel adding to an imposed thickness of 65 km. Here, we revisit their analysis by allowing the thickness of the viscoelastic channel to be independent of the elastic lid thickness. The rationale is that the lower crust could be granulitic and rather cold due to the effect of underthrusting of India on the thermal structure and might therefore not be part of the low viscosity channel.
Late Pleistocene to Holocene paleoshorelines are widespread all over Tibet (e.g., Gasse et al., 1991;Ahlborn et al., 2016;Avouac et al., 1996;Chen et al., 2013;England et al., 2013;Hudson et al., 2015;Lee et al., 2009;Rades et al., 2015;Shi et al., 2015). Five major closed-basin lake systems in the central part of the plateau have already been studied to constrain rheological properties of the Tibetan crust (England et al., 2013;England and Walker, 2016;Shi et al., 2015). The Ngangla Ring Tso, Taro Tso, Zhari Nam Tso and Tangra Yum Tso were studied by England et al. (2013) whereas Shi et al. (2015) analyzed the shorelines around Siling Tso (Fig. 1a). The paleoshorelines were measured using Shuttle Radar Topographic Mission elevations combined with Google Earth imagery of shorelines and some kinematic GPS measurements around the Zhari Nam Tso (England et al., 2013). These lakes are among the largest in Tibet and filled up to 150-200 m above their present level in the Holocene. A large deformation response to unloading might therefore have been expected, but no conspicuous deformation signal was found at any of them. The shorelines are distorted from horizontality by no more than ∼10 m. These observations place important constraints on the rheology of the Tibet crust, but the previous studies leave room for some reanalysis. These studies neglected the potential effects of surrounding glacial unloading, the simulated lakes level histories were very simplified and the small datasets were too sparse to detect short-wavelength distortion. We therefore revisit this problem with more complete load scenarios and an augmented dataset of shoreline elevations that we produced using satellite images.

Geomorphology of Holocene paleoshorelines
We focused on the best-preserved, presumably most recent shorelines, which also have the best-constrained loading history. Those formed on gentle alluvial slope, such as the ridges seen in Fig. 1b, are not very durable and probably of Holocene age in the context of Tibet. The highest stand of the Zhari Nam Tso preserved in the morphology has been found to postdate the Last Glacial Maximum (e.g., Kong et al., 2011). Field observations show that they commonly form couplets, separated in height by ∼0.5-1 m, which could reflect a seasonal effect (England et al., 2013). The elevation of a particular shoreline can vary laterally but generally within less than 0.5 m (Shi et al., 2015). Therefore, the intrinsic variability of shoreline elevation is estimated to be of the order of ∼1 m.
We used ∼2.5 m ground resolution satellite images to pick the highest well-preserved shorelines visible on alluvial surfaces around Zhari Nam Tso and Siling Tso ( Fig. 1b and Fig. 2). Their elevations is estimated by sampling the global DEM ALOSWorld 3D-30 m (AW3D30) released in 2015 by the Japan Aerospace Exploration Agency (JAXA) using a bicubic interpolation. We selected showing well preserved shorelines marking the Holocene highstand of Zhari Nam Tso (white arrows). The sequence of paleoshorelines at lower elevation recorded the regression of the lake from the mid-Holocene to its present level. (c) Estimated lake level variation of Ngangla Ring Tso (Hudson et al., 2015), Taro Tso (Lee et al., 2009), Zhari Nam Tso , and Tangra Yum Tso (Ahlborn et al., 2016;Rades et al., 2015) for the last 25 ky. Each curve is normalized to the highest stand. only paleoshorelines located on gentle slopes, based on our visual assessment, to minimize elevation errors due to the possible misregistration of the images to the DEM (Fig. 1b) Fig. 3a and b show the distribution of the difference between the raw and the filtered elevations data which were smoothed using a 10 km square sliding window. These approximately normal distributions are used to characterize the noise on shoreline elevations. We estimate the 67% uncertainties on shoreline elevation to ∼1.5 m at both Siling Tso and Zhari Nam Tso. The distribution of elevations smoothed at the 10 km scale reveal a variability (∼6 m for each datasets) that exceeds this uncertainty. Given that we do not expect the noise on shorelines elevation to correlate at length scales larger than 10 km, we interpret this pattern as a real deformation signal. The significance of this signal is not immediately intuitive however. It does not show the typical "bowl shape" uplift pattern centered on the centroid of the water load which would be expected from a simple elastic rebound. The central part of the Observed (color-coded circles) and modeled elevation at present of the paleoshorelines marking the Holocene highstand around Siling Tso (a) and Zhari Nam Tso (b). The color scale is centered on their mean restored elevation. The models assume an elastic lid overlying a viscoelastic channel with a rigid base. We selected models representative of the best fitting solutions. For Siling Tso, the elastic lid thickness is Te = 2 km, the viscoelastic channel is L = 9 km thick and the viscosity is η = 1.8 × 10 18 Pa.s. For Zhari Nam Tso the elastic lid thickness is Te = 6.4 km, the viscoelastic channel is L = 5.9 km thick and the viscosity η = 5 × 10 19 Pa.s. The model prediction for Zhari Nam Tso, takes into account the deformation associated to the surrounding lakes (Ngangla Ring Tso, Taro Tso and Tangra Yum Tso). The model prediction for the whole region including the 4 lakes is represented in Fig. 4d. The data indicate deformation of small amplitude at a wavelength of about 100 km. Note also the counterintuitive observation of subsidence nearer the paleolake centers, which is most obvious for Siling Tso. The models are able to produce small amplitude deformation at a wavelength consistent with the observations as well as subsidence nearer the paleolake centers. Comparison between predicted and observed elevations are shown in Supplementary  Fig. E1.
Siling Tso basin was apparently deflected downward, whereas uplift would have been expected in that area of maximum unloading. This central zone of subsidence is fringed by an outer bulge of uplift which is a few tens of kilometers wide. The Zhari Nam Tso shorelines show a similar pattern, though more complicated probably due the more contorted geometry of the paleolake and the interference with the nearby lakes. We characterize the wavelength of the deformation signal using the standard deviation of the paleoshoreline elevations within a sliding square window of varying size between 0 to 100 km ( Fig. 3c and d). We use the filtered data (obtained by averaging within a 10 km × 10 km wide sliding window) as this filtering helps comparison with the model predictions which have no noise. Windows with less than 20 data (chosen arbitrarily) were discarded to filter out poorly constrained values. The standard deviation of both datasets increases rapidly with the window size up to ∼70-80 km for Siling Tso and Tso ∼60 km for Zhari Nam Tso and levels off for larger windows. It suggests that the deformation signal has a wavelength of the order of ∼60-80 km, smaller than the wavelength expected from the rhe-ological models derived from previous studies (Shi et al., 2015;Doin et al., 2015) (Fig. 3c and d).

Lakes level variation over the Holocene
A number of paleoclimatic studies have documented the timing of the highstand and of the regression of the lakes analyzed in this study (Ahlborn et al., 2016;Chen et al., 2013;Hudson et al., 2015;Lee et al., 2009;Rades et al., 2015;Shi et al., 2017). Siling Tso ( Fig. 1d) reached its highstand probably at ∼10 ka and started to regress progressively to its present-day level at ∼4 ka, possibly due to weakening of the monsoon (Shi et al., 2017). The three lakes around Zhari Nam Tso (Ngangla Ring Tso, Taro Tso and Tangra Yum Tso) followed a slightly different history (Fig. 1c). They also reached their highstand after the Last Glacial Maximum and maintained a highstand during the early Holocene climatic optimum ( Fig. 1c) when rainfall was more abundant than at present but started to regress earlier than Siling Tso at ∼8.5 ka (Hudson et al., 2015). This time evolution is somewhat different from the abrupt regression at ∼5 ka assumed by England et al. (2013). The The histograms show the distribution of elevations of the raw data (green) and the filtered data (blue) obtained by averaging within a 10 km × 10 km sliding-window. The histogram of the differences between the filtered and unfiltered data (red) is the one considered to characterize the uncertainties on the paleoshoreline elevation for Siling Tso (standard deviation ∼1.5 m) and Zhari Nam Tso (standard deviation ∼1.4 m). The lower plots show the normalized standard deviation calculated within a sliding square window of 0 to 100 km width, of the post-highstand elevation change derived from the paleoshorelines elevation (black symbols). The data were filtered by averaging within a 10 km wide sliding-window to remove high frequency noise. This filtering helps comparison with the model predictions which have no noise. Error bars show the uncertainty on the calculated standard deviation at the 68% confidence level (1-σ ). Windows with less than 20 data (chosen arbitrary to avoid poorly constrained values) were discarded. The standard deviation increases rapidly with the window size and levels off for a window size exceeding ∼60 km. This pattern indicates that the deformation signal has a wavelength of the order of 60 km. Prediction from one of the best models is shown for both Siling Tso and Zhari Nam Tso (Model 1, same as Fig. 2). Model 2 corresponds to an elastic model with an elastic thickness Te = 10 km over an inviscid medium at the lower end of the elastic models of England et al. (2013). Model 3 corresponds to a model with the best set of parameters determined by Doin et al. (2015); an elastic lid of thickness Te = 30 km overlying a viscoelastic channel of viscosity η = 2 × 10 18 Pa.s and thickness L = 35 km over a rigid base. The wavelength associated with models 2 and 3 is much larger than 60 km, as a result the standard deviation increases nearly linearly with the window size over the range of tested values contrary to what the data show.
lack of data for the early Holocene and Late Pleistocene makes it difficult to estimate when the lakes reached their highstand. The available data show disparities that might reflect a variable contribution from ice melting depending on the particular setting of each basin. For simplicity, we assume that all four lakes around Zhari Nam Tso followed the same time evolution represented by a post LGM highstand plateau of variable duration, t h , and a gradual regression starting at 8.5 ka (Fig. 1c).
To estimate the load induced by the five paleolakes, we extracted the contour lines corresponding to their mean highstand elevation (Ngangla Ring Tso: 4864 m, Taro Tso: 4606 m, Zhari Nam Tso: 4751 m, Tangra Yum Tso: 4741 m and Siling Tso: 4597 m). The difference between the highstand and the local ground elevation within each basin provide an estimate of the spatial distribution of the drop of surface load that resulted from the lake regression.
This estimate can be corrected from the post-regression rebound with a few model iterations.

Variation of glacial extent
Well-preserved moraines are observed at the outlet of most valleys carved into the ranges surrounding the lakes. The major glacial advance preserved in Central Tibet probably occurred in the early Holocene (Owen and Dortch, 2014). If so, these mountain glaciers might have contributed to the deformation of the Holocene paleoshorelines. We therefore estimated the surface load variations that would have resulted for the glacial retreat following the maximum extent preserved in the morphology. To do so, we used the Gc2d thin-sheet ice model of Kessler et al. (2006) (Appendix B). This 2-dimensional numerical model simulates the growth of glaciers for a given topography and meteorology. Using an explicit finite difference scheme by solving ice flux and mass conservation equations the model returns the ice elevation distribution through time. We estimated a lower and upper bound of glacial extent, and a medium case ( Supplementary Fig. B.1).

Forward modeling
To model the distortion of the shorelines different simplified representations of the rheology of the crust were considered. A number of analytical solutions allow calculating the elastic or visco-elastic deformation of a layered planar earth model submitted to a time-dependent vertical load at the surface. We first consider an elastic thin plate over an inviscid fluid to estimate the asymptotic deformation response expected after complete viscous relaxation. We use the analytical solution of Brotchie and Silvester (1969) (Equation C.2, Appendix C). The viscoelastic response is calculated using the approach of Nakiboglu and Lambeck (1982) (Equation C.3 and C.4, Appendix C). It assumes an elastic lid over a channel with a Maxwell rheology and a rigid base. The impact and significance of rigid base assumption is discussed below. These analytical solutions can be used to describe the response to any instantaneous variation of surface load. The load history is simulated by discretizing the time variations of surface load as a series of step functions.
We run the models forward based on the assumed lake level history and then compare the observed and predicted distortions of the paleoshorelines elevation. Given a surface load history, which is estimated based on the present topography and the assumed lake level variations, the models allow predicting vertical displacements that would have affected an initially undeformed horizontal shoreline at the onset of the lake regression.
As a first step, we adopt the same approach as England et al. (2013) for the purely elastic case. We vary the elastic thickness Te from 1 to 40 km and compare the measured difference of elevation between the maximum and the minimum elevation of the paleoshorelines (∼6 m) with the range of vertical displacements. This approach makes sense if the spatial variations of paleoshorelines elevations reflect measurement errors rather than a true pattern of vertical displacement. If the spatial pattern is a real deformation signal, it is then more appropriate to compare the observed and the modeled distortions at the location of the data. The models can be tested by retrodeforming the paleoshorelines, i.e., by subtracting the predicted uplift from the present elevation at each data point. The best models are those which best restore the paleoshorelines to horizontality. A natural goodness of fit criterion is therefore the standard deviation of the retrodeformed elevations.
A model output is the estimated mean elevation of the paleoshoreline at the time of deposition ( Z 0 ). This information can then be used to correct the initial estimated load. We initially started with estimating the load based on the elevation of the paleoshorelines above the present lake level. This estimate ignores the distortion of the paleoshorelines. The load can be re-estimated based on the mean elevation of the paleoshoreline above the retrodeformed topography. The model can thus be adjusted iteratively until it is self-consistent. Fewer than 5 iterations are sufficient to obtain a convergence within 0.1% with the purely elastic models. As a result of these iterations 10 to 20% of additional uplift  (α1+α2) is predicted. However these iterations are not used in the inversions procedure described below due to their computational cost and their second order effect on the surface response to load regression.

Goodness of fit criterion
We define here the goodness of fit criterion used to quantify the discrepancy between the model predictions and the observations. The analytical calculation returns vertical displacement values M relative to the initial horizontal highstand at elevation Z 0 which is a priori unknown. The predicted distortion M(x, y) should equal Z obs (x, y) − Z 0 . It follows that, for each location Z obs (x, y) − M(x, y) is the estimated initial highstand elevation Z 0 .
The best model is the one that restores best the paleoshoreline to horizontality, so the one which minimizes Z obs (x, y) − M(x, y) − Z obs − M , where is the arithmetic average. Thus the best fitting model corresponds to the one that minimizes the dimensionless reduced Chi-square where n is the number of observations, p the number of varying parameters and σ the standard deviation.
The best estimate of the highstand mean elevation Z 0 is Z obs − M that is thus a model output. A model that fits the data within uncertainties yields a χ 2 r value of the order of unity. Chi-squared statistics can then be used to estimate the uncertainties on the model parameters.

Inversion
The adjustable parameters in our inversions are, the thickness Te of the elastic lid, the thickness L and viscosity η of the underlying viscous channel and the highstand duration t h (for Zhari Nam Tso only). We explore the space of model parameters using a Monte Carlo method with the built-in matlab slicesample function (Neal, 2003). This procedure results in a higher density of samples in regions of lower misfit.

Elastic flexure
We first consider the case of an elastic layer of thickness Te over an inviscid fluid (Fig. 4a). The surface load can be related to the lakes alone or the lakes and the surrounding glaciers. This model gives an estimate of the elastic thickness Te which is necessary to support the load with elastic stresses only. As viscous support is ignored, Te estimated in this way should be considered as an upper bound. For simplicity we assume that loading by glaciers and lakes is synchronous. This is not realistic as these two types of loads must have had a different time history. The calculation does, however, quantify the possible effect of the glaciers on the shoreline distortion. Fig. 4a compares the maximum distortion derived from the mapped paleoshorelines around Zhari Nam Tso, with the distortion predicted in response to the lakes regression with, or without the effect of the glaciers retreat. Fig. 4b shows the spatial distribution of vertical rebound for Te = 10 km, at the lower end of the values proposed by England et al. (2013). As expected, we see local maxima of the rebound near the centers of the paleolakes. Neighboring lakes do not interfere much in that case. They start to influence each other for larger values of Te. The flexural parameter is about 23 km in that case and it increases as Te 0.75 (see relationships between Te and β in Table 1). The observation of no more than ∼6 m of distortion requires Te ≥ 35 km, a lower bound consistent with the results of England et al. (2013) and Shi et al. (2015). The lower bound is even larger value if the effect of glaciers retreat is included (Fig. 4a). Glacial unloading has little effect for Te < 10 km because most of the surrounding glaciers are >185 km away from the paleoshorelines. By contrast, if Te >10 km, the glaciers could potentially have an effect on the shoreline distortion. We conclude that the influence of the glaciers is probably very small and does not help explain the short wavelength and low amplitude distortion of the paleoshorelines.

Viscoelastic flexure
A load history is required in viscoelastic models. For Zhari Nam Tso and the 3 lakes around we assume a linear regression from 8.5 ka to present (Fig. 1c). We use the Zhari Nam Tso example to assess the effect of the duration of the highstand, t h , which could be of the same order of magnitude as the relaxation time. The time when the lakes reached their highstand is varied between Fig. 5. Results of the inversion of the Zhari Nam Tso highstand paleoshoreline assuming an elastic lid over a viscoelastic channel with a rigid base. The load due to the surrounding lakes (Ngangla Ring Tso, Taro Tso and Tangra Yum Tso) is taken into account. The 2-dimensional slices into the 4 parameters inversion show the reduced chi-square, as defined in the text (Equation (1)), as a function of the thickness and viscosity of the viscoelastic channel. The range of tested elastic thickness increases from left to right (0-10 km, 10-20 km and 20-30 km). From top to bottom the range of highstand duration increases: (a) 1.5 to 4.5 ky (highstand reached between 10-13 ka), (b) 4.5-8.5 ky (highstand reached between 13-17 ka) and (c) 8.5-11.5 ky (highstand reached between 17-20 ka). 10 ka and 20 ka, which gives a highstand duration t h of 1.5 ky to 11.5 ky. In absence of reliable time constraints on the glacial history, we assume that the glaciers grew from 26 ka to 24 ka and that they retreated as the lakes were gaining volume, presumably supplied by glaciers melting. Our results, detailed below, show that the duration of the highstand does not trade off with any other model parameters. So we did not include the duration of the highstand as a parameter in the case of Siling Tso. For the Holocene history of Siling Tso (Fig. 1d), we assume a linear rise of the lake level from 10.5 ka to 9.5 ka and a regression from 4 ka to present. The duration of the highstand is fixed to 5.5 ky. Finally, the Siling Tso is assumed to have risen by 10 m approximately from 2000 to 2006 based on Doin et al. (2015).
We first compare the paleoshorelines observations with the predictions from the preferred model of Doin et al. (2015) derived from the deformation response to the modern rise of Siling Tso (Fig. 4c). We did not use the original raw InSAR data as it would have entailed redoing the analysis in its entirety including the implementation of a non-standard procedure to separate the signal from the noise. So our analysis assumes that the bestfitting model of Doin et al. (2015) (Supplementary Fig. F.1) is a valid model within the parameter space that was explored in this study. This model has an elastic lid of thickness Te = 30 km overlying a viscoelastic channel of thickness L = 35 km and of viscosity η = 2.10 18 Pa.s over a rigid base. At Zhari Nam Tso it predicts a distortion of less than 6 m, a value consistent with the paleoshorelines measurements. The reduced Chi-square corresponding to this model, without the glaciers, is ∼2.3. Like the purely elastic model, this viscoelastic model fails to fit the observed short-wavelength distortion pattern. The wavelength associated with this model is much larger than ∼60-80 km, as a result the standard deviation of the distortions increases nearly linearly with the window size  The Monte Carlo search for the best set of parameters starts from an initially random set of parameters within a defined range: η = 10 17 -10 21 Pa.s, L = 1-50 km, Te = 1-30 km and for Zhari Nam Tso t h = 1.5-11.5 ky. Fig. 5 displays the misfit (reduced Chi-square, χ 2 r , Equation (1)) for all the realizations associated to the inversion of Zhari Nam Tso paleoshorelines as a function of η and L, with panels corresponding to binned values of Te or t h . The better models, with reduced χ 2 r close to 1.82 fall in a domain defined by an elastic thickness Te <10 km, a viscosity η of 10 18 Pa.s to 10 20 Pa.s and a channel thickness L < 10-20 km. The plot shows a strong trade-off between the channel thickness L and viscosity η: a thicker channel requires a higher viscosity to maintain a comparable fit to the data. Such a trade-off is expected as the relaxation time for channel flow is τ ∼ 2 ηR 2 ρ gL 3 (with R being the radius of the load or the elastic flexural parameter whichever is larger, L the channel thickness, η the viscosity of the viscous medium and ρ its density). The highstand duration t h does not trade-off with any other model parameters in the range of tested values (1.5 ky to 11.5 ky). Fig. 6 shows the result of the inversion of the paleoshorelines around Siling Tso (labeled "mid-term", Fig. 6a) and of the modern ground deformation (labeled "short-term", Fig. 6b). In the later case we use as an input in the inversion the vertical displacements predicted by the best model of Doin et al. (2015) which fits all the InSAR data used in the analysis within their uncertainties ( Supplementary Fig. F.1, data provided as electronic sup-plement). Thus we try to determine models parameterized with 3 variables (elastic thickness, channel thickness and channel viscosity) equivalent to the preferred model of Doin et al. (2015) which has only two independent parameters has an elastic lid of thickness Te = 30 km overlying a viscoelastic channel of viscosity η = 2 × 10 18 Pa.s and thickness L = 65 − 30 = 35 km over a rigid base. At the millennial time scale, the results are similar to those from Zhari Nam Tso (Fig. 5). The plots in Fig. 6a show a strong trade-off between the channel thickness and viscosity: a thicker channel (from 10 km to 20 km) requires a higher viscosity (from 10 18 Pa.s to 10 19 Pa.s) to maintain a comparable fit to the data (χ 2 r ∼6). However a small trade-off between the channel viscosity and the elastic thickness is found at the Holocene time scale. The inversion results at the decadal time scale are similar to those of Doin et al. (2015), showing that the better models (normalized χ 2 r ∼1) have a rather large elastic thickness, Te ∼30 km, a relatively low viscosity, η < 10 19 Pa.s, and a channel thicknesses L > 20 km. In this case, a clear trade-off between Te and η shows that a thin elastic lid (Te ∼ 10 km) requires a viscosity one order of magnitude lower than for equivalent elastic thicknesses Te ∼30 km (Fig. 6b) to maintain a comparable fit to the data (normalized χ 2 r ∼1). Tso the model has: η = 1.8 × 10 18 Pa.s, L = 9 km and Te = 2 km. For Zhari Nam Tso it has: η = 5 × 10 19 Pa.s, L = 6.4 km, Te = 5.9 km and t h = 5.23 ky. Given the small elastic thicknesses of these models, the glaciers have insignificant impact on the paleoshorelines distortion. The comparison between the measured  Supplementary Fig. D1. The load is cylindrical and increases linearly from 16 ka to 12 ka before present. It reaches a maximum of 140 m of equivalent water thickness, stays constant from 12 ka to 8 ka and decreases linearly to 0 from 8 ka to present. The radius (R) of the load is 60 km. Surface displacements relative to the horizontal initial stage at 16 ka are calculated assuming an elastic lid with thickness of 6 km overlying a viscoelastic channel with thickness of 6 km and viscosity of η = 10 18 Pa.s, η = 10 19 Pa.s, or η = 10 20 Pa.s (from left to right respectively).
(a) Surface vertical displacement along a radial cross section starting at the center of the cylindrical load calculated at 8 ka, 4 ka and at present. (b) Difference between the displacements at 8 ka and 4 ka (light blue) and between 8 ka and present (dark blue).
and predicted elevations of paleoshorelines shows in fact a rather poor fit with misfits mostly larger than the estimated 1.5 m uncertainty on the shoreline elevations as reflected the large reduced Chi-squares values ( Supplementary Fig. E1). These models do however predict distortions of the paleoshorelines with a short wavelength consistent with the observations, though of lower amplitude (for Zhari Nam Tso particularly). These models also predict subsidence at the center of the paleolakes, where the surface unloading is maximum. So these models are able to reproduce qualitatively some key features of the data but they fall short of fitting them quantitatively. Possible causes for the misfits include: DEMs error with a correlation scale larger than 10 km; incorrect hypothesis that the highest preserved shorelines around each lake are synchronous and were initially horizontal; incorrect model due to spatial variations of visco-elastic properties; incorrect load estimate because of changes of the topography and redistribution of mass by erosion or sedimentation during the regression of the lakes; sub-crustal deformation. Whatever the cause, the constraints on the rheology of the crust derived above probably hold anyway.

Significance of the rigid base boundary condition
The boundary condition in these calculations is a rigid base, implying support from the medium below the viscoelastic channel as discussed in England et al. (2013). As a result, the shoreline distortions become negligible when the channel thickness tends to 0 km. The best fitting models include this end-member and all require some minimum coupling between the upper elastic lid and the rigid base at the Holocene time scale. Assuming a rigid base is not realistic and subcrustal deformation could in fact have affected distortion of the paleoshorelines. We have therefore carried out forward tests to evaluate the possible contribution of this effect, using the viscoelastic code from Bills et al. (1994) which allows modeling a response of any vertically stratified viscoelastic Earth. Fig. 7 compares the surface deformation in the case of a subcrustal viscoelastic support to the rigid base approximation. The models both assume an elastic lid of thickness Te = 2 km overlying an L = 9 km thick viscoelastic channel of viscosity η = 1.8 × 10 18 Pa.s. One model has a quasi-rigid base (modeled here with a viscosity of 10 25 Pa.s) (Fig. 7a) and the other a 35 km thick sub-crustal elastic lid over a viscoelastic half-space of viscosity 10 18 Pa.s (Fig. 7b). The difference between these two models (Figs. 7c and 7d) is a low amplitude (<2 m) and long-wavelength (>100 km) signal. The observation of limited deformation at the scale of the paleoshorelines footprint thus requires coupling of the upper crust with a strong sub-crustal viscoelastic lid at the millennial time scale. We haven't explored further the constraints placed on subcrustal rheology as there is probably a large trade off between viscous and elastic support which cannot be easily resolved with the data considered in this study.
If criterions of fit were the distortion range as in England et al. (2013), there would be a strong trade-off between the channel properties and substrate viscosity. A low substrate viscosity of 10 21 Pa.s or less would then imply a higher channel viscosity than the values inferred from our inversions. However such models would not be able to produce the short-wavelength deformation of the paleoshorelines. This trade-off is much reduced if, as done in this study, the difference between the observed and predicted elevations at the actual location of the paleoshorelines measurements is minimized. Fig. 9. Synthetic test demonstrating how the apparent time dependent rheology deduced from the observations can result from a biviscous Burger rheology. We calculated the ground deformation due to a 100 km wide cylindrical time-varying load using the viscoelastic code from Bills et al. (1994). The model assumes a 5 km elastic lid over a 60 km thick viscoelastic layer with a transient viscosity of 10 18 Pa.s (dashed line) and a long-term viscosity of 10 20 Pa.s, overlying an elastic half space. Three synthetic datasets corresponding to either a long-term (∼ Zhari Nam Tso), a mid-term (∼ Siling Tso) or a short term (∼ present-day Siling Tso) scenario were produced. Two scenarios mimic the post Late Glacial Maximum history of lake transgression and regression observed at Zhari Nam Tso and Siling Tso. It assumes a highstand from 12 ka to 8 ka (similar to Zhari Nam Tso) or from 10 ka to 4 ka (similar to Siling Tso). The other loading history mimics the recent transgression of lake Siling Tso. It assumes a transgression of 10 m over 10 yr (present-day Siling Tso). The synthetic displacements are then inverted using the same methodology as the one used to invert the real observations. Results of the inversion of the long-term (a), the mid-term (b) and the short-term (c) scenarios are shown as 2 dimensional slices into the 3 parameters space for different ranges of elastic thickness.

Trade-off between viscous and elastic support of surface loads
The short-term deformation response to Siling Tso lake level variations can be fitted equally well with dominantly either elastic or viscous support. This is clearly seen in the trade-off between the viscosity η and the elastic thickness Te (Fig. 6b). By contrast, no such trade-off is seen in the results from the inversion of the paleoshorelines (Figs. 5 and 6). Some insight is gained by considering purely elastic models. A simple model consisting of an elastic lid over an inviscid fluid could explain the insignificant deformation of the shorelines as England et al. (2013) have found. However, the small distortion (<6 m) of the paleoshorelines requires Te ∼35 km. Such a large elastic thickness is considerably higher than the Te ∼ 10 km estimate deduced from gravity in the studied area (Braitenberg et al., 2003). It is therefore likely that, at the millennial time scale of the paleoshorelines deformation, surface loads are at least partially supported by viscous stresses.

Origin of the central downward distortion and upward bulge
We investigate here the mechanism responsible for the short wavelength deformation of the paleoshorelines and the lower elevations of the paleoshorelines nearer to the center of the lakes despite being closer to the maximum unloading. These two nonintuitive features are actually seen in the deformation pattern predicted by the best fitting models (Figs. 2 and 3). To gain insight we analyze the viscoelastic response predicted by our model in the simple case of a single cylindrical load with a load history similar to the one assigned to the Zhari Nam Tso (see simulations in Appendix D.1). Fig. 8 illustrates how viscous effects influence the wavelength, sign and amplitude of the vertical displacements in space and time. The lake transgression induces a central zone of subsidence fringed by a zone of uplift that diffuses away as a result of channel flow. The opposite happens during regression and the two patterns interfere. For large enough channel viscosities and a progressive unloading, the subsidence induced by the lake highstand goes on during the early stage of surface unloading and can Fig. 10. Synthetic test demonstrating how the apparent time-dependent rheology deduced from the observations can result from depth variations of viscosity. We calculated the ground deformation due to a 100 km wide cylindrical time-varying load using the viscoelastic code from Bills et al. (1994). The model assumes a 5 km elastic lid over a 60 km thick stratified viscoelastic body overlying an elastic half-space. Crustal viscosities vary between 10 18 Pa.s and 10 21 Pa.s. The minimum viscosity is a mid-crustal depth where the temperature is presumably maximum (Wang et al., 2013). Three synthetic datasets corresponding to either a long-term (∼ Zhari Nam Tso), a mid-term (∼ Siling Tso) or a short term (∼ present-day Siling Tso) scenario were produced as in Fig. 9, and inverted using the same methodology as the one used to invert the real observations. Results of the inversion of the long-term (a), the mid-term (b) and the short-term (c) scenarios are shown as 2 dimensional slices into the 3 parameters space for different ranges of elastic thickness. dominate the uplift induced by the unloading. The model predicts a transition from subsidence to uplift after some time that scales with the viscoelastic relaxation time. It takes then time before uplift compensates the cumulated initial subsidence. Note that this mechanism does not imply that, at present, the ground would be still subsiding near the center of the paleolakes. This effect doesn't happen in the case of an abrupt lake regression as assumed in the study of England et al. (2013). Uplift then starts directly after unloading. The observation of a downward distortion of the paleoshorelines near the center of Siling Tso and Zhari Nam Tso (Fig. 2) thus suggests a long enough relaxation time that it allows the surface to continue subsiding after the lake started regressing.

Reconciliation of results inferred from decadal and millennial time scales
The analysis of Doin et al. (2015) and our own modeling results ( Fig. 6) suggest an equivalent elastic thickness Te ∼30 km at decadal time scale larger than our estimate of Te < 5-10 km derived at the millennial time scale from the paleoshorelines (Figs. 5 and 6 and Supplementary Fig. E.2). We also find that, at the decadal time scale, a viscosity of ∼10 18 Pa.s lower than at the millennial time scale. The domains of acceptable model parameters do not overlap (Figs. 5 and 6a relative to 6b). The effective rheology of the Tibetan crust, when represented by an elastic lid over a finite viscoelastic medium thus appears time-dependent. This behavior could reflect that the effective viscosity is actually time-dependent as could happen for example with a non-linear (stress-dependent) rheology or a Burgers body rheology. A Burgers body, which consists of a Maxwell element in series with a Kelvin element, is a common model used in postseismic studies (e.g., Bürgmann and Dresen, 2008). Such a model has a short-term viscosity, associated to the Kelvin element, and a long-term viscosity associated to the Maxwell element. A depth-varying viscosity would also imply a time-dependent effective viscosity in response to a surface load. We therefore test if these alternative models can reconcile the results obtained at the decadal and millennial time-scale. In these calculations we use the modeling approach of Bills et al. (1994), which allows different short-term and long-term viscosities and depth variations.
We produce synthetic data at the decadal and millennial time scale with simplified loading history similar to the ones assumed in the actual data analysis. We next invert the synthetic data us-ing the same inversion procedure as above. A cylindrical load of 50 km in radius and a height of 140 m, 60 m and 10 m were chosen to represent the lake load for the millennial (Zhari Nam Tso "long-term" and Siling Tso "mid-term") and decadal (Siling Tso "short-term") cases respectively. We assume two long-term load histories. A 'long term' one is similar to that inferred at Zhari Nam Tso: the load corresponding to the highstand is applied from 12 ka to 8 ka and removed afterwards. The "mid-term" one corresponds to Siling Tso: the load is applied from 10 ka to 4 ka and removed afterwards. At the decadal timescale we assume 10 m of transgression over 10 yrs, which is similar to the Siling Tso case (Doin et al., 2015). For a given viscosity profile, we generate three synthetic datasets representing the three time scales, which are then inverted to find the best equivalent simple viscoelastic model, consisting of an elastic lid over a viscous channel as assumed above. We test a Burgers body model (Figs. 9). In that case two Maxwell elements in parallel are used to produce an equivalent Burgers body rheology (Müller, 1986). We also test a multilayered model (Figs. 10).
The model of Fig. 9 contains a 60 km thick viscoelastic layer with a transient viscosity of 10 18 Pa.s and a steady state viscosity of 10 20 Pa.s. The multilayered model of Fig. 10, assumes a Maxwell rheology in each layer. The model assumes lower viscosity at mid-crustal depth where the temperature has presumably a local maximum due to the temperature inversion caused by the underthrusting of India beneath Tibet (e.g., Wang et al., 2013). The inversion of the synthetic data generated with both models yield low elastic thicknesses at the millennial timescales and higher elastic thicknesses at the decadal cases. This is particularly clear for the stratified model (Fig. 10). At both "long-term" and "mid-term" timescales the equivalent elastic thicknesses is <10 km whereas for a "short-term" decadal timescale the equivalent elastic thickness is >20 km. The apparent viscosities vary significantly.
Compared to the "long term" case, it is lower by a factor 2 in the "mid-term" and by a factor 10 in the "short term" (Fig. 10). These trends are similar to those observed in the inversion of the Siling Tso and Zhari Nam Tso data. We note that the inversion of synthetic data generated with the stratified model reproduces better the trade-off between η, Te and L at the decadal and millennial time scales.

Conclusion
We used the deformation response to lake level variations at the decadal and millennial time scale to place constraints on the rheology of the Tibetan crust. Our study confirms and expands the results of the previous studies presented by England et al. (2013), Shi et al. (2015) and Doin et al. (2015). Compared to the Lake Bonneville archetype (Bills and May, 1987;Nakiboglu and Lambeck, 1982), the paleoshorelines around the Tibetan lakes studied here show much smaller distortions, despite a comparable lake regression, and a downward deflection of the centers of the lakes instead of an upward deflection. The downward deflection is a signature of channel flow in the crust, but the upper crust must have remained well coupled, at the millennial time scale, to a quasi-rigid sub-crustal lid to explain the small distortion of the shorelines. Reconciling the millennial and decadal deformation response can be achieved with either vertical layering or a non-linear rheology. For example, a Burgers body rheology with a transient viscosity of 10 18 Pa.s and a long term viscosity of 10 20 Pa.s could reproduce the deformation response at both the decadal and the millennial time scale. An example of an alternatively layered model has a lower viscosity (10 18 Pa.s) mid-crust (between ∼10 and 30 km depth) embedded in a higher viscosity crust (>10 20 Pa.s). A viscosity <10 18 Pa.s, as proposed in some studies, could only exist in a thin layer (<10 km). The crustal rheology derived from this study is consistent with the low effective elastic thickness of Tibet, Te ∼ 8 km, derived from gravity studies (Braitenberg et al., 2003). Indeed the short-wavelength deformation of the paleoshorelines requires a small Te (<5 km) of the upper-crustal lid and some degree of relaxation through midcrustal channel flow at the millennial time scale. Only the thin upper crustal lid would be able to support topographic loads at the much longer geological time scale associated with the evolution of topography.