Three-dimensional rheological structure of the lithosphere in the Ordos block and its adjacent area

S U M M A R Y The 3-D rheological structure of the lithosphere in the Ordos block and its adjacent area (32◦– 43◦N, 102◦–116◦E) is obtained using iterative and averaging methods, based on the seismic velocity structure, thermal structure and composition of the lithosphere. The thermal structure is obtained from the newly published surface heat flow data sets. The brittle fracture mechanism is taken into account when the rheological structure is calculated. The strain rate obtained from GPS observation in the area under study is used to calculate the strength and the viscosity of the lithosphere. The results show that there are low-viscosity layers at the bottom of the upper crust, the lower crust and the bottom of the lithosphere both under the Ordos block and its adjacent regions, and the viscosity varies horizontally. The viscosity in grabens surrounding the Ordos block is one to two orders of magnitude lower than that in the Ordos block. The results also show that the rheological structure coincides with the major geological structure, below the major fault systems the viscosity is relatively low. The earthquakes with Ms ≥ 4.0 occurred along the large active faults. Most of the hypocentres of the earthquakes are in the upper and middle crust with high viscosity. The maximum depth of the earthquake is above the viscosity contour line of 1022 Pa s. Numerical studies show that the active geological structures must be taken into account when the 3-D rheological structure is used for geodynamics study.


INTROD U C T I O N
The Ordos block is located between the Tibet plateau and the North China plain. It is a stable geological block surrounded by seismically active deformed zones. There are no active tectonic events such as an uplift of the upper mantle, active Quaternary faults or large earthquakes within the Ordos block, but the surrounding deformed zones are very active and several large earthquakes with magnitude Ms ≥ 8.5 have occurred in the deformation zones (Deng et al. 1999). A variety of geological and geophysical researches have been carried out and a large amount of data has been accumulated in the Ordos block and the surrounding areas. These include seismic exploration for the crust structure , the new compilation of heat flow data (Hu et al. 2000(Hu et al. , 2001, GPS observation and analysis of GPS data (Yang et al. 2002;Li et al. 2003;Wang et al. 2003), seismic event location and accuracy , active faults and seismic activity (Han et al. 2003;Zhang et al. 2003) and deformation and movement of the geological blocks (Fan et al. 2003;Ning et al. 2003). Thus useful data for a rheological study are now available. * Now at: Department of Physics, University of Toronto, Toronto, Canada.
The rheological structure of the lithosphere is very important for the study of geodynamics. Considerable amount of work has been done on the rheological structure of the lithosphere and its application (e.g. Brace & Kohlstedt 1980;Kirby & Kronenberg 1987;Ranalli 1991;Meissner & Mooney 1998). The rheological structure has provided a useful tool for explanation of some geodynamic processes (e.g. Meissner & Kusznir 1987;Ranalli 2000;Meissner et al. 2002), although most of the study on the rheological structure of the lithosphere is determined for a 1-D model because for determining a 3-D rheological structure of the lithosphere more data sets and information are required. Zang et al. (2003) have obtained a preliminary model for the 3-D rheological structure of the lithosphere in North China using a constant strain rate and available data sets at that time. It shows that the preliminary model generally shows some correlation with the geological structure and reflects the main tectonic characteristics of North China. However, at the same time they found that the method for determining the rheological structure needs to be improved further and that more data and constraints are needed for the determination of the rheological structure.
In this paper we will use the newly available data to determine the 3-D rheological structure of the lithosphere in the Ordos block and its adjacent area with the strain rate constraint from the GPS results. Through the comparison between the velocity of the movement at GPS stations in this area related to the Eurasian plate from GPS observation and that obtained from the numerical modelling using this 3-D rheological structure, the application of rheological structure is discussed.

THE VELOC I T Y S T RU C T U R E A N D C O M P O S I T I O N O F T H E L I T H O S P H E R E
The study area for the Ordos block and its adjacent area in this paper is 32 • -43 • N, 102 • -116 • E (shown in Fig. 1a). Fig. 1(b) shows the distribution of active faults and focal mechanism solutions of earthquakes with Ms ≥ 5.0 around the Ordos block. The focal mechanisms solutions are for earthquakes that occurred after 1966 because dense station network was set up in China in 1966. Most of the faults are strike-slip faults. Some normal faults and thrust faults which formed before Quaternary period have movement with horizontal component larger than vertical component in Quaternary.
The faults surrounding the Ordos block form four grabens, that is, Yinchuan-Jielantai grabens, Shanxi grabens, Hetao grabens and Weihe grabens. According to Zhang et al. (2003) the North China is divided into several geological blocks. The Ordos block is separated from its adjacent geological blocks by these four grabens as shown in Fig. 1(c).
To determine the 3-D rheological structure of the lithosphere of this region, the 3-D velocity structure, composition and the 3-D temperature distribution of the lithosphere are needed. Zhu et al. (1997) set up a 3-D velocity database for continental China and its vicinity using a variety of deep seismic sounding data and the results of seismic tomography. The crustal velocity structure is based on the data from eight Global Geoscience Transects (GGT), more than 30 seismic exploration profiles and local crust structure studies. The average distance for the original observations varies from area to area, but it is less than 100 km in the study area. The distribution of these seismic exploration profiles can be found in Li & Mooney (1997). The velocity of the P wave at depths of 2, 5, 10, 20, . . . , 400 km is provided in a grid of 2 • × 2 • . We extrapolate this velocity distribution into a grid of 1 • × 1 • and obtain the velocity model for the studied area. According to seismic exploration studies (Li & Mooney 1997;Li et al. 2002), the crust was divided into three layers. The velocity distribution in the layers is shown in Table 1. The depth of the lithosphere is not determined clearly from the seismic study. Therefore, it will be taken from the thermal structure study (Section 3). There is an advantage for us to take the thermal thickness as the lithosphere thickness because the rheological properties of the mantle part of the lithosphere depend strongly on the temperature. The parameters of the 3-D layering model of the lithosphere used in this paper is shown in Table 1.
The composition of the lithosphere, especially of the crust, is very complicated. The composition of the lithosphere and their rheological properties were studied by some authors (e.g. Kirby & Kronenberg 1987;Ranalli & Murphy 1987;Ranalli 1997). Christensen & Mooney (1995) have presented models for seismic structure and composition of the continental crust, based on the data set of seismic refraction profiles and laboratory measurement of seismic velocity of the commonly observed lithologies of surface exposure. The composition and seismic properties of the mantle and crustal rocks exposed along several Geoscience Transects that pass through this area are also studied (Kern et al. 1996). According to the previous studies and the seismic structure in the Ordos block and its adjacent area, the composition of the lithosphere for the study area is given in the third column in Table 1. It should be pointed out that there are other possible choices, but in order to make the calculation possible we only select the minerals and rocks, which have creep parameters from experimental study.

THER M A L S T RU C T U R E O F T H E L I T H O S P H E R E
Temperature is the most important factor for determining the rheological structure of the lithosphere. The data and method of geothermal gradient estimation from observations are described briefly here.
In recent years, new heat flow data have been published in the third compilation of heat flow data in China (Hu et al. 2000(Hu et al. , 2001. The average surface heat flow distribution in a grid of 1 • × 1 • obtained from this new data set and used in this study is shown on Fig. 2(a). All temperature calculation from heat flow values are based on thermal conduction. We did not take into account any convection. Zang et al. (2002b) have discussed two kinds of heat production models for determining the thermal structure of the lithosphere in North China. One is that the crust is divided into three layers as shown in Table 1, and the heat production in each layer is constant. The other is that the crust is divided into two layers. In the upper layer the heat production decreased with depth exponentially until it attains the heat production in the lower part of the crust which is assumed to be a constant, taken as 0.4 µW m −3 (Chi & Yan 1998;Jaupart et al. 1998). They found that the thermal structure of the lithosphere in North China obtained using the second heat production model is more reasonable than that obtained using the model of three layers with constant heat production. So the two layer heat production model is used in this study. In the upper part of the crust, the heat production is A(z) = A 0 exp(−z/D), z is the depth in km, A 0 is the heat production rate near the surface at each grid, D is the characteristic thickness of the heat production layer (Pollack & Chapman 1977). A 0 and D are determined with the following equations, where A l is heat production rate of the lower part of the crust and taken as 0.4 µW m −3 (Jaupart et al. 1998), Z u is the thickness of the upper crust where A(Z u ) = A l , A is the average heat production rate at each grid. We assumed that A is proportional to surface heat flow in each grid, so A = A N /q 0 · q 0 = f · q 0 , A N = 1.25 µW m −3 and q 0 = 60 mW m −2 is the average heat production rate and average surface heat flow in North China, q 0 is the surface heat flow in each grid.
The 3-D temperature distribution of the lithosphere in the Ordos block and its adjacent area was calculated using iterative and running mean method (Zang et al. 2002b) based on the 3-D velocity structure (Zhu et al. 1997) and the heat production model mentioned above. Firstly 1-D temperature distribution with depth at the centre of each grid is calculated using iterative method. The equation and boundary are active in Quaternary and the horizontal component of motion is larger than the vertical one, otherwise the sense of offset is un-kown (Han et al. 2003;Zhang et al. 2003). The earthquakes occurred after 1966. The focal mechanism solutions are from Zhang et al. (1990) for earthquakes occurring during 1966-1976 and from CMT of Harvard for earthquakes from 1976 to 2004. The focal mechanism solutions are lower hemisphere projection. The black area is for compressional one. (c) Map of the study area. The names of the geological blocks are shown in bold italic. The dash-dotted line are boundary of the geological blocks. The Ordos block is separated by four grabens which consist mainly of faults in different geology age . Solid lines are the location of the cross-sections studied in this paper. S1   Note: a The velocity in the lithospheric mantle is not defined, and the lower thermal boundary is taken as the lower boundary of the lithosphere (see the text). b The creep parameters are for Weertman's formula σ = (ε/A) 1/n exp(E/n RT ) and from Ranalli (1997). condition are as follows, The depth of the lower boundary of the thermal lithosphere is determined from the constraint T = 1200 + 0.5z (Chapman 1986), here 0.5 is in • C km −1 and z is in km. Then the 3-D temperature distribution of the lithosphere is obtained by interpolating and running mean method.
The lower boundary of the thermal lithosphere is obtained by using the cubic spline smoothing method. It should be noted that the smoothing of the lower boundary of the lithosphere does not change the temperature distribution inside the lithosphere, so it does not change the rheological structure inside the lithosphere and only changes the shape and depth of the lower boundary of the lithosphere. Fig. 2(b) shows the depth contour of the lower boundary of the thermal lithosphere obtained in this study. Fig. 2(c) shows the temperature distribution along profile BB' in Fig. 1(c). It can be seen that the depth of the thermal lithosphere is from 70 to 200 km in the area under study. The gradient in the lithosphere thickness is large in some areas, although the lower boundary of the lithosphere is smoothed (Fig. 2b). The large gradient of the thickness of the lithosphere is due to the large gradient of the temperature, which is determined by surface heat flow and the crust structure. This will be discussed in Section 6.

RHEO L O G I C A L D E F O R M AT I O N M E C H A N I S M S
There are several deformation mechanisms in the lithosphere, including frictional sliding, brittle fracture and creep. At a certain temperature and pressure condition for any composition, one of them will become dominant. At present there are two problems in the study on the rheological deformation mechanisms: one is that different empirical relationships were used in the calculation of the strength of frictional sliding, and some of them lack the experiment basis; another is that the brittle fracture was ignored and only the frictional sliding and creep were taken into account. Zang et al. (2002a) pointed out that the reason for the neglect of brittle fracture mechanism in previous studies may be due to the lack of an appropriate criterion or formula for brittle fracture, because most results of brittle fracture experiments are based on the samples in the centimetre scale and their brittle fracture strengths are stronger 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 78  71  68  61  56  51  65  66  66  68  79  57  52  55   78  74  69  63  58  55  63  68  68  61  49  55  55  58   76  73  69  65  61  62  67  62  62  57  53  55  59  50   73  71  68  65  63  63  63  61  59  56  54  54  55  47   70  69  67  65  63  62  61  60  57  55  54  53  50  47   66  66  65  63  62  61  59  58  56  54  53  51  49 Zang et al. 2002aZang et al. , 2003. The calculation of the strength is for one-dimension model. The lithosphere structure, the heat flow density are the average of those for the Ordos block. The dashed line is the result considered the brittle fracture mechanism, while the solid line is the result without considering the brittle fracture mechanism. The values h1, h2 and h3 are the depths of the bottom of the upper crust, the middle crust and the lower crust, respectively. Brittle region (a) means that the dominant deformation mechanism is the frictional sliding in this region and brittle region (b) means in this region brittle fracture is dominant.
than that of frictional sliding. The experimental results show that the brittle fracture strength strongly depends on the size of the sample (Paterson 1978;Shimada 1993). According to the study of Shimada (1993) the strength of the brittle fracture of the rock sample in metre scale can be used to represent the strength of the brittle fracture  Shimada (1993); gabbro, basalt and peidotite are from Zang et al. (2002a). K, n, l are dimensionless. K, n of the granite is from Ohnaka (1973); grabbro from Wang et al. (1989); basalt from Zhang et al. (1985); peidotite from Mogi (1965) , T s and l of all the columns come from Wang et al. (1989).   Wang et al. (2003). The black arrows are the velocity vectors. The error ellipses represent 50 per cent confidence interval. (b) Distribution of principal axes of strain rate in the centre of grid of 1 • × 1 • obtained by Shen et al. (2003) from extrapolation and interpolation of the GPS observation. In this figure, T represents Tibet, NCP represents the North China Plain. The dash-dotted lines are the same as in Fig. 1(c). The scale of velocity and the strain rate are in the left up corner in each figure respectively. (c) Maximum shear strain rate distribution along profile BB . Solid line is the observed strain rate from Shen et al. (2003), the dotted line is the constant strain rate.
in the lithospheric condition. Fig. 3(a) shows the variation of the strength of frictional sliding and brittle fracture with depth. It can be seen that the brittle fracture strength for specimens in the centimetre scale is stronger than that of the frictional sliding and the brittle fracture can be ignored. The strength of the brittle  Fig. 1(c). The strength envelopes are calculated using the formulae (5), (6) and (7) with velocity structure mentioned in Section 2, and strain rate as shown in Fig. 4(b). The values h1, h2 and h3 in each subgraph represent the bottom of the upper crust, the middle crust, and the lower crust respectively. (b) 3-D viscosity structure of the lithosphere in the Ordos block and its adjacent area. The position of the five interfaces is along longitudes of 108.5 • E, 115.5 • E, latitudes of 32.5 • N, 39.5 • N and at depth of −60 km, respectively; The view azimuth is 140 • , and the angle of depression is 30 • . fracture for specimens in the metre scale is stronger than that of the frictional sliding at the shallow part of the crust, but becomes weaker than that of the frictional sliding in the deeper part of the crust (the transition depth depends on temperature and pressure), in this case the frictional sliding is dominant in shallow part of the crust and the brittle fracture mechanism will become dominant in the deeper part of the crust (below the transition depth). Fig. 3(b) shows the variation of average strength of the lithosphere for Ordos block with depth. It can be found that when the brittle fracture mechanism is taken into account the rheological strength in brittle region becomes obviously weak and the dominant deformation mechanism changes from frictional sliding to brittle fracture in the brittle region of the lower part of the crust and mantle. Therefore, most of the present rheological results not only increase the magnitude of strength, but also do not coincide with experimental results and real rheological mechanisms. Actually, the brittle fracture mechanism should be substituted for the frictional sliding mechanism in the deeper part of the lithosphere (Zang et al. 2002a). According to the previous discussions, the deformation mechanisms considered in this paper are as follows.  Zang et al. (2002aZang et al. ( , 2003 have chosen Byerlee's law (Byerlee 1978) as the criterion for the frictional sliding because this law was obtained by summarizing a large amount of experimental data and is used widely. The Byerlee's law (Byerlee 1978) is expressed as following, (3) Zang et al. (2002aZang et al. ( , 2003 gave explicit expressions of Byerlee's law for three kinds of stress states, for example, normal, thrust and strike-slip fault. We analyse the fault motion and the earthquake mechanisms in the lithosphere of the Ordos block and its adjacent area. The results show that most of active faults are strike-slip (horizontal movement is dominant). The axis of the minimum principal stress and maximum principal stress are almost horizontal, and the middle principal stress is almost vertical. This can be seen from Fig. 1(b). Therefore, the frictional sliding formula under strike-slip fault condition are used to calculate the strength of frictional sliding in this paper. The explicit expression of Byerlee's law for strike-slip fault, let

Frictional sliding mechanism and the stress state in the Ordos block and its adjacent area
where z is in km and σ 1 − σ 3 is in MPa.
For a quantitative estimation we take an intermediate value of k, that is, assuming k = 0.5, σ 2 = (σ 1 + σ 3 )/2 (Ranalli 1997 We hope to focus our attention on determining 3-D rheological structure. So the effect of pore pressure of fluid was not considered and assumed to be zero, because the situation of the fluid in the lithosphere is unknown and the pressure of fluid changes with position (e.g. Kohlstedt et al. 1995;Ranalli 2000).

Brittle fracture mechanism
The results obtained in laboratory and in situ show that the strength of the brittle fracture in the lithosphere can be studied with the experimental data of large-scale specimens (nearly in the metre scale) (Shimada 1993). Based on the work of Shimada (1993), Zang et al. (2002aZang et al. ( , 2003 extrapolated the results of small samples using the results of the large sample as a constraint and obtained an experiential formula of brittle fracture strength and the relative parameters. The brittle fracture strength is as following, where C 0 is the uniaxial compressive strength of rocks (metre scale),K , n, l, T s are the constants and shown in Table 2, σ c is the confining pressure. In this paper formula (6) and the parametres in Table 2 were used to calculate the brittle fracture strength of rocks.

Creep mechanism
Creep is the main deformation in the lithosphere. Weertman's power law is widely used to describe the creep in the lithosphere because the temperature and pressure are not very high in the lithosphere and the grain size is bigger (Weertman 1970;Meissner & Kusznir 1987;Ranalli 1997), whereε is the strain rate, T is the absolute temperature, R is the gas constant, E is the activation energy, A and n are the rheo- logical parameters independent of temperature and pressure. The columns 3-5 in Table 1 show the creep parameters of some representative rocks in the lithosphere. The creep parametres of the representative rocks are from small scale and generally show a large scatter.

RHEO L O G I C A L S T RU C T U R E I N T H E O R D O S B L O C K A N D I T S A D J A C E N T A R E A
The lithosphere seems to behave as a non-linear rheological material, but it can be described by an effective viscosity (Kirby 1985;Meissner & Kusznir 1987;Meissner & Mooney 1998), where σ is the rheological strength calculated andε is the strain rate. It can be seen that the effective viscosity is the function of strain rate (hereafter it will be called simply as viscosity and written as η). With few exceptions, the rheological strength or viscosity of the lithosphere was usually determined using a constant strain rate, usually 10 −14 s −1 or 10 −15 s −1 (e.g. Ranalli 2000;Meissner et al. 2002;Zang et al. 2003), but it is not useful for describing the real rheological behavior. In this study, the strain rate obtained from the GPS data (Yang et al. 2002;Shen et al. 2003) is used to determine the rheological strength and the viscosity. The GPS velocity field for the Chinese continent is obtained by analysing GPS data from the Crust Motion Observation Network of China (CMONOC) from 1998 to 2001, particularly the data from the regional networks (nearly 1000 sites) of CMONOC observed in 1999 and 2001. Fig. 4(a) shows the GPS velocity field at GPS stations obtained by Wang et al. (2003), and Fig. 4(b) shows the strain rate distribution in the centre of grid of 1 • × 1 • in the study area. The strain rates and rotation rates in each grid cell were derived from the GPS velocity field by Shen et al. (2003) using a localized regression fitting method (Shen et al. 1996). At any given spot and its vicinity the velocity field is assumed to be a linear function of space, and represented by a combination of the on spot velocity, strain rate and rotation rate. To fit the model the velocity data at each site are reweighted according to the distance between site and the spot where the strain rates are to be evaluated. The weighting function is in the form of exp(−L 2 /φ 2 ), where L is the distance mentioned above and φ is a 'smoothing' factor. To optimally model the φ is chosen based on in situ strength of the data (Shen et al. 1996). Fig. 4(c) shows the maximum shear strain rate distribution along profile BB' (Fig. 1). The maximum value of the strain rate is about 5.2 × 10 −15 s −1 , the minimum value is about 2.2 × 10 −16 s −1 . From Figs 4(b) and (c) it can be seen that the strain rates obtained from GPS data vary from place to place.
Using the formulae of (5), (6), (7), with the data of the velocity structure, composition, temperature, pressure in the lithosphere and strain rate from GPS observation (Yang et al. 2002;Shen et al. 2003), the rheological strength of the Ordos block and its adjacent area was estimated. In the calculation of the strength, the iterative method was used to calculate the strength distribution with depth in each centre point of each 1 • × 1 • grid. For example, Fig. 5(a) shows the strength envelopes at two points along an E-W transect through grabens and the Ordos block. Finally the 3-D rheological strength structure of the lithosphere in the Ordos block and its adjacent area was obtained using the running mean method. Then the corresponding viscosity structure can also be obtained using the formula of (8). Fig. 5(b) shows the 3-D viscosity structure of the lithosphere in the Ordos block and its adjacent area. It can be seen that the rheological structure is layered in the vertical direction, but the thickness and the depth of each layer varies horizontally.
For detail illustration, six cross-sections of the viscosity structure of the lithosphere (AA -FF , as shown in Fig. 1) are shown on Fig. 6. Fig. 7 shows the distribution of the viscosity on the horizontal slices at depth of 5, 15, 25, 35 and 45 km, respectively. The thermal thickness of the lithosphere is also shown on Fig. 6. It can be seen from Figs 6 and 7 that there is a low-viscosity layer near the surface and three obvious low-viscosity layers at the bottom of the upper crust, the lower crust and the bottom of the lithosphere below both the Ordos block and the adjacent grabens.
It also can be seen from Figs 6 and 7 that the viscosity varies horizontally at different depths. From Fig. 6 it can be found that the Ordos block has higher viscosities by about one to two orders of magnitude than that of the surrounding grabens at some depths. Furthermore, the viscosity in the interior of the block is not uniform. The western block has much higher viscosity than the eastern block.

The relationship between the rheological structure and geological structure and earthquake distribution
In this section, we will study the relationship between the rheological structure and the geological structure and earthquake distribution. The positions of the major geological units are indicated at the top of each cross-sections in Fig. 6. It can be seen from Fig. 1 that the cross-sections AA , BB and CC' pass through the Yinchuan-Jielantai grabens, the Ordos block, the Shanxi grabens and the North China plain from the West to the East, except cross-sections AA passes through the Hetao grabens in the East instead of passing through Shanxi grabens. The cross-sections DD , EE and FF pass through the Weihe grabens, the Ordos block and Hetao grabens from the South to the North, except the cross-sections FF passes the Shanxi grabens in South instead of the Weihe grabens. The Yinchuan-Jielantai grabens, Shanxi grabens, Hetao grabens and Weihe grabens all consist of several big fault systems and are very active Han et al. 2003). This can be seen from Fig. 1. It can be found from Fig. 6 that the lithosphere is thinner and the viscosity is lower relatively in the Yinchuan-Jielantai grabens, Shanxi grabens and Weihe grabens than in the Ordos block. The viscosity is also relatively low in the Hetao grabens, but it is not very obvious. This shows that the lithosphere in the grabens is relatively weak and seismically active, the lithosphere in the Ordos block is strong and stable. Inside the Ordos block the viscosity is lower in the eastern part of the Ordos block than in the western part. The crosssections AA and FF are close to the boundary, so the viscosity distributions are not as simple as that in other cross-sections. In general, the rheological structure has a close relation with the major geological structure in the study area. Fig. 8(a) shows the distribution of the earthquakes with Ms ≥ 4.0 in the Ordos block and its surrounding area. The earthquakes used occurred from 1985 to 2004 because the local and national seismic networks were set up and the hypocentre are determined more accurately. The hypocentre are determined by the Sub-Center of China National Digital Seismic Network using local and national seismic networks. The errors in depth are less than 5 km Yang et al. 2003). It can be found from Fig. 8(a) that the earthquakes mainly occurred in four grabens. According to historical recording, more than five earthquakes with magnitude ≥ 8.0 have occurred in these grabens since AD 1000, and all the large earthquakes are shallow with large surface ruptured zone (Han et al. 2003). Some earthquakes occurred in the North China plain, Alashan and Tibet. A few occurred in the Ordos block. Fig. 8(b) shows the distribution of earthquakes with Ms ≥ 4.0 along the cross section BB . The earthquakes occurred within 50 km of both sides of the cross-sections BB , but the earthquakes that occurred in Alashan and Tibet are excluded because the seismic observation network in Alashan and Tibet is less dense, so the uncertainty of the depth is large. It can be seen in Fig. 8b that most of the earthquakes are located in the upper and middle crust, but a few earthquakes occurred in the lower crust. Nearly all of the earthquakes are located above the viscosity contour line of 10 22 Pa s. Fig. 8c shows the maximum depths of the earthquakes with Ms ≥ 4.0 occurred in the grid of S1, S3 and S4 from 1985 to 2004. It can be seen in Fig. 8(c) that the earthquakes occurred in the upper and middle crust in the grids S1, S3 and S4. This shows that the upper and middle crust in the grabens are in the brittle regime, consistent with our model.

The influence of active faults on the rheological structure
We have actually made an assumption that the effect of the faults are the same in any place and direction when we calculated the rheological strength of the crust using frictional sliding. Is this assumption valid for the geological block with big active fault system, such as the Ordos block and its adjacent area because the big faults are basically orientated in certain directions (Fig. 1b)? In this section we will test this assumption using numerical modelling. Two 3-D numerical experiments are carried out using a computer program of the finite The strength envelopes and the maximum depths of the earthquake at three points of S1, S3 and S4 as shown in Fig. 1c. The strength envelopes are calculated using the formulae (5), (6) and (7) with velocity structure mentioned in Section 2, and strain rate as shown in Fig. 3(b). The dots with error bar in depth are the maximum depths of the earthquakes from 1985 to 2004 at the three points. The material in the two numerical models is also the same and is the Maxwell body, but the viscosity for the Maxwell body in the two models is different. In Model 1, the viscosity is obtained from the 3-D rheological structure obtained in this study (Section 5). In Model 2, the viscosity in the grabens is reduced to about 20 per cent of its original one used in Model 1, but the viscosity in another area is the same as in Model 1. In the first numerical experiment the velocities at site of GPS stations are modelled using Model 1 and compared with the observed ones as shown on Fig. 9(a). In the second numerical experiment the velocities at site of GPS stations are modelled using Model 2 and also compared with the observed ones as shown on Fig. 9(b). It should be pointed out that in the numerical experiments only the standard GPS observation stations are used because the observations in those stations are stable and the uncertainty of the results is small. It can be seen from Fig. 9(a) that the velocities obtained from the modelling are basically consistent with the observed ones. This shows that the 3-D rheological structure obtained in this paper is reasonably good, but it can also be found that the differences between velocities obtained in the numerical experiment and those observed are very obvious on some GPS sites. When the big fault systems are taken into account (the viscosity in the grabens is reduced), the differences are reduced, which can be seen in Fig. 9(b). It indicates that the active fault oriented in certain direction has quite major influence on the rheological structure. It should be pointed out that we only considered the effect of the fault here, but there are other factors can affect the rheological structure, such as the topographic or subsurface loads or bending stresses. These factors must be taken into account when the rheological structure of the lithosphere is used for geodynamic study.

Comparison of the rheological structures obtained from the constant strain rate and the strain rate from GPS observation
The 3-D rheological structure of the lithosphere is obtained based on the effective viscosity. The definition of the effective viscosity is analogous to that for Newtonian materials (linear viscosity), but it should be noticed that the effective viscosity is not only a function of temperature, pressure and material parameters, but it is also a function of stress (or strain rate). So the strain rate (or stress) is a very important factor for viscosity determination. Conventionally, the effective viscosity (hereafter, call it as viscosity again) is determined using constant strain rate, such as 10 −14 or 10 −15 s −1 . The viscosity or the rheological strength of lithosphere obtained under this assumption is useful for comparative purposes, but it is not reasonable to use them for numerical modelling or other geodynamical studies. In this study the viscosity is determined as described in Section 5 using the strain rate obtained from GPS observation. The observed strain rate from GPS should be more realistic than an assumed constant value. Figs 10(a)-(d) show the rheological strength and viscosity obtained both using the constant strain rate equal to 10 −15 s −1 and observed strain rate on the cross-sections BB . It can be seen that the general pattern of strength and viscosity distribution is not much different because the variation of observed strain rate across the region is small and the main characteristics of these distribution is mainly controlled by the distribution of temperature. The changes of the rheological strength and viscosity vary with place and depth. For quantitative estimation, Fig. 10(e) and Fig. 10(f) show the viscosity distribution at the point S2 inside the Ordos block and S3 inside the Shanxi graben, respectively (Fig. 1c). The maximum change of viscosity is about 2.2 × 10 23 Pa s at S2; while at S3 the viscosity is hardly changed because the observed strain rate obtained from GPS is nearly the same magnitude as the constant strain rate taken as 1.0 × 10 −15 s −1 . If the constant strain rate is taken as 1.0 × 10 −14 s −1 , the viscosity at S3 will have a maximum change of 2.0 × 10 23 Pa s. When the observed strain rate is used as a constraint on the determination of viscosity of the lithosphere, the physical meaning of viscosity of the lithosphere becomes more clear and useful. It can be used for qualitative discussion and for quantitative study. Following the development of geodetic methods, it is possible to get more and more accurate strain rate data for better constraint on viscosity distribution of the lithosphere, although how to extrapolate the surface strain rate to depth needs further study.

Effect of the uncertainty of the heat flow data on the lithospheric thickness and rheological structure
In order to analyse how heat flow uncertainties effect the lithospheric thickness and final model, we change the observed heat flow value in some grids and recalculate the lithospheric thickness. It can be seen in Figs 2(a) and (b) that the observed heat flow in Grid 1 (35 • -36 • N, 103 • -104 • E), Grid 2 (32 • -33 • N, 107 • -108 • E) and Grid 3 (37 • -38 • N, 107 • -108 • E) is lower than in their adjacent grids, respectively. The thickness of the lithosphere is thicker in the Grid 2 (130 km) and Grid 3 (150 km) than in their surrounding grids. There are large gradients near the Grid 1. So we chose these three grids as samples to test the effect of uncertainty of the heat flow data on the rheological structure. Firstly, the heat flow values are increased by 10 per cent in the three grids as shown in Fig. 11(a). Then the thickness of the lithosphere is calculated in the three grids simply for 1-D model, using formula (2) and its boundary condition in Section 3. The result shows that the lithospheric thickness without smoothing changed by 19.2, 22.2 and 18.9 per cent in the three grids, respectively.
Secondly, the heat flow values in three grids are increased by 10 per cent and the observed heat flow values in other grids are not changed. Then the 3-D thermal structure of the lithosphere is recalculated using the procedure described in Section 3. The heat flow distribution after extrapolation is shown on Fig. 11(a). The depth contour of the lower boundary of the recalculated thermal lithosphere is shown in Fig. 11(b). Through the comparison of Fig. 2(a) with Fig. 11(a) and Fig. 2(b) with Fig. 11(b), respectively, it can be found that the extrapolated heat flow values in the adjacent grids of the Grid 1, Grid 2 and Grid 3 are also changed, but the variation of the heat flow in the adjacent grids is less than 10 per cent. The lithospheric thickness also is changed. The lithospheric thickness in Grid 1, 2 and 3 is 141 km, 129 km and 140 km, and its variation with the original one is about 9.4, 1.1 and 6.5 per cent, respectively. The large gradients in the boundary of Tibet block (bottom left hand corner in Fig. 11b) still exist although they become smaller. There are two possibilities for this large gradient. One is due to the lack of observed heat flow data in this region and the extrapolated data may have large uncertainty. Another is that the large gradient may reflect the geological structure because the gradients of the thermal thickness of the lithosphere parallel those of the crustal thickness in this region. Both of the two possibilities need further study using more heat flow data. Fig. 11(c) shows the temperature distribution along profile CC obtained from the heat flow distribution shown in Fig. 2(a) and shown in Fig. 11(a), respectively. The profile CC passes through the Grid 1 where the thermal thickness of the lithosphere changed with maximum value of 9.4 per cent. It can be found that the characters of the two temperature distributions in Fig. 11(c) are basically the same in most part of the profile, except the region from 102 • E to 105 • E. The maximum difference of the temperature occurred below the depth of 100 km near 103 • E. It can cause some changes in viscosity structure, but the changes are not remarkable. Fig. 11(d) and Fig. 11(e) show the viscosity distribution along profile CC calculated under these two different temperature structures. It can be seen that the difference of the two viscosity distributions is very small.

Conclusions
(1) The lithosphere in the study area is considered to be rheologically layered vertically. Constraints on layers are provided by the seismic structure and composition of the lithosphere, so the   boundaries of the rheological structure in the crust are basically similar to seismic boundaries. The horizontal variation of the viscosity is mainly controlled by the variation of temperature (the temperature variation is connected with tectonics). The depth of the lower boundary of the thermal lithosphere is constrained by the mantle adiabat line (Section 3), so the horizontal variation of the viscosity in the subcrust is connected with the variation of the lower boundary of the thermal lithosphere.
(2) The magnitude of the viscosity varies from 10 19 to 10 24 Pa s; The western side of the Ordos block has one to two orders of magnitude higher viscosity in the lithosphere than that beneath adjacent grabens.
(3) The results of our numerical simulations suggest that there are some connections between the rheological structure and major faults. Below the major faults the viscosity is relatively low. The earthquakes with Ms ≥ 4.0 occurred from 1985 to 2004 along large active fault systems, and most of the hypocentre of the earthquake are located in the high-viscosity layers in the upper and middle crust.
The maximum depth of the earthquake is above the viscosity contour line of 10 22 Pa s.
(4) The major faults have strong influence on the rheology of the lithosphere. This must be taken into account when the rheological structure is used in geodynamic studies and further rheological studies.

A C K N O W L E D G M E N T S
We are grateful to Prof. Zhu J.S for kindly providing the velocity database and Prof. Shen Z. K., Prof. Yang G. H, Prof. Li Y.X and Wang M. for their bountiful supply of GPS data. We are also grateful to Sub-Center of China National Digital Seismic Network for the earthquake data. We benefitted from constructive rigorous comments and suggestions by two anonymous referees. This work was supported by the Project of Continental Earthquake and Forecast (Grant No. G95-13-04-06) and the National Natural Science Foundation of China (Grant No. 40174023).