Three‐Dimensional Thermal Structure of East Asian Continental Lithosphere

The thermal structure of East Asian continental lithosphere is critical for understanding its diffuse and variable intracontinental deformation in the Cenozoic. Here, we present a three‐dimensional model of the lithospheric thermal structure of East Asia, using the latest thermal conductivity and radiogenic heat production measurements in mainland China. The results show great lateral heterogeneity of lithospheric temperature. The crust of orogenic belts such as the Tibetan Plateau, the Mongolia Plateau, the Tian Shan orogen, and Northeast China, is generally characterized by high temperatures. Other regions of high crustal temperatures include the southeast coast of China and Lake Baikal. The heat sources for the crust also vary. In the Tibetan Plateau, high crustal temperatures result mainly from radiogenic heat production within the thick crust, indicating the influence of continental collision. However, high heat flow from the upper mantle is the main cause of high surface heat flow in the eastern part of East Asia, where the relatively thin lithosphere and hot asthenosphere is related to the subduction of the Pacific Plate. This subduction is also responsible for Cenozoic continental volcanoes in East Asia.

Much of the East Asian continent is formed by the ancient cratonic blocks, including the NCC, the South China Block, and the Tarim Basin. Steady-state heat conduction may be acceptable for calculating lithospheric temperature in these regions. However, steady-state heat conduction may not hold in active tectonic regions. Transient heat perturbation could result from lithosphere thinning. In this case, steady-state heat conduction may still hold for the crust, because it would take tens of millions of years for the perturbing heat to transfer from the mantle to the crust due to the inefficiency of heat conduction. However, in active orogenic belts such the Tibetan Plateau, perturbing heat may result from both lithosphere-asthenosphere interaction, and from crustal thickening and deformation (Craig et al., 2012;Liu & Shen, 1998;McKenzie & Priestley, 2016). Thus, steady-state heat conductive calculations in these areas may only provide some reference temperature.
In addition to the three-dimensional model that considers the variable lithospheric structure, it is also necessary to consider the variation of thermal parameters. For instance, because the lack of measurements, radiogenic heat production was sometimes obtained through the empirical relationship between seismic wave velocity and heat production (Rybach, 1976;Rybach & Buntebarth, 1982). In recent years, the surface heat flow and heat conductivity data in East Asia, especially in mainland China, have grown significantly (Jiang et al., 2019), and the SinoProbe project  has obtained the abundance data of near-surface radioactive elements covering entire mainland China (Wang et al., 2016;Wang, 2015). Here, we present the latest near-surface radiogenic heat production measurements in mainland China. Although the ways in which thermal parameters measured near the surface vary at greater depths in the crust and mantle lithosphere is still uncertain, these newly available measurements provide some first-order constraints on the lateral heterogeneity for models of three-dimensional temperature structure of the lithosphere in the East Asian continent.
In this study, we constructed a three-dimensional finite element model based on Crust 1.0 model (https://igppweb. ucsd.edu/∼gabi/crust1.html) for the East Asian continental lithosphere. Constrained by East Asia's annual average surface temperatures and the upper mantle temperature or thermal lithospheric thickness obtained from seismic velocity inversion (An & Shi, 2006), the lithospheric thermal structure is solved using the three-dimensional steady-state heat conduction equation. The results show large spatial variations of lithospheric temperatures across different tectonic units that are helpful for understanding the diffuse Cenozoic tectonics and intracontinental volcanism in East Asia.

Equations and Numerical Model
In the spherical coordinate system, the three-dimensional steady-state heat conduction equation can be written as: where T is the temperature (K); r, θ, and φ denote radius, colatitude, and longitude, respectively; k is thermal conductivity (W ⋅ m −1 ⋅ K −1 ), and Q represents the heat production (W ⋅ m −3 ). The heat flow q r (W ⋅ m −2 ) in the r (or vertical) direction is The governing Equation 1 is solved in a finite element numerical model, which is established in the spherical coordinate system and considers the influence of the curvature of the Earth (Figure 2). The area covered by the numerical model is 70-136°E and 20-52°N. The vertical positions of the nodes on the surface are obtained by For cases A and B, the bottom depth is 100 km. For case C, the bottom depth is set according to the thermal lithospheric thickness (An & Shi, 2006). interpolation based on the topographic data. The numerical models are divided into the upper crust, middle crust, lower crust and lithosphere mantle according to the Crust1.0 model (https://igppweb.ucsd.edu/∼gabi/crust1. html). Because the sediments are relatively thin (<1 km) in most areas or missing ( Figure S1 in Supporting Information S1), we do not consider it as a separate layer but part of upper crust. However, in some regions, such as the Tarim Basin, the Sichuan Basin, and the Ordos, the sedimentary thicknesses are over 2 km, and potentially may have impacted the results because the observed thermal parameters are from the near-surface, which cannot represent the whole crust. The mesh is divided according to the actual topography and the thickness of each layer. In this way, the heterogeneity of crustal thickness is fully considered in the model ( Figure S1b in Supporting Information S1). The finite element model uses triangular prism elements for meshing. The total number of elements and nodes of the numerical model are 1,755,312 and 924,050, respectively.

Model Parameters
Solving the three-dimensional steady-state heat conduction equation requires data for the thermal conductivity and heat production parameters (Equation 1). In recent years, the heat flow data of the East Asian continent has greatly increased. Thermal conductivity measurements were also carried out at most of the heat flow observation points (Figures 3a and 3c). In mainland China, 1,136 heat flow data points are now available (Jiang et al., 2019). The data have been categorized into four classes (Categories A-D) depending on the quality and reliability. Categories A and B represent reliable heat flow data for regional geothermal studies. As a necessary complement, Category C data has greater uncertainty because of the unavailability of detailed reference information. Category D data are strongly influenced by deep geothermal fluids or near-surface factors. Date from category D only serve as the essential references for the local geothermal background and geothermal resource investigation (Hu et al., 2000;Jiang et al., 2019). Thus, in this study, we eliminated all the category D data (48 points) in the mainland China. In the study region, the heat flow and thermal conductivity data are supplemented by data from the International Heat Commission database (http://www.ihfc-iugg.org/products/global-heat-flow-database); these data have no detailed quality grading.
As part of the SinoProbe project (http://www.sinoprobe.org/), the China Geochemical Baselines (CGB) was launched in 2008, and sampling was completed in 2012. The purpose of CGB is to document the abundance and spatial distribution of chemical elements covering all of China (Wang, 2015;Wang et al., 2016). The CGB with approximately 1,600 reference grids (1° × 1°) covers the whole China. Two sampling sites were allocated to each CGB grid cell. At each site, two samples were taken: one is from a depth of 0-25 cm and a second, deeper sample from a depth greater than 100 cm. A total of 6,617 samples from 3,382 sites have been collected. The baselines data and maps of 76 elements have correlation with geology, mineral resources, climate and human activities. More detail of methodology and data can be found in the references (Wang, 2015;Wang et al., 2016). We used 1,405 grids data in China mainland because some of reference grids have not been completed ( Figure 3e). The grid data (113 grids) with same sampling methodology have been obtained in south Mongolia through international cooperation between China and Mongolia. This data has produced systematic abundance measurements of radioactive elements (uranium, thorium, and potassium) over mainland China and southern Mongolia. Here, for the first time, we present 1,518 (1° × 1° cell) near-surface radiogenic heat production data (Table S1, Figure 3e) based on the abundance measurements of radioactive elements.
The near-surface measurements, including heat flow, thermal conductivity, and radiogenic heat production, provide good constraints for the thermal structure calculation. However, these data are unevenly distributed. Therefore, we constructed the near-surface heat flow, thermal conductivity, and heat production of the East Asian continent through Kriging interpolation. These near-surface observations are interpolated into the surface nodes of the numerical model (Figures 3b, 3d, and 3f). However, the measured thermal parameters are only distributed on the near-surface. Additional constraints are needed for their depth-dependent variations.
The thermal conductivity is estimated from the following equations (Chapman, 1986;Doin & Fleitout, 1996;Jaupart & Mareschal, 1999;Wang, 2001):  , thermal conductivity (c and d), and radiogenic heat production (e and f). The heat flow and thermal conductivity data are from the International Heat Commission database (http://www.ihfc-iugg.org/products/globalheat-flow-database) and data from mainland China (Jiang et al., 2019). The radiogenic heat production data (1,518 1° × 1° cells) is collected from this study based on the radioactive element abundance data. The density used to calculate the heat production in this figure is 2,600 kg/m 3 . where k 0 is the measured thermal conductivity of the surface element in the numerical model. T is the temperature (K) at depth z. For the upper, middle, and lower crusts, b is 1.0 × 10 −3 K −1 , 5.0 × 10 −4 K −1 , and 1.5 × 10 −4 K −1 , respectively; c is 1.5 × 10 −3 km −1 .
The radiogenic heat production in the crust and mantle is determined by the abundance of radioactive isotopes of uranium ( 238 U; 235 U), thorium ( 232 Th), and potassium ( 40 K), which can be described by (Rybach, 1976): where Q 0 is the heat production in near-surface rocks ( W ⋅ m −3 ), ρ is the density (g ⋅ cm −3 ), which is obtained based on the Crust 1.0. C U 、C Th 、C K are the abundance of uranium (μg/g), thorium (μg/g) and Potassium (wt%), respectively. Generally, radiogenic heat production is assumed to decrease exponentially with depth (Lachenbruch, 1970;Lachenbruch & Bunker, 1971), where Q is the heat production at depth z; D (unit: km) is the characteristic depth over which most of the radiogenic heat production element are distributed. Previous studies have shown that the distribution of the three radioactive elements in the crust affects the relationship between heat production and depth (Jaupart et al., 1981). Some research also shows that D decreases with depth (Čermák & Rybach, 1989). Deep drilling in the Baltic Shield (Kremenetsky et al., 1989), the results of Deep Drilling Program of Germany (KTB) (Pribnow & Winter 1997) and scientific drilling in mainland China (He et al., 2009;Qiu, 2002) found no obvious decrease of the heat production with depth in the upper crust, and the heat production strongly depends on lithology. In some ultrahigh-pressure metamorphic zones, radiogenic heat production even increases with depth in the upper crust (He et al., 2009;He et al., 2006). Recent work suggests that there is not a simple relationship between radiogenic heat production and metamorphic grade (Hasterok & Webb, 2017;Hasterok et al., 2018;Perry et al., 2006) and significant radiogenic heat production may remain in the lower crust in cratonic regions. Therefore, Equation 5 may underestimate the heat production of radioactive elements in the lower crust (Ketcham, 1996). Here, we adopted an alternative model for inversion of upper and middle crustal heat production in cases B and C (Čermák & Rybach, 1989;Turcotte & Schubert, 2002): For the upper and middle crust, in the vertical direction, the heat production remains unchanged in the depth less than D 0 . Below the depth of D 0 , the heat production decreases exponentially. In the horizontal direction, D 0 and D 1 can vary for different sites. In the calculation, D 0 and D 1 are given within a certain range. The range of D 0 is the change from the thickness of the sedimentary layer to the thickness of the upper crust based on the Crust 1.0 model. According to previous studies, the average value of D 1 in the continent is ∼10 km, which ranges from 5 to 15 km (Artemieva & Mooney, 2001;Jaupart et al., 2016). In this study, D 1 varies in 0-15 km. D 0 and D 1 values are obtained through fitting the surface heat flow. However, for the lower crust, constant radiogenic heat production, which represents the residuum of melting where radiogenic heat production primarily resides with potassium remaining in feldspar, may be more appropriate (Artemieva & Mooney, 2001;Hasterok & Chapman, 2011). For the lower crust in cases B and C, radiogenic heat production is kept constant at the level at the bottom of the middle crust according to the Equation 6.
As described above, it is likely that the true distribution of radiogenic heat production in the crust does not resemble the model adopted in many areas. If there are layers with relatively high radiogenic heat production at depth within the crust, they could have a significant effect on the crustal temperatures at a particular depth. It is not possible to determine the true distribution of radiogenic heat production with depth at present, but the model used here is based on geochemical principles and observations and is flexible enough to capture broadscale lateral variations in the importance of radiogenic heat production in controlling crustal geotherms. It is intended to highlight these variations and to provide a reference model for future more detailed studies to be compared against.

Boundary Conditions
Solving steady-state heat conduction equations required two boundary conditions. We constructed three numerical models (cases A-C) with different boundary conditions. In case A, the observed heat flow ( Figure 2b) and temperature ( Figure S2a in Supporting Information S1) are imposed on the upper surface as the thermal boundary conditions. The temperature for the upper surface (Ts) of the model domain is the annual average ground temperature obtained from 195 meteorological stations ( Figure S2a in Supporting Information S1; Sun et al., 2013). In case B, the upper surface temperature boundary condition is the same as in case A. The other boundary condition is the temperature at 100 km depth estimated from S-wave velocity (An & Shi, 2006; Figure S2b in Supporting Information S1). Case C is similar to case B, except that the bottom of the model domain is set at the depth of the thermal lithosphere estimated from seismic data (An & Shi, 2006; Figure S3a in Supporting Information S1), with a fixed isentropic temperature (1280°C) (McKenzie & Bickle, 1988). There are also some lithospheric thickness models based on temperature (Priestley et al., 2018; Figure S3b in Supporting Information S1) or seismic velocity (Lith1.0) (Pasyanos et al., 2014; Figure S3c in Supporting Information S1). The spatial pattern of the lithosphere thickness in these models is similar ( Figure S3 in Supporting Information S1). In this study, we selected the thermal thickness of lithosphere from An and Shi (2006). Because the lithospheric thickness from this model is calculated based on the temperature and its regional resolution (1° × 1°) is higher than the global model ( The crustal structure for all cases is the same.
For cases B and C, to make full use of existing observations of near-surface heat flow, the fitting error between the calculated and the observed surface heat flow is minimized by adjusting the exponential attenuation parameters D 0 and D 1 . The values of D 0 and D 1 corresponding to each surface element of the model are obtained in the inversion calculation. A flow chart for calculating the thermal structure of the lithosphere is shown in Figure S4 of Supporting Information S1.

Overview of the Three Thermal Models
Case A uses the surface observables, heat flow ( Figure 3b) and temperature ( Figure S2 in Supporting Information S1), as the thermal boundary conditions. The lateral variations of thermal parameters (Figures 3d and 3f) are considered in this case. Theoretically, this is a preferred choice, and has been used by others (Artemieva & Mooney, 2001;Hasterok & Chapman, 2011). However, we found that the surface observables are not reliable constraints for the 3-D thermal structure of the lithosphere. Figures 4a-4b show the results when radiogenic heat production was not included in the calculation. This means that the surface heat flow is the same as the heat flow throughout the Moho or the lithosphere, which would require a higher temperature in the lower lithosphere than in cases where contribution from radioactive heating to the surface heat flow is considered. Thus, temperatures in case A may be regarded as the upper limit values. The panels show the temperature at 20 km depth and on the Moho surface. The calculated temperature in some areas is unreasonably high. For example, 800°C at 20 km depth in some regions, which would cause significant anatexis. The white regions are where the temperature surpasses 1300°C, generally considered the temperature at the base the lithosphere (Figure 4b). Temperature over 1300°C at the Moho surface implies that the thermal lithospheric thickness is smaller than the crustal thickness ( Figure S1b in Supporting Information S1). If the radiogenic heat production is included in the calculation, the temperature would be decreased due to the reducing heat flow with depth. The value of D = 20 km, which presents characteristic depth of the radiogenic heat production, is assumed for the whole model domain (Equation 5). This value is larger than the commonly used value of D = 10 km for the Precambrian cratons (Artemieva & Mooney, 2001). However, after adding the radiogenic heat production, the temperature in some regions is still higher than 1300°C (Figures 4c-4d). This may imply that the radiogenic heat production is underestimated for part of East Asia using Equation 5 or a higher thermal conductivity.
The areas where the Moho temperature is larger than 1300°C are mainly in active Cenozoic tectonic zones, such as the Tibetan Plateau, the Tian Shan, and Cenozoic volcanoes. They are generally consistent with previous study (Artemieva, 2006). There are numerous factors influencing temperature estimation in these regions. First, in these areas, there are not enough surface heat flow and conductivity data to provide reliable constraints (Figures 3a  and 3c). Some low-quality surface heat flow observation (Category C) in these areas may be elevated by underground fluids. For instance, the Tibetan Plateau appears as a high heat flow area after interpolation (Figure 3b). This would lead to high estimation of the lithospheric temperature in case A. Secondly, the steady state assumption may be not suitable for these tectonically active areas. In these areas, the lithosphere may not have reached steady state. The selection of thermal conductivity is not well constrained and may also impact the results. In order to avoid the over influence of the near surface measurements, it may be more reasonable to combine surface measurements with some deep temperature constraints.
In cases B and C, the upper mantle temperature obtained from seismic wave velocity anomalies is used to constrain deep temperature. The difference of case B and case C is the bottom depth of the model domain. In case B, a constant depth of 100 km is used. The temperature at the 100 km depth is from the results of An and Shi (2006). Whereas in case C, the thermal lithospheric thickness is imposed on the bottom surface. Furthermore, the observed surface heat flow is also used to constrain the distribution of radiogenic heat production, which has large uncertainty in its depth-dependent variation. Because there are not enough heat flow and thermal conductivity observation in the tectonically active areas (Figures 3a, 3c, and 4b), a larger tolerance error (30%) is set in the process of fitting to observed surface heat flow ( Figure S4 in Supporting Information S1) in these areas, especially for the Tibetan Plateau. Whereas for other regions, 20% tolerance error is set, which is consistent with the error of the high-quality observed surface heat flow (Categories A and B; Hu et al., 2000). A constant bottom depth of 100 km of the numerical model in case B is larger than the thermal thickness of the lithosphere in some areas ( Figure S3a in Supporting Information S1). This leads to some temperature differences in the mantle lithosphere in comparison with case C, which just includes the lithosphere. But the spatial pattern of the temperature in cases B and C is similar. The consistent patterns provide a useful reference for us to understand the thermal state of the lithosphere. The following results are based on case C; the results of case B are similar and are shown in Supporting Information S1.

Heat Flow
Through inverse calculation, temperature and heat flow can be obtained for cases B and C. In order to use the surface heat flow to constrain the calculations, we divided the model surface into 1° × 1° regional cells. The calculated heat flow of case C is compared with the average observed heat flow in the regional cell. There are 473 regional cells with observed heat flow data (Figure 5a). We calculated the surface heat flow (Figure 5b) and fitting error (Figure 5c) through inversion. Statistical analysis shows that the fitting error basically obeys the normal distribution (Figure 5d). The root mean square error (RMSE) of all regional cells with heat flow observation is 15.7 mW/m 2 . In comparison, the observational error of the surface heat flow is 20% (Hu et al., 2001). Figure 5. Comparison between the observed (a) and the calculated surface heat flow (b) of case C in the 1° × 1° regional cells. Fitting error (c) is calculated as (observed data-calculated data)/observed data × 100%. Statistical analysis (d) of calculated and observed surface heat flow difference. There are 473 1° × 1° regional cells in the study area and the Root Mean Square Error is 15.7 mW/m 2 for observed and calculated data.
The calculated heat flows in 343 regional cells are with a fitting error within 20%, accounting for 72.5% of the total data ( Figure 5d). The cells with a fitting error greater than 20% are mainly distributed in the Tibetan Plateau and Mongolia Plateau, where larger tolerance error (30%) is set (Figure 5c) due to the lack of measurements and their poor quality. The calculated heat flows are generally lower than the observed values in these areas. However, if the observed heat flow data in the regions of larger tolerance error is not considered, the RMSE gets to 13.8 mW/m 2 in 423 regional cells.
The calculated surface heat flow (Figure 6a) in case C shows several positive anomalies, including the eastern and western Himalayan syntaxes, Lake Baikal, North China Plain, Northeast China, and Southeast China. Compared with the heat flow on the Moho surface (Figure 6b), it can be seen that the proportion of heat flow (Figure 6c) from the mantle in eastern China, including the North China Plain, Northeast China, and Southeast of China, is relatively high. In the Tibetan Plateau, the main contribution of surface heat flow comes from radiogenic heat production in the thick crust. The heat flow from the mantle generally contributes less than 40% of the surface heat flow (Figure 6c). However, for the Precambrian cratons, such as the Siberian Craton, the Indian Plate, Tarim, the NCC and Yangtze block, heat flow from the mantle is all more than 40% of the surface heat flow (Figure 6c). Moreover, in northeast China, North China Plain, and southeast coast of China, heat flow from the mantle contributes more than 50% of the surface heat flow. However, it is emphasized that there is a large uncertainty in the Tibetan Plateau and other tectonic active areas, where uncertainty may raise from assumption of steady-state heat conduction and poor surface heat flow measurements.
The calculated surface heat flow ( Figure S5 in Supporting Information S1) in case B is very similar with that of case C. In 473 1° × 1° regional cells, the RMSE for observed and calculated heat flow is 16.0 mW/m 2 in case B. If the areas of large tolerance error are excluded, in 423 1° × 1° regional cells, the RMSE is 14.1 mW/m 2 for the observed and calculated heat flow. These results show that the heat flow from case B ( Figure S6 in Supporting Information S1) and case C ( Figure 6) have the same pattern and insignificant difference.

Temperature
The calculated temperature in case C at different depths shows great lateral variations. In the crust less than 10 km deep, temperature in the Tibetan Plateau is higher than that of other parts of East Asia (Figures 7a-7b). At a depth of 20-40 km (Figures 7c and 7d), in addition to the Tibetan Plateau and the Tian Shan, high-temperature anomalies also appear in places around the Ordos block, Northeast China, Mongolia Plateau, and Lake Baikal. At the Moho (Figure 7f), the temperature is relatively low, ranging between 500°C, and 700°C, beneath the Precambrian cratons, such as Tarim, NCC, South China Block, and the Indian Plate. In contrast, high temperature (>800°C) is found beneath orogens or regions adjacent to plate boundary zones, including the Tibetan Plateau, Mongolia Plateau, Northeast China, and Trans-North China orogeny. The mantle lithosphere beneath the east part of East Asia, including Northeast China, North China Plain, and the Cathaysia block, shows high temperature at 60 km depth (Figure 7e). It is similar to the temperature of the upper mantle at a depth of 100 km estimated from seismic velocity ( Figure S2b in Supporting Information S1).    (Figures 8c-8f).
We have also compared the temperature results between case B and case C. The pattern of temperature either in the lateral (Figures 7 and S7 in Supporting Information S1) and vertical (Figures 8 and S8 in Supporting Information S1) directions is similar. However, the crustal temperature of case C beneath the stable blocks, such as the Indian Plate, Tarim, and South China Block, is lower than that of case B. This difference results from the different bottom thermal boundary condition. In case B, the bottom surface is set at 100 km, which is deeper than the thermal lithospheric thickness ( Figure S3a in Supporting Information S1) in some areas, and the assumption of steady-state thermal conduction may not apply. Thus, case C is the preferred model for the East Asian continental lithosphere.

The Exponential Attenuation Parameter D
The exponential attenuation parameters D 0 and D 1 of the radiogenic heat production (Equation 6) are inverted for case C (Figure 9). Their values indicate the degree of enrichment of radioactive elements in the crust. The results show that relatively large D 0 and D 1 values (i.e., higher radioactive heating) are mainly found in the Tibetan Plateau, the Ordos, the Mongolia Plateau, and Northeast China. High values are also found in Tian Shan, Shandong Peninsula, and Changbai Shan. In these areas, radiogenic heat production is a major contributor to the high crustal temperature. Some of large D 0 and D 1 values in these areas may be an artifact of a lithosphere thicker than assumed in the model, which reduces mantle heat flow. Thus, an artificially larger D 0 and D 1 of radioactive elements would be required to fit the observed surface heat flow. The D 0 and D 1 values in the calculation also depend strongly on the surface heat flow, which has large uncertainty in the Tibetan Plateau and Mongolia Plateau because of sparse observation and the questionable steady state assumption (Figure 4).

Uncertainty of the Model
This study uses the three-dimensional steady-state heat conduction equation, fully using the available observational thermal data. The model considers the spatial heterogeneity of the parameters and the influence of the Earth's curvature. However, the model has significant limitations.
The near-surface heat flow, thermal conductivity, and heat production are unevenly distributed (Figure 3). Especially in the Tibetan Plateau and Mongolia Plateau, the data are insufficient to constrain the calculation (Figure 4). For example, surface heat flow and thermal conductivity observation in the southwestern and southeastern regions of the Tibetan Plateau are insufficient, which directly affect this area's temperature results. In addition, the uncertainty of thermal parameters has a certain impact on the model results.
Assuming steady-state heat conduction, a change of 5% of the surface heat flow and a 20% change of the average crustal heat production can cause 30°C-50°C and 50°C-70°C variation at 50 km depth, respectively (Artemieva & Mooney, 2001). Upper mantle temperature constrained from seismic data is also associated with large uncertainties. In this model, the temperature of the upper mantle at a depth of 100 km is used as a boundary condition (case B) or an isentropic temperature (1280°C) is set as the bottom of the thermal lithosphere (case C). In case B, the bottom surface may be in the asthenosphere in some regions, therefore the assumption of steady-state heat conduction does not apply, and larger errors may be introduced. Furthermore, the temperature from An and Shi (2006) cannot provide a shallower or crustal temperature boundary condition, because the method of temperature derived from S-wave velocity is established for the composition of the upper mantle (Goes et al., 2000). In case C, the thickness of the thermal lithosphere also has uncertainty ( Figure  S9 in Supporting Information S1). The error of S-wave velocity will cause the temperature variation according to this method (Goes et al., 2000). For instance, S-wave velocity variation of ∼0.15 km/s can cause ∼100°C temperature variation at 130 km depth (An & Shi, 2006). In this study, we used the temperature at 100 km depth in case B. However, the maximum temperature error at 100 km depth can be as high as 150°C ( Figure S2c in Supporting Information S1). The temperature uncertainties caused by the bottom boundary condition get to 0-150°C. This error is not larger than the temperature difference of different blocks. Moreover, in case C, the error of the thermal lithospheric thickness is not more than 30 km ( Figure S9d in Supporting Information S1). Comparison of cases B and C shows little difference in the crustal temperature. Importantly, the relative temperate pattern is the same. For this reason, the relative temperature obtained in this study is more reliable than its actual value.
In the selection of thermal parameters, it is assumed that the lateral variation of thermal parameters observed in the near surface also applies to the entire lithosphere (Equation 3). Such assumptions may not be valid, because the lithology of the mantle lithosphere or even the crust may not be directly related to the near surface material, especially for thermal parameter data from sedimentary layers. In case C, using the thermal conductivity formula from Jaupart and Mareschal (1999) (which did not consider the lateral heterogeneity of thermal properties) instead of Equation 3 (corresponding to taking k 0 = 3.0 W/(m.K) in Equation 3) produces no major changes of the results ( Figure S10 in Supporting Information S1). Temperatures in some regions increase slightly, about 0-50°C. These are the results of using relatively higher thermal conductivity of the mantle comparing to the near surface observation (Figures 3c-3d). Thus, multiple constraints in the inversion calculations provide relatively reliable results. Although the thermal conductivity and radiogenic heat production both affect the results, the results tend to be consistent under the constraints of surface heat flow and temperatures at the top and bottom surfaces.
As mentioned above, all these results are calculated based on assumption of the steady-state heat conduction equation, which is reasonable in Precambrian cratons (Mareschal & Jaupart, 2004). However, most orogens in East Asia have been formed since the India-Asian collision. Due to the strong tectonic activity in some areas of the East Asian continent, it may not be suitable for the steady-state heat conduction assumption. Especially for some intracontinental subduction zones, such as the Himalayan collision zone and the Pamir (Chen et al., 2017;Schneider et al., 2013;Tilmann & Ni, 2003), the thermal state is transient, and the results may provide only a reference for the lithospheric thermal structure. Furthermore, thermal perturbations from both the bottom of the lithosphere and the crust could affect the assumption. If the thermal perturbation comes from the asthenosphere, the steadystate heat conduction assumption is acceptable in some regions because thermal conduction is inefficient in the lithosphere. The character length (L) of conduction is = √ , where κ is the thermal diffusivity (generally, κ = 10 −6 m 2 /s) (Turcotte & Schubert, 2002). For a thermal perturbation at the base of the lithosphere 50 Ma ago (the time of initial Indo-Asian collision), L = ∼40 km. In other words, the thermal signal has not propagated to the crust for a thick lithosphere. However, if the thermal perturbation occurred in the crust (McKenzie & Priestley, 2008;Wang et al., 2013) or fluids are involved (Nábělek & Nábělek, 2014), the assumption of steady-state heat conduction does not apply.
Despite these limitations, the model provides a useful framework of the thermal structure of the East Asian lithosphere. A comparable result is the Curie depth. Previous studies have shown the global Curie depth inversed from the 2-arc-minute resolution Earth Magnetic Anomaly Grid Maus et al., 2009;Xiong et al., 2016). Assuming the Curie isothermal surface temperature is 580°C (Bhattacharyya & Leu, 1975), we compared the Curie point depth obtained by case C with that obtained from the magnetic anomaly data inversion (Figure 10). In some areas, the Curie depth obtained by the two methods is comparable (Figures 10a and 10b), such as in the Indian Plate, the Tarim Basin, and the Yangtze block. However, the difference is relatively large in the southern Tibetan Plateau, Northeast China, and NCC. Besides the lack of magnetic anomaly data (Meyer et al., 2017), part of the reason may be thinner lithosphere beneath these shallower Curie depth areas. This comparison is based on a simple relationship between the temperature and the Curie depth. However, the Curie depth estimated from magnetic data has large differences among different studies. Recent study show that Curie depth differs by as much a ±20 km among different models (Figures 10a  and 10c; Gard & Hasterok, 2021). Furthermore, the Curie temperature of different magnetic minerals varies Figure 10. Comparison between the Curie depth (a and c) and the 580°C isotherm depth in this study (b). (a) Curie depth inversed from the 2-arc-minute resolution Earth Magnetic Anomaly Grid (http://geomag.org/) Maus et al., 2009); (b) The 580°C isotherm depth of this study; (c) The Curie depth based on the equivalent source dipole method (Gard & Hasterok, 2021). greatly. For instance, the Curie temperature of magnetite (Fe 3 O 4 ) is approximately 580°C, while it would be lower when Ti content of titanomagnetite (Fe 2-x Ti x O 3 ) increases (Tanaka et al., 1999). The Curie temperature of partially serpentinized ultramafic bodies, which have metal alloys as the prime magnetic source material, ranges from 620°C to 1100°C (Haggerty, 1978). Thus, the Curie depth from the magnetic anomaly data was not imposed on this numerical model to constrain the results.

Thermal Structure and Tectonic Evolution
The thermal structure of the lithosphere affects its rheology, which controls the deformation of the lithosphere. The results show that high temperature anomalies in the crust are found in Cenozoic orogenic belts, such as the Tibetan Plateau, Tian Shan, Mongolia Plateau, and the peripheral orogenic belts of Ordos (Figures 7a-7d  and 7f). Surface heat flow in these regions is mostly contributed by radiogenic heat production in the thick crust (Figure 6c). The major cratons in East Asia, such as the Tarim, North China Plain, Sichuan Basin, and Indian Plate, show low crustal temperatures. This regional variation is consistent with seismic tomography (Bao et al., 2015;Wei & Zhao, 2020). For these old Precambrian cratons, the contribution of radiogenic heat production to the crustal temperatures is relatively small (Figure 6), consistent with the small D 0 and D 1 values in these regions ( Figure 9). Thus, crustal heat production accounts for a relatively low proportion of the surface heat flow. The results are consistent with the ideas that the cratonic regions may contain relatively little radiogenic heat production due to crustal thickening in their early history (McKenzie & Priestley, 2008Sandiford et al., 2002). Because this process would lead to partial melting of the lower crust and generate granulitic compositions which are relatively depleted in radiogenic heat elements. The existence of lower crustal earthquakes in these cratonic areas can be interpreted to be the result of an anhydrous granulitic lower crust allowing brittle failure at relatively high temperature (Sloan et al., 2011;Zhang et al., 2014).
The NCC has different lithospheric temperature in its eastern and western parts (Figures 7c-7f). The high temperature of the mantle lithosphere in the eastern NCC is mainly caused by lithospheric thinning in the Mesozoic (Chen et al., 2009;He, 2015). The lithospheric temperature of the west part of the NCC, such as the Ordos, is relatively low. Similar internal variations are found in the South China Block. Its southeastern part, the Cathaysia block, has a hot mantle lithosphere. Its northwestern part, The Yangtze block, especially the Sichuan Basin, shows low temperature (Figures 8d and 8e). This suggests that the subduction of the Paleo-Pacific Plate still has a significant impact on the East Asian continent at present. During the subduction rollback (Zhu & Xu, 2019), the rising asthenosphere heated the eastern part of the East Asian continental lithosphere. According to our results, the crust in these regions, such as the North China Plain and South China Block, is still cold. However, for Northeast China, the high mantle temperature caused by the subduction slab may heat the crust because the temperature on the Moho surface is relatively high in that region (Figures 7f and 8f). This is consistent with the results of seismic tomography, which shows low velocity anomalies in the upper mantle He, 2019;Zhao et al., 2017).
As the largest Cenozoic continent-continent collisional orogen, the Tibetan Plateau shows high crustal temperature anomalies (Figure 7). The high crustal temperature in the eastern and western Himalayan syntaxes, Qiangtang terrane, and Songpan-Ganzi terrane corresponds to lower velocity anomalies in the crust . Our results indicate that high temperature regimes in the crust of southeastern Tibetan Plateau are not connected at depth (Figures 7c-7f). But considering the limitations of our model, the results cannot exclude the model of widespread lower crustal flow (Royden et al., 2008). The high temperature anomaly is concentrated on the northeastern boundary of the Tibetan Plateau (Figure 7e), where the Tibetan Plateau has continuously grown laterally (Sun & Liu, 2018). The results in depth-profiles (Figures 8a-8c) show that high-temperature anomalies are found in the mantle lithosphere beyond the stable blocks adjacent to the Tibetan Plateau, such as the Tian Shan to the north of the Tarim, northeastern edge of the Ordos. High mantle lithospheric temperature is consistent with the low seismic velocity anomalies (Li et al., 2018). However, whether or not these results indicate long-range lateral mantle flow (Liu et al., 2004) needs further studies.

Late Cenozoic Volcanoes in East Asia
In the Late Cenozoic, continental volcanoes were widely distributed in East Asia. However, the dynamics of their formation are not clear (Ivanov et al., 2015;Petit & Déverchère, 2006). Their formation is commonly related to the India-Eurasian collision (Molnar & Tapponnier, 1975;Tapponnier & Molnar, 1976) or the Pacific-Eurasia convergence Ward et al., 2021;Zhao et al., 2017). The Late Cenozoic volcanoes on the East and Northeast Asian coast are mainly related to the subduction and stagnation of the Pacific Plate (Northrup et al., 1995;Zhao et al., 2007). However, for the volcanic belts in the Mongolia Plateau and Baikal rift, the mechanism of their formation is controversial (Thybo & Nielsen, 2009;Zhao et al., 2006). Some workers believe that the long-range effect of the India-Eurasian collision activated the strike-slip faults in the southern part of the Mongolia Plateau, which led to the formation of Baikal Rift and the formation of volcanoes in Mongolia Plateau (Gao et al., 1994;Tapponnier & Molnar, 1976). Others have suggested that the volcanoes in the Mongolia Plateau are related to a mantle plume (Chen et al., 2015;Windley & Allen, 1993;Zorin et al., 2003). Recent studies have shown that this mantle plume may be formed because the Paleo-Pacific Plate (Izanagi Plate) sunk into the lower mantle (Ivanov et al., 2015;Kimura et al., 2018;Zhang et al., 2017;Zorin et al., 2006).
Our results suggest that the high-temperature anomalies of the Tibetan Plateau are mainly in the crust. In contrast, the hot mantle lithosphere in east China is likely related to the Pacific subduction or the stagnation of the subducted slabs (Figures 8e-8f). The high crustal temperature anomalies of the Tibetan Plateau are limited by the boundary of Tian Shan-Qilian Shan-Longmen Shan. There is no evidence that the thermal perturbation has extended to central Mongolia or the Baikal rift. Furthermore, geophysical data and evidence from xenoliths in volcanic rocks indicate no evidence for a large-scale lower mantle upwelling and crustal thinning beneath the Baikal rift and central Mongolia Plateau (Ionov, 2002;Thybo & Nielsen, 2009), which implies that a mantle plume does not drive them. Thus, our study shows that thermal perturbation associated with the subduction of the Pacific plate is the dominant cause for the continental volcanoes in East Asia. The isolines of the Benioff zone depth of the Pacific Plate and the frontier of the stagnant slab (Liu et al., 2017;Zhao, 2021) are spatially correlated with the distribution of volcanoes along the Baikal rift-central Mongolia Plateau-North China Plain (Figures 1 and 8). This indicates that the asthenospheric disturbance caused by the subducted and stagnated Pacific Plate may have affected the Mongolian Plateau and Lake Baikal, and is not be limited to the Northeast China as previously thought (Xu et al., 2020).

Conclusions
Using the newly available measurements of surface heat production that cover much of mainland China and southern Mongolia, we developed a three-dimensional model of lithospheric thermal structure for the East Asian continent. This model is a major improvement over the exiting thermal model (An & Shi, 2007;Sun et al., 2013). Major results include the following.
1. The crust of the Cenozoic orogenic belts in East Asia shows high-temperature anomalies, especially in areas with strong Cenozoic tectonism, such as the Tibetan Plateau, the peripheral area of the Ordos block, and Lake Baikal. 2. The Tibetan Plateau has a high temperature crust, mainly resulting from high radiogenic heat production in the thickened crust. The bounding blocks, including the Tarim, Ordos and Sichuan Basin, show a relatively cold lithosphere. High temperature anomalies are found in the Tian Shan, northeastern margins of the Ordos, and central Mongolia. 3. Subduction of the Paleo-Pacific and Pacific Plates has a significant impact on the thermal structure of the eastern blocks of the East Asian continent. The distribution of Late Cenozoic volcanoes in East Asia correlates with high temperature anomalies in the mantle lithosphere, indicating that Pacific subduction is the main cause for continental volcanism in East Asia.