Stratiﬁed Data Reconstruction and Spatial Pattern Analyses of Soil Bulk Density in the Northern Grasslands of China

: The spatial pattern of soil bulk density in the grasslands of northern China largely remains undeﬁned, which raised uncertainty in understanding and modeling various soil processes in large spatial scale. Based on the measured data of soil bulk density available from soil survey reports from the grasslands of northern China, we constructed a soil Stratiﬁed Pedotransfer function (SPTF) from the surface soil bulk density. Accordingly, the stratiﬁed bulk density data of soil vertical proﬁle was reconstructed, and the estimation of soil bulk density data in horizontal space was performed. The results demonstrated that the soil bulk density of the grasslands of northern China was typically high in the central and northwestern regions and low in the eastern and mountainous regions. Mean soil bulk density of the grasslands was 1.52 g · cm − 3 . According to geographical divisions, the highest soil bulk density was observed in the Tarim basin, with mean soil bulk density of 1.91 g · cm − 3 . Conversely, the lowest soil bulk density was observed in the Tianshan Mountain area, with mean soil bulk density of 1.01 g · cm − 3 . Based on data obtained on various types of grasslands, the soil bulk density of alpine meadow was the lowest, with a mean soil bulk density of 0.75 g · cm − 3 , whereas that of temperate desert was the highest, with mean soil bulk density of 1.80 g · cm − 3 . Mean prediction error, root mean square deviation, relative error, and multiple correlation coe ﬃ cient of soil bulk density data pertaining to surface layer (0–10 cm) in the grasslands of northern China were 0.018, 0.223, 16.2%, and 0.5386, respectively. The approach of employing multiple data sources via soil transfer function improved the estimation accuracy of soil bulk density from stratiﬁed soils data at the large scale. Our study would promote the accurate assessment of grassland carbon storage and ﬁne land characteristics mapping.


Introduction
Soil, an essential component of terrestrial ecosystem, is fundamental for the vegetation survival. Its physical and chemical properties affect the growth of plants, as well as restrict their productivity [1,2]. Soil bulk density, defined as the soil mass (or weight) per unit volume of undisturbed soil column [3], is an important physical property of soil that has a critical impact on soil permeability, infiltration, water-holding capacity, solute transport, and soil erosion resistance [4,5]. Hence, it quantitatively characterizes the ecological functions of soil and is one of the major indicators for evaluating the

Soil Bulk Density Survey and Soil Sampling
Two soil datasets were collected in the present study (Table 1). One dataset was obtained from the "Investigation on the Degradation and Causes of Grassland Resources in Key Pastoral Areas of Temperate Grassland in China" project from 2013 to 2015. Part of the soil bulk density data in dataset was applied in SPTF model validation. Axillary data (environment data, remote sensing data and vegetation data) was applied in spatial interpolation. The surface soil bulk density survey and soil sampling of grassland of northern China were conducted together with the grassland vegetation sampling survey. The establishment of grassland vegetation sample sites was based on the distribution characteristics including grassland type, utilization mode and intensity, and road traffic conditions among others. The survey route with relatively wide coverage was set up in the survey area, and the landscape information of grassland vegetation along the along the survey route was always collected using GPS (Global Positioning System) positioning instrument (American MAGELLAN, eXplorist 610, position errors less than 15 m, position errors less than 15 m) and camera (Canon, EOS 2000D, Japan, Oita). Simultaneously, one grassland vegetation community survey site of 100 m × 100 m should be set up every 50 km along the survey route. However, the setting of grassland sample site was determined in accordance with the specific conditions of grassland type, utilization mode and intensity, and site conditions. For example, grassland sample site was set for several kilometers in mountainous regions, whereas it was set at intervals of 80 km in plain regions with relatively simple type of grassland. The present study included 6716 grassland landscape information collection sites and 587 grassland survey sites. The field survey of grassland was completed within a period of 3 years from 2013 to 2015, and the sample survey and sampling of grassland were performed from July to September every year. The longitude and latitude, elevation, utilization mode, and landscape photos of each sample site were recorded. Three 1 m × 1 m sites were

Soil Bulk Density Survey and Soil Sampling
Two soil datasets were collected in the present study (Table 1). One dataset was obtained from the "Investigation on the Degradation and Causes of Grassland Resources in Key Pastoral Areas of Temperate Grassland in China" project from 2013 to 2015. Part of the soil bulk density data in dataset was applied in SPTF model validation. Axillary data (environment data, remote sensing data and vegetation data) was applied in spatial interpolation. The surface soil bulk density survey and soil sampling of grassland of northern China were conducted together with the grassland vegetation sampling survey. The establishment of grassland vegetation sample sites was based on the distribution characteristics including grassland type, utilization mode and intensity, and road traffic conditions among others. The survey route with relatively wide coverage was set up in the survey area, and the landscape information of grassland vegetation along the along the survey route was always collected using GPS (Global Positioning System) positioning instrument (American MAGELLAN, eXplorist 610, position errors less than 15 m, position errors less than 15 m) and camera (Canon, EOS 2000D, Japan, Oita). Simultaneously, one grassland vegetation community survey site of 100 m × 100 m should be set up every 50 km along the survey route. However, the setting of grassland sample site was determined in accordance with the specific conditions of grassland type, utilization mode and intensity, and site conditions. For example, grassland sample site was set for several kilometers in mountainous regions, whereas it was set at intervals of 80 km in plain regions with relatively simple type of grassland. The present study included 6716 grassland landscape information collection sites and 587 grassland survey sites. The field survey of grassland was completed within a period of 3 years from 2013 to 2015, and the sample survey and sampling of grassland were performed from July to September every year. The longitude and latitude, elevation, utilization mode, and landscape photos of each sample site were recorded. Three 1 m × 1 m sites were set up along the diagonal line of each grassland sample site for evaluating the grassland community index, including data of main species and species in the sample site, the mean height of grassland, coverage, aboveground biomass of species, and underground biomass [38]. In the present study, the surface soil bulk density survey and soil sample sampling were conducted in the three sample sites used in the grassland community index survey. Surface soil bulk density was sampled using cutting ring method. Three samples were collected along the diagonal line in each sample plot. The sampling site was in the center of 0-10 cm soil layer. The cutting ring knives were numbered and brought back to the room. The samples were dried to constant weight in an oven at 105 • C. The measured weight converted into standard unit (g·cm −3 ). The specifications of the cutting ring knife were 100 cm 3 . Physical and chemical analyses for the soil samples were performed using soil drilling method. In the sample sites used for the grassland community index survey, soil drills with 33 cm inner diameter were used to stratify the soil layers of 0-10 cm, 10-20 cm, 20-30 cm, and 30-40 cm, respectively. Each sample was assessed with 3 drills, following which they were mixed with layers according to the sample number. The samples were packed in plastic bag and brought indoor. The soil samples were spread and dried, mixed evenly, and analyzed using quadruple method. After grinding, the samples were screened using 0.3-mm separate nylon mesh bags and were packed into small plastic bags, numbered, and sent to the Botanical Research Center of the Chinese Academy of Sciences for analysis of total carbon and organic carbon content. The carbon content was determined using dry burning method, and the instrument was vario macro cube (German, Manchester).
The other dataset was collected from 143 grassland sample sites involved in the Pilot Project of the Chinese Academy of Sciences for the establishment of a Pedotransfer function model of stratified soil bulk density (Table 1). It was used to construct stratified model. Sampling methods of soil bulk density and soil analysis samples were based on the guidelines for the Investigation of Carbon Fixation Status, Rate and Potential of Grassland Ecosystem in China. A 100 cm (length) × 50 cm (width) × 100 cm (depth) trench was excavated on the side of the sample by trench method. According to the soil profile, stratified samples were obtained from 0-5 cm, 5-10 cm, 10-20 cm, 20-30 cm, and 30-50 cm layers, and repeated sampling was performed for five samples. This method is extremely time-consuming and laborious to be implemented in the large-scale survey of the northern temperate grasslands. The soil samples were spread out and dried, evenly mixed, and analyzed indoors using quadruple method. After grinding, the samples were screened using a 100-mesh sieve and packed into small plastic bags, numbered, and sent to the Botanical Research Center of the Chinese Academy of Sciences for analysis of total carbon and organic carbon content.

Spatial Data Collection and Processing on Geographical Elements
We proposed a new soil stratified Pedotransfer function from the surface soil bulk density, which help construct the missing stratified data in a large scale, and used MWRM (multi-factor weighted regression model) model [39] to interpolate soil bulk density of temperate grassland of northern China.
Above that, auxiliary data including remote sensing data, climatic data, topographic data, and resource distribution of grassland types, were required (Table 1).
(i) Remote sensing data: normalized difference vegetation index (NDVI) data from mid-August 2013 to mid-August 2015 of northern temperate grassland were obtained from MOD13A2 of the United States Geological Survey, which was tiled and projected via MRT (MODIS Reprojection Tool, https://lpdaac.usgs.gov/tools/modis_reprojection_tool/) application with a spatial resolution of 1 km. (ii) Climatic data: according to the observation data of meteorological centers (stations) from 2013 to 2015 provided by the National Meteorological Administration, a professional software for meteorological interpolation, ANUSPLIN (http://fennerschool.anu.edu.au/research/products/ anusplin-vrsn-44), was used to estimate spatial interpolation of annual mean temperature, annual rainfall, accumulated temperature (>10 • C), humidity (Ivanov humidity) [36] and other spatial variables, with a spatial resolution of 1 km. (iii) Altitude data: data was collected from Digital Elevation Model (DEM) in National Earth System Science Data Sharing Platform (www.geodata.cn), with a spatial resolution of 30 m. (iv) Spatial distribution of grassland resource map: based on the vector data for distribution mapping of grassland types obtained from the national grassland survey in 1980, which is calibration value for compiling and revising the spatial map of grassland of 2010 using TM (Thematic mapper) data, with a spatial resolution of 20 m. (v) Matching and extraction of geographic element data of grassland sample sites: based on the longitudinal and latitudinal coordinates of 587 grassland sample sites in the northern temperate grassland region, the geographic element spatial data mentioned above were loaded in ArcGIS, and the values for 6 elements, including altitude, NDVI, annual mean temperature, annual rainfall, accumulated temperature (>10 • C), humidity, and surface soil bulk density of those 587 sample sites, were extracted. A dataset (two-thirds of them totaled 397 samples data as modeling data, and the remaining 190 samples data were used as test data) consisting of these six ecological factors, surface soil bulk density, and soil organic carbon content was constructed for statistical, regression, and spatial grid analyses.

Statistics and Regression Analysis
All collected data related to soil bulk density and organic carbon content in sample sites were analyzed using SPSS 20.0 software, including correlation analyses for soil organic carbon content with soil depth (K SOC ), correlation analyses for soil bulk density with soil depth (K SBD ), regression analyses for association between surface soil bulk density and geographical factors, as well as correlation analyses for coefficient of variation in vertical profile of soil bulk density (K SBD ) and coefficient of variation in vertical profile of soil organic carbon content (K SOC ).

Construction of Stratified Soil Bulk Density from Pedotransfer Function
According to the Figure 2, the soil bulk density dataset of 143 grassland sample sites from "Strategy Pilot Project of Chinese Academy of Sciences" was used for stratified model after K S normality test (p < 0.01) which was performed using SPSS 20.0 software, and 117 verified sample sites were eventually obtained for statistical analyses after excluding sites with obvious anomalies. Based on the stratified data of soil bulk density in 117 verified sample sites, the variability coefficients of vertical profile for soil bulk density and soil organic carbon content in each sample site, named K SBD and K SOC , respectively, could be obtained by using linear regression analysis (Equation (1)). Furthermore, 397 grassland sample sites from "Investigation on the Degradation and Causes of Grassland Resources in Key Pastoral Areas of Temperate Grassland in China" were used to constructed spatial model (Equation (2)). Based on the SPTF Equations (1) and (2), soil bulk density in different soil layers (0-10cm, 10-20cm, 20-30cm, 30-50cm) was estimated.
In the equation, a and b are model constant terms; SBD(x) is the soil bulk density of a soil layer. When x = 0, it is the soil bulk density of the surface layer (0-10 cm); when x = 1, it is the soil bulk density of 10-20 cm; when x = 2, it is the soil bulk density of 20-30 cm; when x = 3, it is the soil bulk density of 30-50 cm; and K SBD is the variability coefficient of vertical profile for soil bulk density.
ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 7 of 16 In the equation, a and b are model constant terms; SBD(x) is the soil bulk density of a soil layer. When x = 0, it is the soil bulk density of the surface layer (0−10 cm); when x = 1, it is the soil bulk density of 10−20 cm; when x = 2, it is the soil bulk density of 20−30 cm; when x = 3, it is the soil bulk density of 30−50 cm; and KSBD is the variability coefficient of vertical profile for soil bulk density.

Spatialization of Soil Bulk Density
The MWRM model [39] was established using the surface soil bulk density data obtained from the grassland survey, including 6 geographical factors of elevation, annual mean temperature, annual mean rainfall, accumulated temperature (≥ 10°C), humidity, and NDVI. Here, we used MWRM model to interpolate the spatial stratified pattern of soil bulk density ( Figure 2). We got the spatial grid data (1 km × 1 km) of surface soil bulk density (0-50 cm) according to the SPTF equation (1) and (2). After that, the stratified spatial pattern of soil bulk density (0-10cm, 10-20cm, 20-30cm, 30-50cm) was estimated by MWRM model with auxiliary data and spatial grid data.

Model Validation and Accuracy Evaluation
In order to evaluate the estimation accuracy of sample points and model results, the average prediction error (MPE), root mean square error (RMSE), relative error (E), and multiple correlation coefficient (R 2 ) were applied. In equation (3)-(6), , , , and represent the measured value, predicted value, average predicted value, and average value, respectively, and n is the number of

Spatialization of Soil Bulk Density
The MWRM model [39] was established using the surface soil bulk density data obtained from the grassland survey, including 6 geographical factors of elevation, annual mean temperature, annual mean rainfall, accumulated temperature (≥ 10 • C), humidity, and NDVI. Here, we used MWRM model to interpolate the spatial stratified pattern of soil bulk density ( Figure 2). We got the spatial grid data (1 km × 1 km) of surface soil bulk density (0-50 cm) according to the SPTF Equations (1) and (2). After that, the stratified spatial pattern of soil bulk density (0-10cm, 10-20cm, 20-30cm, 30-50cm) was estimated by MWRM model with auxiliary data and spatial grid data.

Model Validation and Accuracy Evaluation
In order to evaluate the estimation accuracy of sample points and model results, the average prediction error (MPE), root mean square error (RMSE), relative error (E), and multiple correlation coefficient (R 2 ) were applied. In Equations (3)-(6), P i , M j , M j , and T represent the measured value, predicted value, average predicted value, and average value, respectively, and n is the number of samples. A total of 190 sample points was used to model validation, to compare with average value of different layers of soil bulk density (0-10 cm, 10-20 cm, 20-30 cm, 30-50 cm).

Spatial Distribution and Patterns of Stratified Soil Bulk Density of Grassland in Northern China
Based on the data of surface soil bulk density available from soil survey reports, with the geographic element spatial data mentioned above including altitude, NDVI, annual mean temperature, annual rainfall, accumulated temperature (>10 • C), humidity, the spatial distribution of surface soil (0-10 cm) soil bulk density (Figure 3a) was established using MWRM model [40]. A total of 190 sample points was used to model validation, to compare with average value of different layers of soil bulk density (0-10 cm, 10-20 cm, 20-30 cm, 30-50 cm). The spatial distribution of soil bulk density in the northern grasslands were established by extracting and removing the non-grassland geographical units, such as desert and farmland.

Spatial Distribution and Patterns of Stratified Soil Bulk Density of Grassland in Northern China
Based on the data of surface soil bulk density available from soil survey reports, with the geographic element spatial data mentioned above including altitude, NDVI, annual mean temperature, annual rainfall, accumulated temperature (>10 °C), humidity, the spatial distribution of surface soil (0−10 cm) soil bulk density (Figure 3a) was established using MWRM model [40]. A total of 190 sample points was used to model validation, to compare with average value of different layers of soil bulk density (0-10 cm, 10-20 cm, 20-30 cm, 30-50 cm). The spatial distribution of soil bulk density in the northern grasslands were established by extracting and removing the non-grassland geographical units, such as desert and farmland. Based on the validation analysis for 190 test sites (Figure 4) in the northern grassland region, the mean prediction error (MPE), root mean square error (RMSE), relative error, and R 2 of the surface soil bulk density (0−50 cm) was 0.018, 0.223, 16.2%, and 0.5386, respectively. The testing results indicated that the spatial distribution of surface soil bulk density of grassland surface layer (0−10 cm) using MWRM model [39] is highly reliable and accurate, which could be used to calculate the spatial data of surface soil bulk density. Based on the validation analysis for 190 test sites (Figure 4) in the northern grassland region, the mean prediction error (MPE), root mean square error (RMSE), relative error, and R 2 of the surface soil bulk density (0-50 cm) was 0.018, 0.223, 16.2%, and 0.5386, respectively. The testing results indicated that the spatial distribution of surface soil bulk density of grassland surface layer (0-10 cm) using MWRM model [39] is highly reliable and accurate, which could be used to calculate the spatial data of surface soil bulk density.  According to the spatial pattern of stratified soil bulk density (Figure 3a, b, c, d), the spatial variation of soil bulk density in different soil layers is consistent, which is high in central and northwest in grasslands of northern China (mainly distributed the temperate desert and temperate desert steppe), and low in the east and southwest (mainly distributed the temperate meadow grassland and montane meadow).
Regarding statistical analysis in Figure 5, the value for soil bulk density in grasslands of northern China, accounting for 61.58% of the total grassland area, was normally >1.5 g·cm -3 . Class I soil bulk density (0-0.8 g·cm -3 ) was sporadically distributed in Tianshan Mountain and Altai Mountains, and the main grassland type was alpine meadow, accounting for 4.54% of the total grassland area. Class II soil bulk density (0.8-1.2 g·cm -3 ) was concentrated in Tianshan Mountain, Altai Mountains, Da Hinggan Ling Prefecture and Amir-Kunlun-Altun Plateau Region, and the main grassland type was montane meadow, accounting for 13.16% of the total grassland area. Class III soil bulk density (1.2-1.5 g·cm -3 ) was concentrated in the eastern of Inner Mongolia Plateau and Hulunbuir Plateau, and the main grassland type was temperate meadow grassland and montane meadow, accounting for 20.74% of the total grassland area. Class IV soil bulk density (1.5-1.8 g·cm -3 ) was concentrated in the western of Inner Mongolia Plateau, northern of the Junggar Basin and Loess Plateau and edge of Alxa Plateau, and the main grassland type was temperate steppe and temperate desert steppe, accounting for 36.03% of the total grassland area. Class V soil bulk density (>1.8 g·cm -3 ) was concentrated in the southern of Junggar Basin, edge of Tarim basin and central of Alxa Plateau, and the main grassland type was temperate desert, accounting for 25.55% of the total grassland area. In addition, with the increase of soil depth, the proportion of class V soil bulk density grassland has increased. Compared with the surface layer (0-10 cm), the total area of class V soil bulk density grassland increased by 21.27% in 30-50 cm soil layer. Meanwhile, the proportion of class I, II, III and IV soil bulk density (0-1.8 g·cm -3 ) grassland has decreased, and the total area of 30-50 cm soil layer decreased by 3.3%, 5.06%, 7.48%, 5.43%, respectively, compared with surface soil (0-10 cm). According to the spatial pattern of stratified soil bulk density (Figure 3a-d), the spatial variation of soil bulk density in different soil layers is consistent, which is high in central and northwest in grasslands of northern China (mainly distributed the temperate desert and temperate desert steppe), and low in the east and southwest (mainly distributed the temperate meadow grassland and montane meadow).
Regarding statistical analysis in Figure 5, the value for soil bulk density in grasslands of northern China, accounting for 61.58% of the total grassland area, was normally >1.5 g·cm −3 . Class I soil bulk density (0-0.8 g·cm −3 ) was sporadically distributed in Tianshan Mountain and Altai Mountains, and the main grassland type was alpine meadow, accounting for 4.54% of the total grassland area. Class II soil bulk density (0.8-1.2 g·cm −3 ) was concentrated in Tianshan Mountain, Altai Mountains, Da Hinggan Ling Prefecture and Amir-Kunlun-Altun Plateau Region, and the main grassland type was montane meadow, accounting for 13.16% of the total grassland area. Class III soil bulk density (1.2-1.5 g·cm −3 ) was concentrated in the eastern of Inner Mongolia Plateau and Hulunbuir Plateau, and the main grassland type was temperate meadow grassland and montane meadow, accounting for 20.74% of the total grassland area. Class IV soil bulk density (1.5-1.8 g·cm −3 ) was concentrated in the western of Inner Mongolia Plateau, northern of the Junggar Basin and Loess Plateau and edge of Alxa Plateau, and the main grassland type was temperate steppe and temperate desert steppe, accounting for 36.03% of the total grassland area. Class V soil bulk density (>1.8 g·cm −3 ) was concentrated in the southern of Junggar Basin, edge of Tarim basin and central of Alxa Plateau, and the main grassland type was temperate desert, accounting for 25.55% of the total grassland area. In addition, with the increase of soil depth, the proportion of class V soil bulk density grassland has increased. Compared with the surface layer (0-10 cm), the total area of class V soil bulk density grassland increased by 21.27% in 30-50 cm soil layer. Meanwhile, the proportion of class I, II, III and IV soil bulk density (0-1.8 g·cm −3 ) grassland has decreased, and the total area of 30-50 cm soil layer decreased by 3.3%, 5.06%, 7.48%, 5.43%, respectively, compared with surface soil (0-10 cm).

Comparison of Stratified Soil Bulk Density in Different Geographic Regions
Considering the analyses of main geographical distribution ( Figure 6 and Table 2), the distribution of stratified soil bulk density in different geographical regions showed a trend from low to high. Mean soil bulk density of grassland in the grasslands of northern China was 1.52 g·cm -3 , the maximum soil bulk density is 2.23 g·cm -3 while the minimum soil bulk density was 0.14 g·cm -3 . The lowest soil bulk density value was observed in Tianshan Mountain area, with a mean soil bulk density of 1.01 g·cm -3 and surface soil bulk density was 0.94 g·cm -3 . Conversely, the highest soil bulk density value was observed in the grassland of Tarim basin, with a mean soil bulk density of 1.91 g·cm -3 and surface soil bulk density of 1.85 g·cm -3 . The results showed that the high value of soil bulk density mainly distributed in plain and basin, which was arid climate, and the low value of soil bulk density mainly distributed in mountain and plateau, which was humid climate.

Comparison of Stratified Soil Bulk Density in Different Geographic Regions
Considering the analyses of main geographical distribution ( Figure 6 and Table 2), the distribution of stratified soil bulk density in different geographical regions showed a trend from low to high. Mean soil bulk density of grassland in the grasslands of northern China was 1.52 g·cm −3 , the maximum soil bulk density is 2.23 g·cm −3 while the minimum soil bulk density was 0.14 g·cm −3 . The lowest soil bulk density value was observed in Tianshan Mountain area, with a mean soil bulk density of 1.01 g·cm −3 and surface soil bulk density was 0.94 g·cm −3 . Conversely, the highest soil bulk density value was observed in the grassland of Tarim basin, with a mean soil bulk density of 1.91 g·cm −3 and surface soil bulk density of 1.85 g·cm −3 . The results showed that the high value of soil bulk density mainly distributed in plain and basin, which was arid climate, and the low value of soil bulk density mainly distributed in mountain and plateau, which was humid climate.

Comparison of Stratified Soil Bulk Density in Different Geographic Regions
Considering the analyses of main geographical distribution ( Figure 6 and Table 2), the distribution of stratified soil bulk density in different geographical regions showed a trend from low to high. Mean soil bulk density of grassland in the grasslands of northern China was 1.52 g·cm -3 , the maximum soil bulk density is 2.23 g·cm -3 while the minimum soil bulk density was 0.14 g·cm -3 . The lowest soil bulk density value was observed in Tianshan Mountain area, with a mean soil bulk density of 1.01 g·cm -3 and surface soil bulk density was 0.94 g·cm -3 . Conversely, the highest soil bulk density value was observed in the grassland of Tarim basin, with a mean soil bulk density of 1.91 g·cm -3 and surface soil bulk density of 1.85 g·cm -3 . The results showed that the high value of soil bulk density mainly distributed in plain and basin, which was arid climate, and the low value of soil bulk density mainly distributed in mountain and plateau, which was humid climate.

Comparison of Stratified Soil Bulk Density in Different Grassland Types
According to the data of grassland types in the grassland of northern China [36], soil bulk density data of main grassland types were extracted, and statistical analysis of different types of stratified soil bulk density was conducted (Table 3 and Figure 7). The result revealed that the soil bulk density of different grassland types varied significantly, among which the alpine meadow had the lowest soil bulk density, with mean soil bulk density of 0.75 g·cm −3 and surface soil bulk density of 0.68 g·cm −3 . The soil bulk density of temperate desert had the highest value which was 1.80 g·cm −3 and mean surface soil bulk density was 1.72 g·cm −3 . The results showed that the high value of soil bulk density mainly distributed in desert which was arid climate, and the low value of soil bulk density mainly distributed in steppe and meadow which was humid climate. bulk density of different grassland types varied significantly, among which the alpine meadow had the lowest soil bulk density, with mean soil bulk density of 0.75 g·cm -3 and surface soil bulk density of 0.68 g·cm -3 . The soil bulk density of temperate desert had the highest value which was 1.80 g·cm -3 and mean surface soil bulk density was 1.72 g·cm -3 . The results showed that the high value of soil bulk density mainly distributed in desert which was arid climate, and the low value of soil bulk density mainly distributed in steppe and meadow which was humid climate.

Relationship Between Grassland Type, Organic Carbon Content and Soil Bulk Density
The change of soil bulk density is not only influenced by soil properties, such as texture, compactness, and land use type [2,41], but also by soil physical and chemical properties, such as clay content, water content, silt content, depth, PH value, and organic carbon content [42]. However, the most important factor influencing the change of soil bulk density is organic carbon content [43,44]. The content of soil organic carbon is comprehensive result of organic matter input and mineralization, and, for the grassland ecosystem, the input of organic matter mainly relies on the photosynthesis of herb and litter of small shrub. However, the ability of organic matter input for different types of grassland shows difference, causing differences in organic carbon content and thus affecting the soil bulk density.
On the other hand, the formation and distribution of soil types are adapted to the bioclimatic zone, and the soil bulk density of different grassland types varies widely. For grassland in the study area, there are three main types of grassland: temperate steppe, temperate desert steppe, and temperate meadow grassland. Three compatible soil types, namely chestnut soil, brown soil, and chernozem soil, were formed. The soil physical and chemical properties of the three soil types, such as texture and water content, are significantly different, and these differences affect the distribution of vegetation in turn. In general, soil porosity of brown soil is greater than that of chernozem soil, resulting in increased soil hardness and reduced soil water holding capacity, which will increase soil bulk density. Céspedes [41] found that the increase of soil bulk density would significantly affect the accumulation and holding capacity of soil organic carbon content. The change of grassland type leads to the increase or decrease of soil bulk density, which affects the soil organic carbon content also.

Accuracy Analyses of Soil Bulk Density Data
According to the comparative analysis of soil bulk density research results of different regions by different scholars, the soil bulk density stratification data estimated in this study is close to the results of various scholars (Table 4), which basically reflects the distribution of grassland soil bulk density spatial distribution pattern in the temperate steppe region of northern China. But there are significant differences in individual regions. In this study, we collected the relevant studies of seven main geographical areas in the study area. However, because of the difficulties and uncertainties of soil bulk density sampling, along with the variation and heterogeneity of sample site settings, the prediction models for soil bulk density established by different scholars were rather diverse, leading to differences in soil bulk density measurements. Zhou [37] employed the multivariate weighted regression model (MWRM) [39] for the prediction of 198 sampling sites in Yili Xinjiang and obtained extremely different results. This difference might be attributed to the fact that Zhou [37] did not estimate the soil bulk density based on the soil depth of Yili area. Other studies have shown that there are significant differences in soil bulk density between different soil depth ranges in the vertical direction in the Yili region [45]. Moreover, because of the special topography of Yili region, the soil layer contains more gravel, which increases the variability of soil bulk density in the vertical direction. Zhang [45] found that the presence of gravel can increase soil bulk density by approximately 16%, which indicates that the difference in soil composition of 0-30 cm in Yili area will cause significant changes in soil bulk density. Cheng [46] mainly measured the soil bulk density of 0~30cm in the montane meadow of the northern slope of the Tianshan Mountains, and the number of samples was small, which did not fully reflect the actual situation of the Tianshan Mountains soil. Wang [47] performed soil bulk density analyses at 27 sampling sites of Hulunbeir grassland, and the results were significantly different from those obtained in the present study. Although this difference may be attributed to the limited number of sampling points that may not fully reflect the overall soil characteristics of the region, it is also noteworthy that Wang [47] selected samples of desertification gradient from Hulunbeir region without sampling other types of grassland. This difference may also be affected by the boundary of the study area, soil heterogeneity, and so on.

Application and Problems of Pedotransfer Function
China has a vast territory, complex and changeable natural geographical environment, and the relationship between soil bulk density and other soil properties in different soil types could be considerably diverse. Therefore, most studies have fully considered the territoriality of soil when applying a Pedotransfer function model for soil bulk density [50][51][52]. Moreover, based on the classification of soil system, Han [50] grouped the existing soil data of China and established a polynomial model for each soil subtype. Although this approach improved the accuracy of the model for prediction to a certain extent, it was not suitable for complex soil environment and did not consider the variation trend of soil bulk density in vertical profile. On the other hand, the models proposed in other studies, such as the previous empirical models, did not consider factors including topography, climate, and parent materials [53,54]. Therefore, it is necessary to construct more accurate soil bulk density Pedotransfer function for large-scale prediction of soil bulk density.
In the present study, a stratified Pedotransfer function model for soil bulk density in the grasslands of northern China was established using the complete soil bulk density profile data from a grassland survey conducted by the Pilot Project of Chinese Academy of Sciences. The impact of regional differences and vertical profile on soil bulk density was completely considered, in combination with climate and topography among other factors. The accuracy of soil bulk density prediction was improved to a certain extent.

Conclusions
In this study, based on the sampled data, we simulated the spatial pattern of soil bulk density 0-50 cm in northern China grassland using MWRM model and stratified soil Pedotransfer function. The result indicated that: (i) The soil bulk density in grassland of northern China was high in the central and northwestern regions and low in the eastern regions. The mean soil bulk density of grassland was 1.52 g·cm −3 . (ii) Soil bulk density shows spatial heterogeneity, and its distribution trend is consistent with grassland type. (iii) According to the validation results (MPE = 0.018, RMSE = 0.223, relative error = 16.2%, R2 = 0.5386), our Pedotransfer function provides insights into the construction of stratified data for soil bulk density estimation at a large scale, which raised accuracy of soil organic matter calculation.