A bound on the viscosity of the Tibetan crust from the horizontality of palaeolake shorelines

Palaeoshorelines around four large lakes in central Tibet record a latest-Pleistocene-to-Holocene high stand during which the lakes were ﬁ lled 150 – 200 m more deeply than they are at present. GPS measurements of shoreline elevations around Zhari Namtso show that they are horizontal to within 7 2 m at the 2 (cid:1) s level. Measurements of height made by combining Shuttle Radar Topographic Mission elevations with Google Earth imagery of shorelines around Zhari Namtso, Tangra Yumtso, Taro Tso, and Ngangla Ringtso show that all the palaeoshorelines are horizontal within measurement uncertainty. Support of the lake loads by elastic stresses can explain the horizontality of the shorelines only if the equivalent elastic thickness of the crust exceeds 15 – 25 km. The observations are more plausibly explained by support of the lake loads through viscous stresses in the middle to lower crust. This support requires that the viscosity of the middle to lower crust is at least 10 19 – 10 20 Pa s. These values are consistent with estimates from studies of post-seismic relaxation after large earthquakes of the region and are higher, by two orders of magnitude, than would permit signi ﬁ cant lateral ﬂ ux of material through a channel in the middle to lower crust. & 2013 The Authors. Published by Elsevier B.V.


Introduction
There is little agreement about the strength of the continents, beyond the obvious fact that the lithosphere as a whole is strong enough to support the lateral variations in density associated with structures such as mountain ranges (e.g. England and Houseman, 1986;Molnar, 1991, 1997;Molnar, 1991;Molnar and Lyon-Caen, 1988). This strength is variously suggested to lie within the upper crust, the crust as a whole, the lithospheric mantle, or combinations thereof (e.g. Brace and Kohlstedt, 1980;Bürgmann and Dresen, 2008;Burov and Watts, 2006;Chen and Molnar, 1983;Jackson, 2002;Jackson et al., 2008;Maggi et al., 2000;Molnar, 1991;Sonder and England, 1986). It has been suggested that parts of the crust beneath the Tibetan plateau form a low-viscosity (∼10 16 Pa s) channel that accommodates significant horizontal flux of material (e.g. Clark et al., 2005;Clark and Royden, 2000;Royden et al., 1997;Morgan, 1985, 1987). If this suggestion is correct, then the deformation of the upper crust is qualitatively different from deformation of the lower lithosphere and deduction, from surface observations, of forces acting on the lithosphere as a whole will be difficult or impossible (e.g. Molnar, 1988). In contrast, estimates based on post-seismic deformation following large earthquakes suggest that the viscosity remains relatively high (≳10 18 Pa s) throughout the crust in Tibet (Hilley et al., 2005(Hilley et al., , 2009Ryder et al., 2011;Wen et al., 2012;Yamasaki and Houseman, 2012). Here we present an analysis of the response of the Tibetan crust to loading by palaeolakes, on the time scale of kiloyears, which supports the inferences drawn from post-seismic relaxation on the time scale of years.
Google Earth imagery with data from the Shuttle Radar Topographic Mission (SRTM, Farr et al., 2007). Close agreement between these two techniques allows us to extend our analysis, using the remotely sensed data alone, to the three other large palaeolakes in the region: those surrounding Tangra Yumtso, Taro Tso and Ngangla Ringtso (Figs. 1 and 3).

Palaeoshorelines of Zhari Namtso
A well-preserved suite of palaeoshorelines surrounds Zhari Namtso, from the level of the lake itself to a height of about 140 m above the present-day water level (e.g. Fig. 2). We inspected these shorelines in the field at the locations shown in Fig. 3. On slopes with low gradient (≲101), such as those found along the eastern margin of the lake and in small inlets elsewhere, shorelines are typically preserved as ridges of shingle, which form barriers to drainage with marshlands formed between them. On steep slopes (≳301), shorelines are usually preserved as notches cut into bedrock, though the precise position of the break in slope may be obscured by outwash. Many of the shorelinesincluding the highest levelappear to form couplets, separated in height by 0.5-1 m, which may indicate a seasonal variability in water height; this inferred seasonality is not apparent on steeper slopes where the shoreline is preserved as a single erosional notch.

Precise elevations of palaeoshorelines
Elevations of the highest palaeoshorelines were measured using differential GPS. Measurements were made during two weeks in September 2006, using two Ashtech Promark II GPS receivers, one acting as a fixed base station and the other as a kinematic rover. The base stations were sited on bedrock, where available, or on large boulders firmly embedded in sediment, and the maximum horizontal distance between successive stations was 25 km. The base stations were leapfrogged by occupying simultaneously, at the end of each day, the current day's base station and the base station chosen for the following day. L1 phase measurements, obtained at 10-s intervals, were processed using the RTKLIB version 2.4.1 (Takasu, 2011), using fixed IGS final orbits and fixed Earth orientation parameters from IERS Bulletin B. Estimates of total ionospheric electron contents from global processing were used to correct for delays to the L1 phase, and the delay from the dry atmosphere was calculated from the parameters of Saastamoinen (1972). The day-to-day repeatability of base-station locations was within ∼1 m. The GPS processing delivered heights above the reference ellipsoid, which were converted into orthometric heights by subtracting at each elevation the interpolated geoid height from EGM2008 (Pavlis et al., 2008).
Profiles across palaeoshorelines were generally measured on slopes between 51 and 251 because the ridges of shingle formed by shorelines on more gentle slopes were sometimes difficult to identify in the field, and the shorelines on steeper slopes were often obliterated by colluvium. At each measurement site a profile was recorded along the base of the shoreline notch and one to three transects were recorded perpendicular to the shoreline. The heights reported here are those of the base of the shoreline notch, cross-checked against changes in slope on the transects. The locations of the stations are shown in Fig. 3a, and in the accompanying Google Earth KMZ file (Supplementary material).
The elevations of all the highstand shorelines we visited are shown in Fig. 4a; the mean elevation of the shorelines is 4750.7 m, with a standard deviation of 1.1 m; no highstand shoreline lies outside the range from 4748 m to 4753 m.

Elevations of palaeoshorelines derived from remote sensing
We use Google Earth satellite imagery and digital elevation data from the Shuttle Radar Topographic Mission (SRTM) to extend our data set. The palaeoshorelines around Zhari Namtso that we  (Farr et al., 2007). Present lakes are shown in blue; their maximum extents at the time of high stands discussed here are shown in cyan. NR: Ngangla Ringtso; TT: Taro Tso; ZN: Zhari Namtso; and TY: Tangra Yumtso. Boxes show extents of Fig. 3a and b.
In the inset, the red box shows the location of the main figure; L, S and N mark the positions of other lakes referred to in the text: Longmu Tso and Sumxi Tso (Kong et al., 2007), Siling Tso (Shi et al., 2012) and Nam Tso (Wallis et al., 2010) respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) visited in the field are clearly visible on the satellite images; multiple palaeoshorelines show as lineations that can be followed laterally for a kilometre or more. The tops of these flights of shorelines are marked by a contrast in the colour and texture of the surface, corresponding to the sharp transition between areas that were submerged during highstand conditions, which now preserve flights of regressional shorelines, often with exposures of lake sediment and marshland between the beach ridges, and areas that were above the maximum lake level, where exposures consist of bedrock or alluvium with no modification by wave action (see Fig. 2). Each of the locations where we measured the lake shorelines in the field coincides, to within about 10 m in the horizontal, with such a visual marker on the satellite images.
The sharp transition marking the top of the flight of shorelines is visible at many places around Zhari Namtso, even in regions of low surface slope where the shorelines were difficult to pick unambiguously on the ground. We therefore investigated the feasibility of measuring shoreline heights by picking geographical locations from the satellite images, and determining the elevations of those locations from SRTM data. In order to minimize height errors due to errors in horizontal position, and because height errors are smaller in regions of low relief (Farr et al., 2007), we restricted our picks to regions where the magnitude of the local surface slope was lower than 51, identifying 58 such places. We fit bi-cubic surfaces to the SRTM elevation data within 300 m of each of these 58 points, to estimate their elevations. To avoid errors resulting from possible short-wavelength noise in the SRTM data, we rejected any solution in which the RMS misfit to the smooth surface exceeded 3 m; this resulted in the loss of fewer than 5% of the points in each of the lakes studied. The distribution of elevations of the palaeoshorelines of Zhari Namtso estimated in this way is shown in Fig. 4b; their mean value is 4750.8 m, with a standard deviation of 1.0 m. The agreement between this mean and the mean elevation of shorelines measured by kinematic GPS is coincidental; differences are small in comparison with the uncertainties of either data set.
We used the same remote-sensing techniques to measure the heights of well-preserved sets of palaeoshorelines, similar in appearance to those seen around Zhari Namtso, around the three other lake systems shown in Fig. 1: Tangra Yumtso (Fig. 3a), Ngangla Ringtso (Fig. 3b), and the lakes of Taro Tso and Zabuye Tso, which are surrounded by a single highstand shoreline (Fig. 3b). The topographic slopes around Tangra Yumtso are steeper than around the other lakes and we used a cut-off of 101 for the slope of surfaces on which we measured heights; otherwise, the procedures were the same. The heights of these shorelines were Tangra Yumtso, 4740.6 71.8 m (106 sites, Fig. 5a); At first sight, the small scatter in heights of palaeoshorelines is surprising, because the nominal accuracy of SRTM data is between 6 and 8 m (Farr et al., 2007). The mean elevation errors in the SRTM are, however, usually estimated at a continent-by-continent level (e.g. Farr et al., 2007;Rodríguez et al., 2006), and absolute errors vary considerably from place to place. Comparison of SRTM data in Tibet with elevations derived from kinematic GPS indicate that the errors there are as low as 2 m (Rodríguez et al., 2005, Fig. 3), and the error in height difference between two locations within a couple of hundred kilometres from one another is smaller than the absolute height error (Rodríguez et al., 2006). Because our measurements are restricted to locations where the slopes are less than 51, the errors should be lower than the regional averages (Farr et al., 2007); in addition, we have used a smooth surface fitted through 49 elevation measurements in estimating the height of each point, reducing errors due to local noise in the DEM. A Google Earth KMZ file (Supplementary material) shows the locations of all the points measured, and contours corresponding to the mean elevation of each shoreline, and to elevations 4 m above and below ( 72 m in the case of Zhari Namtso). (a) Zhari Namtso and Tangra Yumtso. Larger yellow circles show the locations at which the elevations of the highest palaeoshorelines around Zhari Namtso were measured using kinematic GPS. Smaller circles show the locations at which the elevations of the highest palaeoshorelines around Zhari Namtso and Tangra Yumtso were measured using Google Earth imagery and the SRTM data set. The town of Coqên is indicated by the letter C. The mean elevation of the highstand palaeoshorelines around Zhari Namtso is 4751 m (Fig. 4a); the mean height for Tangra Yumtso is 4741 m (Fig. 5a). (b) As (a), for Ngangla Ringtso, Taro Tso and Zabuye Tso. The mean height of the palaeoshorelines around Taro Tso and Zabuye Tso is 4606 m (Fig. 5b) and, of those around Ngangla Ringtso, 4864 m ( Fig. 5c). Colour scale for topography is as in Fig. 1. The additional water depth at the time of high stand is the difference between the height of the palaeoshorelines at each lake, and the present elevation of the land, or lake, surface lower than those shorelines. White lines show shorelines of present lakes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Loading by the palaeolakes
We used the positions of the high shorelines discussed above to quantify the extent of the palaeolakes and the loads that they represented (Fig. 3). If the crust of Tibet responded to the loads in a purely elastic manner, then the timing of loading and unloading would be immaterial; only the shapes and sizes of the loads would matter (Section 3.1). If, however, the loads were supported in part by viscous stresses, then the amount of deformation would depend on the time interval over which the load was emplaced, and the degree to which that deformation is detectable at the present day would depend upon the time elapsed since the loads were removed.
Cosmogenic dating of palaeoshoreline deposits corresponding to the 4751-m highstand level of Zhari Namtso yields ages of 14 kyr ( 10 Be) and 15 kyr ( 26 Al) (Kong et al., 2011). (The elevation of this site is given as 4761 m by Kong et al. (2011), but comparison of their site photograph with Google Earth imagery and SRTM data suggests that an elevation of ∼4751 m is more appropriate, see KMZ file.) A shoreline at 4718 m, 100 m above the present level of Zhari Namtso, yields ages of 27 kyr ( 10 Be) and 24 kyr ( 26 )Al) (Kong et al., 2011). The prominent highstand shoreline at 4741 m, around Tangra Yumtso also encircles Tangqung Tso, a small lake that is connected to Tangra Yumtso by a pass at 4610 m. Deposits from shorelines near the highstand at 4741 m around Tangqung Tso yield ages between 20 and 23 kyr with both 10 Be and 26 Al (Kong et al., 2011). Two 14 C dates on shorelines 20 m below the highstand above Zabuye Tso yield ages of 11.2 and 9.9 kyr (Hudson and Quade, 2013), while shorelines 100 m below the highstand record 14 C ages of 22-24 kyr (Qi and Zheng, 1995;Zheng et al., 1989, cited by Yu et al., 2001). Sedimentary and palaeoenvironmental records from Zabuye Tso are consistent with the dates for the shorelines; these records suggest a period of lake filling began at about 28 kybp, with a strong phase of replenishment at 23-22 kybp, and a second period of high lake level, probably relating to increases in glacial melt water following the last glacial maximum, lasting from 16 to 11 ka (Wang et al., 2002). Seven 14 C dates on the highstand shorelines around Ngangla Ringtso yield ages of 10.4-8.6 kyr, while two dates from shorelines 40 m below the highstand give 7.9 and 8.1 kyr (Hudson and Quade, 2013).
A number of lines of evidence indicate that the lake levels we investigate began to drop in middle-Holocene time. The palaeoenvironmental record from Zabuye Tso suggests that the lake began to shrink some time after 10.6 ka, with the most rapid desiccation beginning at 5.6 kybp (Wang et al., 2002), and the lake level in Tangra Yumtso seems to have declined rapidly between about 7 ka and 2 ka (Hao et al., 2012). To the west of the area we consider studies of sediment cores from, and geomorphology around, Sumxi Tso and Longmu Tso (Fig. 1) show periods of high water level from about 11 kybp, terminating in a high stand between 7.5 and 6 kybp, following which the lake level dropped abruptly in perhaps a few centuries (Avouac et al., 1996;Gasse et al., 1991;Kong et al., 2007).
The influences of climate, catchment area, and topography, all of which vary across the Tibetan plateau, preclude construction, from the existing data, of a detailed loading and unloading history for the lakes we study. The conclusions to the analysis that follows are, however, insensitive to the detail of that history. Below we show that our estimates of the viscosity of the Tibetan crust depend either on the duration of the lake loading or on the time that has elapsed since the loading ended. It seems likely from the available data that the duration of the most recent phase of lake high stand was of order 3-10 kyr, although there may well have been fluctuations in levels during that time (e.g. Wang et al., 2002). The decline in lake levels appears to have begun in early-to-mid-Holocene time, and our conclusions do not depend strongly on whether that decline was abrupt, as inferred by Avouac et al. (1996) for Longmu Tso, or took place gradually over the past 5-10 kyr. There is evidence for an earlier phase of lake filling from poorly preserved shorelines at higher elevations than those we have measured. Shorelines close to 4900 m around Tangra Yumtso yield ages of about 220 kyr (Kong et al., 2011), consistent with their poorer state of preservation. Comparable higher and poorly preserved shorelines are also seen above Taro Tso and Zabuye Tso, to about 40 m above the level of the prominent highstand at 4605 m, and a remnant shoreline above the 4864-m highstand at Ngangla Ringtso has been dated, using U-Th, to 211 kyr (Hudson and Quade, 2013). These higher shorelines lie at levels from which the lakes could drain to the north, and form part of a much larger system of lakes. Because this loading took place about 200 kyr before the loading we consider here, it has no influence on our analysis.

Analysis
The palaeoshorelines of all the lakes we investigate here are horizontal to within 74 m at the 2Às level when measured using remote-sensing techniques (Section 2.3, Fig. 5), and to within 72 m at the 2Às level when heights are measured using kinematic GPS (Section 2.2, Fig. 4). In this section we compare the observed shoreline elevations (Figs. 4 and 5) with the calculated response of the earth's surface to loading and unloading by raising and lowering of the lake levels. The elevations of the palaeoshorelines carry no information about the absolute vertical motion of the land surface; they record only the relative vertical motion of one part of the land surface with respect to any other part, which we refer to as distortion.
The palaeoshorelines must have been horizontal when they formed, and there are two classes of explanation for the observation that they are, to within measurement uncertainty, horizontal now. Either the land surface at the locations of shorelines was not measurably distorted by the loading of the lakes, or measurable distortion did occur, but insufficient time has elapsed to allow the crust to recover from that distortion. We analyse these possibilities using a physical model consisting of an elastic upper crust overlying a visco-elastic lower crust. We use the expressions of Nakiboğlu and Lambeck (1982) for a cylindrical load, which allow us both to investigate the scaling, by approximating the lakes as cylindrical loads of the appropriate size, and to perform detailed calculations by discretising the lakes into many separate sub-loads and superposing the solutions (Appendix).
We note immediately that, in the case of a perfectly axisymmetric lake on a flat surface, any shoreline formed during the filling or emptying of the lake would be at a constant distance from the centre of the load and therefore, although it would move vertically in response to the varying load, at no time would it be distorted. Even in the case of non-circular loads, the distortion of shorelines may be small if their shapes resemble the shapes of the loads themselves-as may occur if steep topography forces the shape of the lake to vary little between its highstand and lowstand phases. This consideration alone renders most of the smaller lakes in Tibet unsuitable for our analysis, and may also apply to larger lakes such as Nam Tso, whose palaeoshorelines lie close to their present shores.

Elastic support of lake loads
We first consider the case of an elastic lid of thickness h c that rests on an inviscid foundation (Appendix, Eq. (4)). In this case, the response of the earth's surface to loading and unloading is instantaneous, and the distortion of the palaeoshorelines at the present day would be equal and opposite to the distortion of the land surface at the time of loading. Illustrative calculations (Figs. 6a, 7a, and 8a,d) show depression of the land surface caused by all the lake loads, discretised as described in the Appendix, acting simultaneously on a lid of thickness 10 km. It can be seen that displacements of the palaeoshorelines are determined primarily by the loading of their own lakes, but the neighbouring lakes have secondary influences, particularly in the region between Zhari Namtso and Tangra Yumtso (Fig. 6). We allow the thickness of the elastic lid to vary between 5 and 30 km in steps of 1 km and calculate the depression of the earth's surface caused by these loads, for each site at which we have measured the elevation of a palaeoshoreline. Reversal of this depression following unloading would result in a range of displacement of the sites of palaeoshorelines which we show, lake by lake, in Figs. 6b, 7b, and 8b,e. For example, if the elastic lid beneath Zhari Namtso had a thickness h c ¼10 km, the calculated displacements for the sites where we have measured the palaeoshorelines would lie between +5 and +15 m (Fig. 6a, and marked by red lines in Fig. 6b), so the observable distortion of the palaeoshoreline would be 10 m.
The observed distortion of the palaeoshorelines heights around Zhari Namtso is 72 m at the 2Às level. If this observed range of heights is to be explained by elastic support of the upper crust, then the thickness of this elastic layer must be at least 21 km: the thickness at which the calculated distortion falls to the same level as the observed range (red error bar in Fig. 6b). The distributions of shoreline heights around the other three lakes yield similar conclusions. The thickness of the elastic lid required to explain the shoreline elevations around Tangra Yumtso (Fig. 5a) is 14 km if the shorelines are horizontal to within 74 m (the 2Às level of the remote-sensing measurements), and 22 km if they are horizontal to within 72 m, the 1-s level (Fig. 7). The equivalent thicknesses estimated from shorelines around Taro Tso, and around Ngangla Ringtso (Fig. 5b and c) are 14 and 12 km at the 2Às level, or 25 and 19 km at the 1 À s level (Fig. 8).

Viscous support of lake loads
The alternative mechanism for supporting the lake loads without leading to measurable distortion of the shorelines is viscous stress associated with flow in the crust beneath the elastic lid. This flow takes place on a time scale, τ, that is governed primarily by the radius, A, of the load, and the viscosity, η, of the medium (see Appendix). If the loading phase was much shorter than the viscous response time, τ, then the depression of the earth's surface would have been small in comparison with that attained at full isostatic compensation; equally, if the time since the load was removed is much smaller than τ, then insufficient time will have elapsed to allow distortion of the shorelines established during loading. Note that the Maxwell time (ratio of viscosity to shear modulus) is shorter than ∼300 yr for the viscosities that concern us here, so τ (Eq. (1)) is the time scale that dominates the problem.
We make a first-order estimate of the viscosity by approximating the lakes as cylindrical loads whose effective radii are ∼30 km (Figs. 6-8). The time interval over which the lakes were close to their maximum loading was in the range from ∼3000 yr to ∼10;000 yr (Section 2.4). If we set a minimum bound on the viscosity by assuming that the crust was too viscous to allow measurable distortion during the loading phase, then τ≳3000 yr, and η≳3 Â 10 19 Pa s. By the same argument, if we suppose that the lack of distortion of shorelines arises because the viscous response time is longer than the time since the end of the high stand (∼3-10 kyr, see Section 2.4), then the viscosity of the lower crust must be ≳3 Â 10 19 −10 20 Pa s. In what follows, we assume that the unloading took place instantaneously at 5 kybp but very similar results are obtained if the age of unloading is assumed to be 10 kyr.
We now expand our analysis to determine how the conclusions drawn from the scaling arguments are modified by combining an elastic lid with a viscoelastic layer. We calculate (see Appendix) the distortions of the palaeoshorelines for the lakes under the following conditions. The lakes are assumed to have unloaded instantaneously at 5 kybp, having loaded instantaneously at 10 kyr, 3 kyr, or 1 kyr before that time. We allow the thickness of the elastic lid to vary between 5 and 25 km. The lid rests on a visco-elastic half space; the influence of replacing the half space by a layer of finite thickness is discussed in Section 3.3. The shear elastic modulus is taken to be 3 Â 10 10 Pa throughout, and the viscosity of the half space is allowed to vary between 10 17 and 10 22 Pa s. For all calculations, we use the distributions of loads shown in Fig. 3, discretised as described in the Appendix, and calculate the range of heights that would be observed at the locations of measured palaeoshorelines, for the given combinations of thickness of the lid and viscosity of its substrate.
These calculations show that two distinct regimes lead to undistorted shorelines under the range of loading conditions that seem likely for the late-Pleistocene-Holocene lakes in central Tibet. The first occurs, as would be expected from the calculations for an elastic lid on an inviscid foundation (Figs. 6-8), when the thickness of the elastic lid exceeds a critical value, above which the loads are supported by elastic stresses. In the second regime the lid is below that critical thickness and the small observed range in shoreline heights requires that a significant fraction of the load is supported by viscous stresses in the substrate. The transition between these regimes is sharp: if the elastic lid is even a couple of kilometres thinner than its critical thickness, then a highviscosity substrate is required to support the loads. For Zhari Namtso this critical thickness is 21 km (Fig. 9a-c) and the required viscosity of the substrate is ≳2 Â 10 20 Pa s if the duration of loading is 10 kyr (Fig. 9a), dropping to ≳5 Â 10 19 Pa s if the duration of loading is 1 kyr (Fig. 9c). Equally sharp transitions take place at elastic lid thicknesses of 14 km (2Às level) or 22 km (1 Às level) for Tangra Yumtso (Fig. 9d), at 14 or 25 km for Taro Tso (Fig. 9e) and at 12 or 19 km for Ngangla Ringtso (Fig. 9f). As would be expected from the scaling relation (Eq. (1), and Appendix) the viscosities required to explain the observed distributions of shoreline heights are comparable for all the lakes, because they depend primarily on the duration of loading.

Alternative lower boundary conditions
The analysis we have pursued so far treats the middle-to-lower crust as a half-space of constant viscosity. We investigate the sensitivity of our results to this assumption by restricting flow beneath the elastic lid to a layer of finite thickness, using two endmember boundary conditions for the base of that layer: zero velocity and zero shear traction. The rigid base could represent the uppermost mantle (e.g. Chen and Molnar, 1983), or a lower crust Layer Thickness (km) Fig. 10. Influence of the thickness of the viscous layer below the elastic lid, and its basal boundary condition, on estimate of its viscosity. Curves show viscosity of the layer that is required to keep the distortion of the palaeoshorelines around Zhari Namtso lower than 4 m as the thickness of the layer varies. The solid line corresponds to the case in which the shear traction is zero on the base of the viscous layer; the dashed line corresponds to the case of a rigid boundary condition on the base of the layer. The thickness of the elastic lid is 10 km, and the same loading history is used as in Fig. 9b. Black dot shows the value of viscosity deduced from calculations in which the viscous layer is a half-space (Fig. 9b).  Fig. 9. Range of displacements (equivalent to observable distortion) for sites where elevations of palaeoshorelines have been determined around the lakes shown in Fig. 1. The loads are applied to an elastic lid of thickness h c , overlying a viscoelastic half-space, of viscosity η; the shear modulus is 3 Â 10 10 Pa throughout. (a-c) Calculated ranges for Zhari Namtso; the loads shown in Fig. 3 are applied for 10 kyr (a), 3 kyr (b), and 1 kyr (c), and in all cases the load is removed 5 kyr before present. (d-f) As (b), for Tangra Yumtso, Taro Tso, and Ngangla Ringtso. Grey lines indicate the observed range of shoreline elevations for Zhari Namtso, measured by GPS ( 72 m, Section 2.2) and estimates of the ranges the other lakes, based on the remotely sensed measurements ( 7 2 to7 4 m, for Tangra Yumtso and Ngangla Ringtso; 7 1.5 to 7 3 m, for Taro Tso, Section 2.3). that is strong because of its mineralogy or low water content (e.g. Jackson et al., 2004); the free-slip lower boundary condition could represent middle or lower crust that has low viscosity, due to its high temperature or to the presence of partial melt (e.g. Chen and Molnar, 1983;Clark and Royden, 2000). Fig. 10 shows the influence of replacing the viscous half-space by a layer of finite thickness; calculations are shown for shoreline observations at Zhari Namtso but the results for other lakes are similar. If the viscous layer has zero traction at its base (corresponding to a lower viscosity layer beneath it), then it responds more rapidly to a given load than does a half-space of the same viscosity; in consequence a higher viscosity is required to account for horizontality of the shorelines. For example, if the middle-to-lower crust beneath Zhari Namtso is treated as a viscous layer with no shear traction at its base, the required viscosity is greater than about 3 Â 10 21 Pa s when the layer is 20 km thick, and greater than about 3 Â 10 20 Pa s if it is 50-60 km thick (dashed line in Fig. 10). The reciprocal argument applies when the base of the layer is rigid; in this case the lower bound on the viscosity required to support the loads of the palaeolakes is 4 Â 10 19 Pa s for a layer thickness of 60 km, 2 Â 10 19 Pa s for a thickness of 40 km, and 3 Â 10 18 Pa s for a thickness of 20 km. If the layer is thicker than about 100 km, estimated viscosities with either basal condition are essentially the same as those for the half-space.

Discussion and conclusions
Although the shorelines of the four palaeolakes are horizontal because their loads were supported by a combination of elastic stresses in the upper crust and viscous stresses in the middle and lower crust, the fraction of parameter space in which these two sources of support are comparable in magnitude is small (Fig. 9). It therefore makes sense to consider separately the conditions under which either elastic support or viscous support can account for the horizontality of the palaeoshorelines.

Support by elastic stresses
If the palaeoshorelines are horizontal because the effective elastic thickness of the crust was great enough to support the lake loads without measurable distortion of the surface, then that effective thickness is most tightly constrained by our GPS measurements around Zhari Namtso, which have a standard deviation of 1 m, thus a range of 72 m at the 2Às level. Those observations (Fig. 6) imply that the elastic lid beneath Zhari Namtso must be at least 21 km thick at the 95% confidence level: there is a 5% probability that the true range of shoreline heights is greater than 72 m, allowing an elastic thickness smaller than 21 km. Our measurements around the other lakes are less precise and, therefore, permit a lower equivalent elastic thicknesses. The range in shoreline heights of 74 m at the 2Às level around Tangra Yumtso requires an elastic thickness of at least 14 km (Fig. 7); the range of73 m at the 2Às level around Taro Tso and Zabuye Tso requires an elastic thickness of at least 14 km (Fig. 8b), and the range of74 m around Ngangla Ringtso requires a thickness of at least 12 km (Fig. 8e).
It may well be, however, that the actual range of shoreline heights around the lakes that we measured with remote-sensing techniques alone is no greater than the 72 m (4-m range) that we measured around Zhari Namtso with GPS (Fig. 4). This supposition is supported by two other sets of measurements of palaeoshorelines of large Tibetan lakes using GPS, which also report shoreline heights that are restricted to a range of 72 m (Siling Tso, Shi et al., 2012; and Nam Tso, Wallis et al., 2010; see Fig. 1 for locations). If we use a 4-m range for the shorelines around Tangra Yumtso, Taro Tso and Ngangla Ringtso, the required elastic thicknesses are, respectively, at least 22 km, 25 km, and 19 km (Figs. 7 and 8).
Estimates of the equivalent elastic thickness of the crust in central Tibet based on gravity anomalies are lower than 10 km (Braitenberg et al., 2003;Jordan and Watts, 2005), and consistent with a recent estimate of 7 km for the elastic thickness of the crust in eastern Tibet (Fielding and McKenzie, 2012). InSAR studies of the depth-extent of faulting in large earthquakes of the region constrain most of the moment release to lie in the top 10-12 km of the crust (Elliott et al., 2010;Lasserre et al., 2005;Li et al., 2011;Wang et al., 2007), and the depths of large earthquakes beneath the Tibetan plateau, as determined by long-period body wave modelling, are mostly shallower than 12 km (Chu et al., 2009;Molnar and Chen, 1983;Molnar and Lyon-Caen, 1989;Zhu et al., 2006). Recalling that the depth of earthquake rupture is likely to be larger than the long-term equivalent elastic thickness, both because the top few kilometres of the crust are probably too weak to contribute to the elastic strength, and because earthquakes may rupture downwards into material that creeps during the interseismic period, we conclude that these observations are all consistent with an equivalent elastic thickness for central Tibet that is less than ∼10 km. GPS data have been used elsewhere to estimate the locking depths of faults and hence to constrain the thickness of the seismogenic upper crust. Such a constraint is not possible here, because it would require ∼10 points within a distance of ∼π times the elastic thickness (e.g. Segall, 2010, p. 45), and there are fewer than 20 published GPS velocities within the whole of central Tibet (Gan et al., 2007;Zhang et al., 2004).

Support by viscous stresses
If the effective thickness of the elastic lid in central Tibet was greater than ∼20 km, then horizontality of the palaeoshorelines would give no information about the rheology of the middle and lower crust. If, however, the effective elastic thickness of the crust is only a few kilometres less than the critical thickness required to support the lake loads, then the crust below that lid must have viscosity of 4 10 19 -10 20 Pa s (Fig. 9), or 4 10 20 -10 21 Pa s if it is underlain by a weak substrate (Fig. 10). These bounds on viscosity, which apply on a time scale of 10 3 -10 4 yr, are consistent with the estimates of the viscosity of the crust in Tibet from post-seismic relaxation, on the time scale of years, which lie in the range of 10 18 -10 22 Pa s (Hilley et al., 2005(Hilley et al., , 2009Ryder et al., 2011;Wen et al., 2012;Yamasaki and Houseman, 2012). Our bounds, and the post-seismic estimates of viscosity, are consistent with laboratory experiments, and with post-seismic measurements of viscosity of the middle-to-lower crust in other parts of the world (Bürgmann and Dresen, 2008).
The palaeoshoreline observations do not allow us to place a strong constraint on the thickness of the layer of high viscosity in the crust. They permit a strong viscous layer to persist into the lowermost crust, as could be the case if the lower crust is anhydrous (e.g. Jackson et al., 2004), but they also permit the existence of a layer of relatively low viscosity in the lower crust if the viscosity of the middle crust is high enough (see calculations with a free-slip base to the layer, Fig. 10). Post-seismic estimates of the average viscosity of the whole of the ductile crust yield values of 10 18 -10 20 Pa s (Hilley et al., 2009), ∼10 19 Pa s (Ryder et al., 2011), and 2-5 Â 10 19 Pa s (Wen et al., 2012). Yamasaki and Houseman (2012), solving for depth-dependent viscosity of the crust beneath north-central Tibet, found that the viscosity at 10 km depth lies between 6 Â 10 20 Pa s and 10 22 Pa s, declining to ∼10 18 Pa s at the Moho.
It has been postulated that a channel exists within the Tibetan crust whose viscosity is so low that it accommodates most of the horizontal mass fluxes associated with active tectonics in Tibet, with little deformation of the crust above or below it (e.g. Beaumont et al., 2001;Clark et al., 2005;Clark and Royden, 2000;Royden et al., 1997;Morgan, 1985, 1987). The time constant for flow in such a channel is proportional to its viscosity and inversely proportional to the cube of its thickness (Clark and Royden, 2000;McKenzie et al., 2000); if the channel is ∼10-20 km thick, then its viscosity must be lower than about 10 16 Pa s for it to accommodate significant lateral flux of crust beneath Tibet (Clark and Royden, 2000). We cannot rule out the existence of such a channel, but if it does exist, it is overlain by a much more viscous layer that is some tens of kilometres thick. That high-viscosity layer, and the seismogenic upper crust to which it is coupled, deform in response to buoyancy forces associated with regionalscale (∼10 3 km) variations surface height (England and Molnar, 1997). The observation that active normal faulting within Tibet is confined to the regions with surface height above 4500 m (Elliott et al., 2010;Molnar et al., 1993) suggests that the buoyancy forces are resisted by local deformational stresses (∼300Àkm scale) rather than being dispersed by flow in a low-viscosity channel. Finally, even with the highest estimate of crustal viscosity discussed here (∼10 22 Pa s), the time scale for the crust to respond to regional topographic loads is short compared with the time scales of geological deformation (Bird, 1991;McKenzie et al., 2000, and Eq. (8), Appendix) so the low topographic relief of Tibet can be explained by flow of the entire crust, without the need to invoke large lateral fluxes of material confined within a low-viscosity channel.
where ρ s is the density of the visco-elastic substrate, here taken to be the same as that of the elastic lid. The Kelvin functions Ker, Ber, Kei and Bei, and their primed equivalents, are defined by Abramowitz and Stegun (1970, p. 379). The flexural parameter, λ, is Eq.
(2) was integrated numerically using Gauss-Legendre quadrature (Press et al., 1992), and the Kelvin functions were evaluated using the method of Chang and Jin (1996). To calculate the solutions involving a viscous half-space the layer thickness, L, was set to 1000 km.

A.1. Scaling relations
We may illuminate the results of the detailed calculations by simple scaling relations. We first consider the case in which the load is supported by an elastic lid floating upon an inviscid substrate. The maximum distortion of the land surface in this case is at the centre of the load, and is (Nakiboğlu and Lambeck, 1982) wð0Þ If the load was locally compensated (point-wise or Airy compensation) wð0Þ would be equal to h 0 ρ w =ρ c . This condition is not approached until the radius of the load is approximately twice the flexural parameter, λ, of the lid (Fig. 11a). (The depression of the surface below the level of Airy compensation for A=λ 42:7 is explained by the upward displacement of the substrate beneath the distal bulge of the elastic lid, Nakiboğlu and Lambeck, 1982.) For our purpose it is sufficient to note that if the radius of the load is smaller than the flexural parameter of the lid, the load will be supported primarily by elastic stresses in the lid, and the distortion of the surface will be small in comparison with that expected from local compensation of the load. If, however, A≳2λ, then the distortion of the surface will be comparable to that caused by local compensation of the load. The loads imposed by the lakes may also be supported by shear stresses associated with flow in the lower crust. The response time of a medium of viscosity η to a load of characteristic length A emplaced upon it is where ρ is the density of the medium, and the constant of proportionality depends on the shape of the load (e.g. Turcotte and Schubert, 2002). Fig. 11b shows that the response time of a viscous medium to a cylindrical load is well described by Eq. (7), with the constant of proportionality being ∼π. One other time scale is of interest to us, because we wish to relate our estimates of viscosity, based on kilo-year-scale measurements, to tectonic processes taking place over longer times. The time scale for the decay of a topographic load by pure-shear flow of a viscous layer is where l is the thickness of the layer (McKenzie et al., 2000, Eq. (A.21)).

A.2. Surface distortions due to distributed loads
For detailed calculations of distortions due to distributed loads, we discretise each of the loads due to the palaeolakes shown in Fig. 3 by subdividing each lake into 2.5-by-2.5 km squares. (This length scale is small in comparison with any horizontal length scale in the problem, so discretisation errors are negligible.) We represent the load of each square by a cylinder, concentric with the square, and of the same area, whose weight is equal to that of the water overlying the square. The displacements due to all these loads were calculated using Eq. (2) and summed to give the total displacement at each point of interest.
In each calculation, we adopt a simple loading history. The duration of loading is specified (time t 0 to t 1 , Eq. (2)), and the intensity of loading is taken to be constant within that time interval. Unloading is assumed to be instantaneous at the end of that time interval. The uncertainty introduced by assuming simple phases of loading and unloading is small in comparison with the uncertainties in the time intervals of the phases themselves.

Appendix B. Supplementary material
Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.epsl.2013.05.001. These data include Google Earth data of the palaeoshorelines described in this article. The following KMZ file contains Google Earth data of the palaeoshorelines described in this article. The file also contains the locations at which GPS measurements were made on palaeoshorelines around Zhari Namtso (Section 2.2); the precision of these observations has been degraded to conform with current Chinese law. The approximate locations of the photographs of Figure 2 are shown, as is our estimate of the position of a shoreline whose age was measured by Kong et al., 2011 (see