• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Spatial variability of soil organic carbon in the forestlands of northeast China

    2014-09-06 11:13:31LingLiuHaiyanWangWeiDaiXiangdongLeiXiaojuanYangXuLi
    Journal of Forestry Research 2014年4期

    Ling Liu · Haiyan Wang · Wei Dai · Xiangdong Lei Xiaojuan Yang · Xu Li

    ORIGINAL PAPER

    Spatial variability of soil organic carbon in the forestlands of northeast China

    Ling Liu · Haiyan Wang · Wei Dai · Xiangdong Lei Xiaojuan Yang · Xu Li

    Received: 2013-10-17; Accepted: 2014-03-06

    ? Northeast Forestry University and Springer-Verlag Berlin Heidelberg 2014

    Soil organic carbon (SOC) is an effective indicator of soil fertility and productivity, and it varies spatially and temporally in relation to other soil properties. Spatial variability of SOC in the forestlands of northeast China was characterized using geostatistics. Soil samples at the depths of 0-20 cm, 20-40 cm and 40-60 cm were collected from sixty-three temporary plots to evaluate SOC concentration and density (SOCD) and other soil properties. We analyzed correlations between SOC and soil properties. Soil organic carbon concentrations were high. The total amount of C stored in soil (0-60 cm) was 16.23 kg·m-2with the highest SOCD of 7.98 kg·m-2in topsoil. Soil properties in most cases differed by horizon, suggesting different processes and effects in each horizon. Soil organic carbon had positive relationships with total N, P and K as well as readily available K, but did not show a significant positive correlation with available P. Spatial factors including elevation, slope and aspect affected SOC distribution. Soil organic carbon at 0-60 cm had strong spatial autocorrelation with nugget/sill ratio of 5.7%, and moderate structured dependence was found at 0-20 cm, which indicatedthe existence of a highly developed spatial structure. Spatial distributionsof SOC concentration and SOCD were estimated using regression-kriging, with higher prediction accuracy than ordinary kriging. The fractal dimension of SOC indicated the preferential pattern of SOC distribution, with the greatest spatial heterogeneity and strongest spatial dependence in the northeast-southwest direction.

    northeast China, soil organic carbon, spatial variability, geostatistics

    Introduction

    Soil organic carbon (SOC) is a key determinant of soil quality (Stocking 2003) and represents one of the largest reservoirs of organic carbon in the global carbon cycle that can influence global warming (Lützow et al. 2006). In terrestrial ecosystems, forest soil carbon is the largest carbon sink. Soil organic carbon storage in global forest ecosystems was estimated as 402×1012-787×1012kg, or 25%-50% of the global soil carbon storage (Arrouays and Pelissier 1994). In northeast China, considerable organic carbon is stored in soils because of the dense forest vegetation and the thick accumulation of organic matter. Humic cambisols (Dark brown forest soils, in Chinese soil taxonomy) are typical soils mainly found in warm-temperate, mixed coniferous and broad-leaved forests and cold-temperate coniferous forests in northeast China. The soils were higher in natural fertility and supported dense vegetation. The soils of 40.2×106ha are mainly located in the provinces of Heilongjiang, Jilin and Inner Mongolia (Zhao et al. 2004). These areas have been, therefore, among the most important for timber production and as a sink for organic matter on a global scale. Humic cambisols are also considered to store substantial amounts of organic matter, thereby serving as CO2sinks.

    Geostatistics is an effective method for study of spatial distribution characteristics and their variability. Geostatistical methods quantify spatial distributions based on the spatial scale of the study area, distance between detecting points, and spatial patterns of modeling semivariograms. They have been widely ap-plied to evaluate spatial correlation in soils and to analyze the spatial variability of soil characteristics, such as soil physical and chemical properties (Vieira et al. 2007; Zheng et al. 2009) and biological properties (Fromm et al. 1993). Similar to other soil properties, SOC concentration and stock are highly variable over space (horizontally and vertically) and time (Wigginton et al. 2000). Soil organic carbon density can vary due to intrinsic or extrinsic factors. Intrinsic variability is caused by natural variation (Cambardella et al. 1994) and extrinsic variability results from external factors imposed on sites (Rao and Wagenet 1985). Spatial correlation of SOC can be explored by geostatistical analysis. On a national scale, Lacelle established the soil landscape and SOC description database of Canada, and calculated that SOC storage at depths of 0-30 cm was 72.8×1012kg (Lacelle 1998). Bernoux estimated 0-30 cm SOC pool of 36.4(±3.4)×1012kg in Brazil with soil and vegetation maps (Bernoux et al. 2002). Soil carbon reserves in different soil horizons were studied in New Zealand (Scott et al. 2002) and Congo (Schwartz and Namri 2002). In China, with the aid of GIS spatial analysis, soil organic carbon density (SOCD) and storage were estimated for the entire country and a high degree of spatial variability was documented. Many studies investigated the spatial variability of SOC on a regional scale. Xue et al. (2003) surveyed the Qaidam basin and calculated the SOCD to be 6.11 kg·m-2. Using current 1:1000000 digital soil maps and 736 typical soil profiles in northeast China, organic carbon storage in topsoil and its spatial distribution were estimated, and the average SOCD at 0-1 m depths was 16.13 kg·m-2(Sun et al. 2004). In the tropical area of south China, SOC storage and spatial distribution were studied by combining soil data with a digital elevation model (DEM). Similar studies were conducted in Jiangsu, Guangdong and Anhui provinces (Jiang et al. 2005; Gan et al. 2003; Cheng and Xie 2009). According to previous studies, the initial application of GIS was to soil and land resource survey and evaluation, with focus on assessment of farmland soil degradation and mapping. Forestland, as a high potential land use type for carbon storage, has not yet gained much attention. Moreover, forest SOC storage or SOCD estimation on regional scales has large uncertainties because it is mainly based on soil survey data and literature and lacks detailed and reliable measurement. Ordinary kriging (OK) has been widely used for spatial interpolation in plains, but its predictive capability is limited in regions of complex topography with highly developed spatial conditions such as mountains and hills. Some ancillary variables such as terrain attributes and vegetation indices extracted from DEMs and remote sensing images are available for digital soil maps. Combining maps of soil types with soil properties information from point observations can improve spatial prediction of soil properties compared with prediction from a soil map only (Goovaerts and Journel 1995).

    The forest area in northeast China covers 3590×106ha, accounting for 24% of the total forest area of China, and the total forest stock volume amounts to 3.2×109m3(Department of Forest Resources Management 2010). It is the largest forest area and the transition region between broad-leaved mixed forest to the south and pure coniferous forest to the north. Research in this region has mainly focused on the relationship between forest and climate change and on spatial heterogeneity of organic carbon storage in vegetation layer. However, the forest soil carbon pool has seldom been studied and most spatial interpolation used OK. Information on the spatial distribution of SOC is important for evaluating soil ecological functions and understanding soil carbon sequestration processes. The use of geostatistics in combination with DEM and regression-kriging (RK) can provide such information and help to define and improve the accuracy of representations of the distribution of soil parameters across landscapes. Increased accuracy in the delimitation of SOC content classes could assist in meaningful soil quality prediction and assessment of the impacts of soil redistribution. The objectives of the present study were to evaluate the carbon-related soil properties and terrain attributes, and to determine the spatial variability and SOC distributions in different soil horizons of Humic cambisols in forests of northeast China using geostatistical analysis. This information will provide a scientific basis for sustainable forest management and ecological environment improvement in this region.

    Materials and methods

    Study area

    The study was carried out in Jincang Forest Farm, a typical forest area in the north temperate region of Wangqing County, Jilin Province, China (Fig. 1). The study area was in the eastern part of the low-middle hill region of the Changbai Mountain range (43°19′-43°23′ N, 130°26′-130°37′ E) at elevations of 562 to 950 m. The study area has a north temperate, monsoon climate with hot summers, cold winters, and abrupt changes of temperature. Average annual air temperature is 3.9°C. Average annual precipitation is 547 mm, concentrated from July to September (1970–2010). Soils are characterized by thick and black topsoil, loamy texture, an acidic pH range, large amount of humus and high natural fertility. The predominant soils in the area are Humic cambisols, together with meadow, peaty, swamp soil in some valleys that are developed from residual deposits of parent material such as granite, gneiss and basalt. Due to special landforms, abundant precipitation and changeable climate, rich vegetation resources are formed, including coniferous forest, broad-leaved forest, and mixed coniferous and broad-leaved forest. The major arbor species are Changbai larch (Larix olgensis), spruce (Picea koraiensis), Korean pine (Pinus koraiensis Sieb. et Zucc.), maple (Acer mono Maxim.), ash (Fraxinus mandshurica Rupr.), Asian white birch (Betula platyphylla Suk.) and linden (Tilia amurensis Rupr.).

    Soil sampling and analysis

    In August 2011, a total of 63 temporary circular sample plots (250 m2each plot) were set up in the study area. Within each plot, a soil profile was dug to describe profile characteristics and samples from the “bottom up” were taken with a ring sampler formeasurement of soil bulk density. The surface litter was removed prior to composite soil sampling. Soil samples at depths of 0-20, 20-40, and 40-60 cm were collected with 5 soil cores each well mixed into a composite soil sample, which was further divided into two sets of sub-samples. One set of the sub-samples was air dried and passed through a 2 mm sieve for analysis of soil texture, pH and available P and K; the remaining sub-sample was air-dried, homogenized and passed through a 0.25 mm sieve for analyses of SOC, total N, P and K (Li et al. 2011). Gravel and stones (fractions >2 mm) were sorted out and weighed. Soil texture was measured using the hydrometer method (Wilde et al. 1972). Soil pH was measured (1:2.5 = w/v) with a pH-meter (pH/ISE meter, model 710A; Orion). Soil organic carbon (SOC) concentration was determined by Walkley-Black wet oxidation method (Nelson and Sommers 1996; Bao 2000). Soil total N (TN) concentration was determined by Kjeldahl’s azotometer and total P (TP) by spectrophotometer, following wet digestion in concentrated H2SO4(Bremner 1996). Total K (TK) was determined by flame spectrophotometer following wet digestion in HF-HClO4(Jackson 1958; Knudsen et al. 1982). Available P (AP) was extracted with HCl-NH4F and estimated by spectrophotometer (Bray and Kurtz 1945). Readily available K (AK) was extracted with 1.0 M NH4OAc and quantified by flame spectrophotometer (Carson 1980). Soil bulk density was determined by oven-drying the soil samples collected with a ring sampler at 105°C for more than 8 h. In addition to soil sampling, the longitudes and latitudes of the sampling points were recorded with global positioning system (GPS). Forest type and survey, including environmental information, vegetation and management, for each sampling point were also recorded. A digital elevation model (DEM) of the study area was acquired at 5-m spatial resolution. The DEM and sampling points were geo-referenced to the Gauss Kruger coordinate system with Pulkovo 1942 GK Zone 22 (Fig. 1).

    Fig. 1 Study area location and sampling design.

    Calculation of soil organic carbon density (SOCD)

    Soil organic carbon density (SOCD) is the organic carbon content per unit area. We calculated SOCD in the profile with k layers (Sk, kg·m-2), on a ground-area basis up to 60 cm depth as follows:

    where, Ciis SOC concentration (g·kg-1), Diis soil bulk density (g·cm-3), Eiis soil thickness (cm) and Giis volumetric percentage of the fraction >2 mm (rock fragments). Siis soil organic carbon density in the ith horizon. In the present study, k=3.

    Statistical analysis

    Descriptive statistics

    A descriptive statistical analysis was carried out to determine the mean, minimum, maximum, standard deviation (SD) and coefficient of variation (CV) of the original data variables. The correlations among SOC, other soil properties and topographic factors such as elevation and aspect were also tested using correlation analysis and multiple linear stepwise regression. The significance of the difference between means of each property variable for different soil horizons was analyzed using one-way ANOVA. The statistical software SPSS 18.0 was used for statistical analysis.

    Geostatistical analysis

    The terrain analysis was based on a 5-m resolution digital elevation model (DEM) generated from a 1: 5000 contour topographical map. The elevation, slope in degrees, and aspect expressed in positive degrees from 0 to 360 measured clockwise from north of the study site were generated from the DEM using ArcGIS 9.3 (Esri Company, RedLands, USA). Geostatistical analysis was performed using the software GS+ and ArcGIS 9.3. A sample semivariogram was calculated for each soil variable as follows (Goovaerts 1999; Wang and Shao 2013):

    where, z(xi) and z(xi+h) are sample values at two points separated by the distance interval h, and N(h) is the total number of sample pairs for the lag interval h. For irregular sampling, it is rare for the distance between the sample pairs to be exactly equal to h. Therefore, h is often represented by a distance interval. γ(h) is the sample semivariance for the distance lag h, which reflects spatial relationship between neighbor sampling points.

    The sample semivariogram was calculated and then fitted with a suitable theoretical model chosen from spherical, exponential, Gaussian, Mate′rn or Stein’s Mate′rn (Cressie 1993). The semivariogram models provide information about the spatial structure as well as the input parameters for kriging interpolation. Spatial dependence of soil variables can be evaluated by three important parameters: nugget (C0), partial sill (C1) and range in the semivariogram plot. The nugget is defined as the variability at a scale smaller than the sampling interval and/or sampling and analyticalerror; the partial sill reflects the amount of spatial structural variance; and the range expresses the distance at which the semivariogram stabilizes around a limiting value. The spatial variation distance increasing to a constant value (sill) is known as range of spatial dependence (Cambardella et al. 1994), and it is affected by the sampling scale. The distribution distances of sampling points closer to the range are spatially correlated and the distances greater than the range are independent. In addition, the fractal dimension could be used to determine the spatial variability of soil properties. When the fractal dimension values are small, soil properties have weak spatial dependence and simple distribution, or vice versa (Burrough 1981, 1983; Wang and Shao 2013).

    Regression-kriging (RK) uses traditional regression methods that model the trend or drift in conjunction with geostatistical methods that model the random residuals from the regression (Hengl et al. 2004; Rivero et al. 2007). In RK the trend and residuals from the regression are fitted separately and then summed afterwards. It considers influence of geographic factors in the complex topography and simulates both spatial distribution trend and uncertainty. The trend item was separated by establishing regression equations between the auxiliary variable and target variable, and interpolation was applied to the residual with ordinary kriging (OK). The predicted value of the target variable was then obtained by spatial overlay between the trend item of regression prediction and OK value of the residual.

    Accuracy evaluation

    The prediction accuracy could be tested by comparing observed values with predicted values for soil nutrients, soil organic carbon (SOC) and soil organic carbon density (SOCD) from sample plots. Mean prediction error (M) and root mean square error (R) of prediction were adopted as the evaluation indices, and R was used to analyze the relative accuracy improvement value by comparing regression-kriging (RK) with ordinary kriging (OK) (Lal 2000; Kucharik et al. 2001).

    where, n is a sample point of test set, Zoiand Zpiare, respectively, observed values and predicted values, ROKis root mean square error of prediction using OK, and RRKis the prediction using RK. The prediction accuracy of RK is higher than OK when the R1value is positive, and the larger the value, the greater the accuracy. In this study, 10 soil sampling points were randomly extracted from 63 points to test the prediction accuracy of the remaining 53 points.

    Results

    Descriptive statistics of soil organic carbon (SOC) and other soil properties

    Table 1 shows the descriptive statistics of soil organic carbon (SOC) concentration and soil organic carbon density (SOCD). The mean of soil characteristics indicated the central tendency, and standard deviation (SD), coefficient of variation (CV) and skewness were used as the estimates of variability from each site. The average SOC concentration was 42.32 g·kg-1in the topsoil, and it decreased with depth until reaching 14.35 g·kg-1at 40-60 cm. The CV was low for SOC concentration in 0-20 cm (CV<10%), indicating low variability of SOC in topsoil. However, at soil depths of 20-40cm and 40-60 cm, SOC concentration showed high CVs (67% and 80%). Accordingly, the SOCD was 7.98 kg·m-2for the topsoil, and it declined gradually with increasing depth to 5.23 kg·m-2at 20-40 cm and 3.06 kg·m-2at 40-60 cm. Although there was a decreasing trend with increasing depth, the subsoil also contained a considerable amount of C, and the total amount of organic C stored to the depth of 60 cm was 16.23 kg·m-2, with a CV of 52%.

    Table 1: Descriptive statistics of soil organic carbon (n=63)

    The descriptive statistics for other soil properties in the study area are shown in Table 2. Soil pH had a mean of 5.41-5.57 with similar SDs of 0.26-0.30 in three soil layers. The CVs of total N (TN), total P (TP) and available P (AP) all exceeded 50%, which indicated considerable spatial variation in these properties within the study area. The concentrations of TN, TP and AP had means of 1.07 g·kg-1, 0.38 g·kg-1and 25.49 mg·kg-1, respectively, at 0-60 cm depth for all sampling points, and means and SDs varied by soil horizon. Mean values for all of these properties peaked at 0-20 cm and declined with increasing depth.

    Soil properties other than AP varied significantly by soil horizon. Concentrations of total N (TN), total P (TP), total K (TK) and readily available K (AK) declined with increasing depth in the following rank order: 0-20 cm > 20-40 cm > 40-60 cm withsignificant differences between each horizon (p<0.05); this trend was, however, not followed by available P (AP). In general, the typical distribution of soil chemical properties suggested higher variation between horizons and better nutrient status in the topsoil.

    Table 2: Descriptive statistics of soil properties in different soil horizons

    Correlation analysis The correlation matrix of selected soil properties was shown in Fig. 2. Soil organic carbon (SOC) concentration was significantly and positively correlated with total N (TN), total P (TP), total K (TK) and readily available K (AK), but not with available P (AP). Soil TN was positively related with SOC concentration (r = 0.795, p<0.01), but its correlation with AP was weak. Correlations of AP with TN and TP were weak, while correlation of AP with AK became stronger with increasing depth. In each horizon, there were strong positive correlations between soil properties. In addition, TN, TP, TK and AK had significant positive correlations with each other. The difference in SOC concentration was caused by natural organic matter distribution. To a certain extent, soil nitrogen levels affect SOC concentration through soil microorganisms. Different concentrations of P and K can be caused by different effects of organic matter, different pH values, and mineral composition.

    To the depth of 60 cm, the significant correlations between soil organic carbon (SOC) concentration and topographic attributes indicated that SOC concentration tended to be greater at higher elevation. The correlation between SOC concentration and elevation was very significant and positive (r =0.532, p<0.01). A very significant positive correlation was also found between soil organic carbon density (SOCD) and elevation (r = 0.565, p<0.01), but there was a very significant negative correlation between SOCD and slope (r = -0.451, p<0.01).

    Fig. 2: Correlation matrix between soil properties.??significant at probability level of 0.01;?significant at probability level of 0.05. TN is total N; TP is total P; TK is total K; AP is available P; AK is readily available K; SOC is soil organic carbon.

    There were very significant correlations between SOC concentration and SOCD (Table 3), and between these variables and cosine of slope aspect (p<0.01). SOC had, however, weak correlation with sine of slope aspect (which showed the degree of facing north), indicating that the more the aspect faced north the larger the SOC concentration. It tended to be greater on shady slopes at higher elevations. SOC in the topsoil had the same correlations with elevation and cosine of aspect, but this was not the case at depths of 20-40 cm and 40-60 cm.

    Table 3: Correlations between SOC concentration as well as SOCD and spatial factors (n=63)

    level of 0.05. SOC is soil organic carbon; SOCD is soil organic carbon density. cosα=cosine of aspect; sinα=sine of aspect. Aspect represents the degree of facing north.

    Multiple linear stepwise regression of soil organic carbon (SOC) concentration and soil organic carbon density (SOCD) on topographic factors yielded the following equations. Elevation and cosine of aspect (cosα) emerged as the optimal predictors of SOC and SOCD.

    Equations for 0-60 cm were:

    where, SOCis soil organic carbon concentration, SOCDis soil organic carbon density, and H is elevation. The coefficients of determination for fitting SOC and SOCD were 0.344 and 0.385, respectively, p<0.01 in both cases indicating good fits for both equations (6) and (7).

    Equations for 0-20 cm were:

    Geostatistical analysis

    Spatial variability of soil organic carbon (SOC) and semivariogram models

    Geostatistical parameters of soil organic carbon (SOC) are shown in Table 4. After interpolating the residual with ordinary kriging (OK), the log-transformed data of SOC concentration at 0-20 cm depth was fitted well by a Gaussian model having a nugget/sill ratio (NSR) of 5.7%, a range of 0.294 km, and very strong spatial autocorrelation. The large spatial range indicated a strongly structured regional pattern of SOC. The Gaussian model was best fitted to sample semivariograms of SOC at 0-60 cm depth with moderate structured dependence (NSR=47.8%). As for soil organic carbon density (SOCD), the NSRat 0-20 cm was 15.4%, suggesting a considerable development of the spatial structure and high spatial correlativity. SOC at 0-60 cm depth was fitted by an exponential model and the spatial structure was moderate to well developed (NSR=49.3%).

    Table 4: Semivariogram models for SOC concentration and SOCD’s residual error

    Accuracy evaluation for different prediction methods

    The regression-kriging (RK) method considers environmental factors and separates residual error to increase prediction accuracy. Fifty-three sampling points were randomly extracted from 63 points to conduct ordinary kriging (OK) interpolation for regression prediction residual error of soil organic carbon (SOC). Fig. 3 shows the spatial distribution of SOC concentration and soil organic carbon density (SOCD) interpolated using RK. At 0-60 cm, the spatial distribution of SOC mainly concentrated on northerly aspects. Taking the digital elevation model (DEM) in comparison, SOC concentration and SOCD were relatively rich at high elevations. Soil organic carbon at 0-20 cm showed that spatial distribution of SOC concentration had a number of high-value centers but was mainly concentrated in high elevation regions. The results suggest that slopes of northerly aspect served as sinks of organic matter and became the largest sources of SOC.

    Fig. 3: Spatial distributions of soil organic carbon (SOC) and soil organic carbon density (SOCD) interpolated by regression-kriging.

    The remaining 10 sampling points were used to test the accuracy of the two prediction methods. Regression-kriging (RK) prediction was much more detailed on spatial variation and more closely approximated spatial distribution condition of the actual value (Table 5). Positive R1values indicated higher accuracy using the RK method.

    Table 5: Comparison of accuracy using different prediction methods

    The fractal dimension values of soil organic carbon (SOC) in different soil horizons are shown in Table 6. The greatest isotropic fractal dimension was found at 0-20 cm, indicating the most complex spatial pattern of SOC and the highest degree of uniformity. The fractal dimension at 20-40 cm was the smallest, indicating a relatively simple spatial pattern and stronger spatial dependence than other soil horizons.

    In terms of slope aspect, the greatest fractal dimension of soil organic carbon (SOC) at 0-20 cm was found in the northwest-southeast (NW-SE) direction, and the smallest in the northeast-southwest (NE-SW) direction, which had the greatest spatial heterogeneity and the strongest spatial dependence, and was also the direction of the longest transect through the study site (Fig. 3). As for the coefficient of determination, NE-SW direction showed a preferential pattern of SOC distribution. At 20-40 cm, the NE-SW direction had the smallest fractal dimension with the strongest spatial dependence; east-west (E-W) direction was the preferential pattern distribution for SOC. The fractal dimension of SOC at 40-60 cm depth was greatest in the NW-SE direction; and smallest in the NE-SW direction, which had greatest spatial heterogeneity. The coefficient of determination showed that NE-SW had the preferential pattern of SOC distribution.

    Table 6: Fractal dimension of SOC in different soil layers

    Discussion

    The average soil organic carbon (SOC) concentration was 27.40 g·kg-1and soil organic carbon density (SOCD) to the depth of 60 cm was 16.23 kg·m-2, which is in agreement with previous studies that estimated SOCD in a range of 16.55-27.02 kg·m-2in this region (Zu et al. 2011; Sun et al. 2004). Because of differences inthe locations of study areas and management activities, SOC varies spatially. In addition, research methods are not standardized, and estimates are not free of uncertainty due to sampling depth and quantity.

    Soil organic carbon (SOC) and soil organic carbon density (SOCD) were higher in topsoil and declined with increasing depth. Distribution of SOC is related to vegetation types and succession, and human disturbance, so organic carbon in topsoil is greatly influenced by these factors. In the forest ecosystem, soil organic matter is mainly from litter and root decomposition. That surface soil can obtain more forest litter leads to the enrichment of SOC. The SOCD is associated with soil bulk density and stoniness (diameters >2 mm). Several studies found that SOCD increased with increasing soil depth (Yang and Wang 2005). However, some studies in forests suggest that the changes in SOC stocks are due to the development of vegetation restoration and succession. For instance, vegetation restoration can increase non-active organic carbon concentration in surface soil (Su et al. 2005). Succession of restored vegetation can also increase SOC and nitrogen concentration (Christian et al. 2001).

    Correlation analysis quantifies the tendency of variables such as SOC, soil properties, and terrain attributes to vary in similar ways. Soil properties have indirect influence on SOC by affecting litter decomposition, microbial activity and organic carbon mineralization rates (Fayez and Esmael 2006). Soil total N (TN) concentration can promote the accumulation of biomass in understory vegetation and litter generation, thereby facilitating SOC accumulation. As a result, SOC, elevation, aspect, TN and total K (TK) were all interrelated, suggesting that the geographical characteristics and soil properties were correlated with SOC. Soil organic carbon (SOC) in all soil horizons had a positive relationship with TN, total P (TP), TK and readily available K (AK), in agreement with previous observations (Mládková et al. 2005). A significant positive correlation between SOC and elevation indicated that SOC concentration tended to be greater at higher elevation. This is inconsistent with studies reporting greater SOC concentrations on the valley floor because of the deposition of water, solutes and sediments down along the slope (Mueller and Pierce 2003; Sumfleth and Duttmann 2008; Wang et al. 2001). Our findings agree with those of Tsuia et al. (2004) and Zhang et al. (2012) who documented increasing SOC with increasing elevation.

    Spatial patterns of soil organic carbon (SOC) suggest that organic matter dynamics in the field were affected by topography. The spatial variability of SOC and soil organic carbon density (SOCD) differed by soil horizon. The spatial distribution of surface SOC demonstrated more developed spatial structure than the other two horizons, implying that topsoil might be more sensitive to environmental disturbances, whereas spatial variability of SOC at depths of 20-40 cm and 40-60 cm was greatly affected by random factors and suggested considerably weaker spatial correlativity. By combining the kriging map with the digital elevation model (DEM), the SOC concentration was related to elevation, indicating greater SOC in the higher terrain area. Similar findings with regard to the importance of elevation to SOC were reported by Geng et al (2011). However, some research yielded the opposite results (Cheng et al. 2009; Zhang et al. 2012). This study showed that predictive capability of ordinary kriging (OK) was sometimes limited in regions of complex topography affected by intense cultivation and land use change. Recently, more studies have demonstrated regression-kriging (RK) can provide secondary information to improve the prediction accuracy (Lark and Webster 2006; Liu et al. 2010). It is possible that other useful ancillary information, such as soil and forest types, can be added in using RK to produce more accurate representations of the spatial distribution of SOC (Hengl et al. 2004). The RK method with elevation as the predictor gave more details in space than the OK method in predicting spatial distributions of SOC variables. In the present study, it was difficult to use the spatial information reliably to establish correlations with the SOC variables by using RK analysis because study area spanned a small range of elevations.

    Conclusions

    Soil organic carbon (SOC) concentration and soil organic carbon density (SOCD) decreased significantly with increasing soil depth. Mean SOC concentration was 27.40 g·kg-1and SOCD in the whole soil profile (0-60 cm) was 16.23 kg·m-2. Correlations between SOC and soil chemical properties varied by soil depth. Soil organic carbon concentration showed significant positive relationships with total N (TN), total P (TP), total K (TK) and readily available K (AK), but not with available P (AP). A significant positive correlation was found between SOC and elevation, as well as between SOC and slope aspect. Soil organic carbon tended to be higher on sunny slopes. In the surface horizon, spatial variation produced by random factors was larger and had very small spatial correlation. With increasing depth, SOC was less sensitive to external influence with smaller spatial variation. The effect of terrain factors on SOC was very complex. Resulting maps showed the spatial distribution of SOC in the study area based on the data of spatial dependency, and they can be used to delineate the areas with different impact of elevation.

    Arrouays H, Pelissier P. 1994. Modeling carbon storage profiles in temperate forest humic loamy soils of France. Soil Science, 157(3): 185-192.

    Bao SD. 2000. Soil and agricultural chemistry analysis. Beijing: China Agriculture Press, p. 495. (In Chinese)

    Bernoux M, Conceicao SCM, Volkoff B. 2002. Brazil’s soil carbon stock. Soil Science, 66: 888-896.

    Bray RH, Kurtz LT. 1945. Determination of total, organic and available forms of P in soils. Soil Science, 59: 39-45.

    Bremner JM. 1996. Nitrogen—total. In: Sparks DL (ed), Methods of Soil Analysis. Madison, WI: SSSA-ASA, pp. 1085-1121.

    Burrough P. 1981. Fractal dimensions of landscapes and other environmental data. Nature, 294: 240-242.

    Burrough P. 1983. Multiscale sources of spatial variation in soil. I. The application of fractal concepts to nested levels of soil variation. European Journal of Soil Science, 34: 577-597.

    Cambardella CA, Moorman TB, Novak JM, Parkin TB, Karlen DL, Turco RF, Konopka AE. 1994. Field-scale variability of soil properties in central Iowa soils. Soil Science Society of America Journal, 58: 1501-1511.

    Carson PL. 1980. Recommended potassium test. In: Dahnke WC (ed) Recommended chemical soil test procedures for the North Central region. North Dakota Agricultural Experimental Station Bulletin 499: 17-18.

    Cheng XF, Xie Y. 2009. Spatial distribution of soil organic carbon density in Anhui Province based on GIS. Sientia Geographica Sinica, 29(4): 540-544. (In Chinese)

    Christian PG, Michael GR, Robert MH. 2001. Tree species and soil textural controls on carbon and nitrogen mineralization rates. Soil Science Society of America Journal, 65: 1272-1279.

    Cressie NAC. 1993. Statistics for Spatial Data. New York: John Wiley & Sons p. 900.

    Department of Forest Resources Management. 2010. The 7th national forest resources inventory and the status of forest resources. Forest Resource Management, 2: 1-8. (In Chinese)

    Fayez R, Esmael A. 2006. Soil microbial activity and litter turnover in native grazed and ungrazed rangelands in a semiarid ecosystem. Biology and Fertility of Soils, 43: 76-82.

    Fromm H, Winter K, Filer J. 1993. The influence of soil type and cultivation system on the spatial distributions of the soil fauna and microorganisms and their interactions. Geoderma, 60: 109-118.

    Gan HH, Wu SH, Fan XD. 2003. Reserves and spatial distribution characteristics of soil organic carbon in Guangdong Province. Chinese Journal of Applied Ecology, 14(9): 1499-1502. (In Chinese)

    Geng GP, Gao P, Lü SQ, Zhang J. 2011. Spatial distribution of soil organic matter and total nitrogen in Matiyu small watershed in hilly area of middle southern Shandong Province. Science of Soil and Water Conservation, 9(6): 99-105. (In Chinese)

    Goovaerts P, Journel AG. 1995. Integrating soil map information in modeling the spatial variation of continuous soil properties. Soil Science, 46(3): 396-414.

    Goovaerts P. 1999 Geostatistics in soil science: state-of-the-art and perspectives. Geoderma, 89: 1-45.

    Hengl T, Heuvelink GBM, Stein A. 2004. A generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma, 120: 75-93.

    Jackson ML. 1958. Soil chemical analysis. NJ, USA:Prentice Hall, Inc. Englewood Cliffs, p. 498.

    Jackson RB, Caldwell MM. 1993. The scale of nutrient heterogeneity around individual plants and its quantification with geostatistics. Ecology, 74(2): 612-614.

    Jiang XS, Pan JJ, Li XL. 2005. Organic carbon density and storage in topsoil of Jiangsu and spatial analysis. Chinese Journal of Soil Science, 36(8): 501-503. (In Chinese)

    Knudsen D, Peterson GA, Pratt P. 1982. Lithium, sodium, and potassium. In: Page AL (eds), Methods of Soil Analysis, Vol. 2. Madison, Wis:American Society of Agronomy, Soil Science Society of America, pp. 225-246.

    Kucharik CJ, Brye KR, Norman JM, Foley JA, Gower ST, Bundy LG. 2001. Measurements and modeling of carbon and nitrogen cycling in agroecosystems of southern Wisconsin: potential for SOC sequestration during the next 50 years. Ecosystems, 4: 237-258.

    Lacelle B. 1998. Canada’s soil organic carbon database. In: Lal R, Kimble JM, Follett RF, Stewart BA (eds), Advances in Soil Science: soil processes and the carbon cycle. CRC Press, pp. 93-101.

    Lal R. 2000. Carbon sequestration in dry lands. Annual of Arid Zone, 39(1): 1-10.

    Lark RM, Webster R. 2006. Geostatistical mapping of geomorphic variables in the presence of trend. Earth Surf Process Landforms, 31: 862-874.

    Li LJ, Zeng DH, Yu ZY, Fan ZP, Yang D, Liu YX. 2011. Impact of litter quality and soil nutrient availability on leaf decomposition rate in a semi-arid grassland of Northeast China. Journal of Arid Environments, 75:787-792.

    Liu DW, Wang ZM, Zhang B. 2006. Spatial distribution of soil organic carbon and analysis of related factors in croplands of the black soil region, Northeast China. Agricultural Ecosystem Environment, 113: 73-81. (In Chinese)

    Liu SL, Li Y, Wu JS. 2010. Spatial variability of soil microbial biomass carbon, nitrogen and phosphorus in a hilly red soil landscape in subtropical China. Soil Science and Plant Nutrition, 56:693-704. (In Chinese)

    Lützow Mv, K?gel-Knabner I, Ekschmitt K, Matzner E, Guggenberger G, Marschner B, Flessa H. 2006. Stabilization of organic matter in temperate soils: mechanisms and their relevance under different soil conditions-a review. European Journal of Soil Science, 57(4):426-445.

    Mládková L, Bor?vka L, Drábek O. 2005. Soil properties and toxic aluminium forms in acid forest soils as influenced by the type of vegetation cover. Soil Science and Plant Nutrition, 51:741-744.

    Mueller TG, Pierce FJ. 2003. Soil carbon maps: enhancing spatial estimates with simple terrain attributes at multiple scales. Soil Science Society of America Journal, 67: 258-267.

    Nelson DW, Sommers LE. 1996. Total carbon, organic carbon, and organic matter. In: Sparks DL, Page AL, etc. (eds), Methods of Soil Analysis, Part 3. Chemical methods. Wisconsin, WI, USA: Soil Science Society of America Book Series, Vol. 5. pp. 961-1010.

    Rao PSC, Wagenet RJ. 1985. Spatial variability of pesticides in field soils: methods for data analysis and consequences. Weed Science, 33(suppl. 2):18-24.

    Rivero RG, Grunwald S, Bruland GL. 2007. Incorporation of spectral data into multivariate geostatistical models to map soil phosphorous variability in a Florida wetland. Geoderma, 140: 428-433.

    Schwartz D, Namri M. 2002. Mapping the total organic carbon in the soils of the Congo. Global and Planetary Change, 33: 77-93.

    Scott NA, Tare KR, Giltrap DJ. 2002. Monitoring land-use effects on soil carbon in New Zealand: quantifying baseline soil carbon stocks. Environmental Pollution, 116: 167-186.

    Stocking MA. 2003. Tropical soils and food security: the next 50 years. Science, 302(5649): 1356-1359.

    Su J, Zhao SW, Ma JD. 2005. Influence of man-made vegetation on carbon pool in southern Ningxia Region in Loess Plateau. Research of Soil and Water Conservation, 12(3): 50-52. (In Chinese)

    Sumfleth K, Duttmann R. 2008. Prediction of soil property distribution in paddy soil landscapes using terrain data and satellite information as indicators. Ecological Indicators, 8: 485-501.

    Sun WX, Shi XZ, Yu DS. 2004. Estimation of soil organic carbon density and storage of northeast China. Acta Pedologica Sinica, 41(2):298-300. (In Chinese)

    Tsuia C, Chena Z, Hsieh C. 2004. Relationships between soil properties and slope position in a lowland rain forest of southern Taiwan. Gerderma, 123: 131-142.

    Vieira VA, Mello CR, Lima JM. 2007. Spatial variability of soil physical attributes in small watershed. Ciência e Agrotecnologia, 31(5): 1477-1485. (In Portuguese).

    Wang J, Fu BJ, Qiu Y. 2001. Soil nutrients in relation to land use and landscape position in the semi-arid small catchment on the loess plateau in China. Journal of Arid Environments, 48(4): 537-550

    Wang YQ, Shao MA. 2013. Spatial variability of soil physical properties in a region of the Loess Plateau of PR China subject to wind and water erosion. Land Degradation and Development, 24(3): 296-304.

    Wigginton JD, Lockaby BG, Trettin CC, Nelson EA, Kolka RK, Wisniewski J. 2000. Soil organic matter formation and sequestration across a forested floodplain chronosequence. In: Proceedings of a workshop held at Clemson University, 12-14 April 1999, pp. 141-155.

    Wilde SA, Voigt GK, Iyer JC. 1972. Soil and plant analysis for tree culture. New Delhi, India: Oxford and IBH Publishing Co., p. 209.

    Xue L, Ma HZ, Cao GC, Sha ZJ, Ou LY. 2003. The application of geographic information systems on the calculation of storage of organic carbon in soils. Ecology & Environment, 12(4):419-422. (In Chinese)

    Yang JY, Wang CK. 2005. Soil carbon storage and flux of temperate forest ecosystems in northeastern China. Acta Ecologica Sinica, 25(11): 2875-2882. (In Chinese)

    Zhang W, Wang KL, Chen HS, He XY, Zhang JG. 2012. Ancillary information improves kriging on soil organic carbon data for a typical karst peak cluster depression landscape. Journal of the Science of Food and Agriculture, 92(5): 1094-1102.

    Zhao J, Zhang JM, Meng K. 2004. Spatial heterogeneity of soil nutrients in black soil, China-a case study at Hailun County. Bulletin of Soil and Water Conservation, 24(6): 53-57. (In Chinese)

    Zhao X, Wang Q, Kakubari Y. 2009. Stand-scale spatial patterns of soil microbial biomass in natural cold-temperate beech forests along an elevation gradient. Soil Biology and Biochemistry, 41: 1466-1474.

    Zheng Z, Zhang FR., Ma FY. 2009. Spatiotemporal changes in soil salinity in a drip irrigated field. Geoderma, 149: 243-248.

    Zu YG, Li R, Wang WJ. 2011. Soil organic and inorganic carbon contents in relation to soil physicochemical properties in northeastern China. Acta Ecologica Sinica, 31(18): 5207-5216. (In Chinese)

    DOI 10.1007/s11676-014-0533-3

    Project funding: The study was jointly supported by Natural Science Foundation of China (No. 31270697), the Fundamental Research Funds for the Central Universities (TD2011-2), State Forestry Administrative public service sector project “Key management techniques for the health of typical forest types in China” (20100400201) and National ‘973’project “Soil carbon stock and its temporal and spatial distribution pattern in natural forests” (2011CB403201).

    The online version is available at http:// www.springerlink.com

    Ling Liu, Haiyan Wang (), Wei Dai, Xiaojuan Yang, Xu Li

    College of Forestry, Beijing Forestry University, Beijing 100083, P. R. China. E-mail: haiyanwang72@aliyun.com

    Xiangdong Lei

    Institute of Forest Resource Information Techniques, Chinese Academy of Forestry, Beijing 100091, P. R. China.

    Corresponding editor: Zhu Hong

    97超视频在线观看视频| 久久精品国产亚洲网站| 99在线人妻在线中文字幕| 极品教师在线免费播放| 国产精品久久视频播放| 韩国av在线不卡| 欧美极品一区二区三区四区| 看黄色毛片网站| 日韩精品中文字幕看吧| 国产三级中文精品| 色综合亚洲欧美另类图片| 欧美日韩中文字幕国产精品一区二区三区| 最近最新免费中文字幕在线| 国模一区二区三区四区视频| 日本爱情动作片www.在线观看 | 少妇高潮的动态图| 成人亚洲精品av一区二区| 亚洲精华国产精华精| 一个人免费在线观看电影| 人妻制服诱惑在线中文字幕| 亚洲成人免费电影在线观看| 黄片wwwwww| 国产亚洲91精品色在线| 色吧在线观看| 亚洲一级一片aⅴ在线观看| ponron亚洲| 国产伦在线观看视频一区| 韩国av在线不卡| 亚洲最大成人av| 国产精品综合久久久久久久免费| 亚洲精华国产精华液的使用体验 | 亚洲国产精品成人综合色| 人人妻,人人澡人人爽秒播| 久久香蕉精品热| 亚洲,欧美,日韩| 中文字幕免费在线视频6| 蜜桃亚洲精品一区二区三区| 直男gayav资源| 国产av一区在线观看免费| 两人在一起打扑克的视频| av女优亚洲男人天堂| 国产黄色小视频在线观看| 日韩中字成人| 国产精品久久久久久久电影| 日本成人三级电影网站| 国产精品一区二区三区四区免费观看 | 亚洲精品影视一区二区三区av| 草草在线视频免费看| 大型黄色视频在线免费观看| 一级黄色大片毛片| 午夜免费男女啪啪视频观看 | 88av欧美| 婷婷六月久久综合丁香| 永久网站在线| 成人国产一区最新在线观看| 日韩精品青青久久久久久| 村上凉子中文字幕在线| 一区福利在线观看| 日韩欧美 国产精品| 波多野结衣高清无吗| 搡老熟女国产l中国老女人| 直男gayav资源| 国产精品人妻久久久久久| 亚洲第一区二区三区不卡| 最好的美女福利视频网| 国产一区二区在线观看日韩| 精品人妻视频免费看| 欧美日韩黄片免| 又粗又爽又猛毛片免费看| 我要搜黄色片| 日韩欧美在线乱码| 看黄色毛片网站| 天天一区二区日本电影三级| 免费av观看视频| 亚洲av成人精品一区久久| 久久久久久久亚洲中文字幕| 黄色欧美视频在线观看| 精品久久久久久成人av| 内地一区二区视频在线| 观看免费一级毛片| 精品久久久久久成人av| 女人被狂操c到高潮| 乱系列少妇在线播放| 亚洲精品久久国产高清桃花| 一边摸一边抽搐一进一小说| 精品一区二区三区视频在线| 久久久久性生活片| 免费看av在线观看网站| 熟女人妻精品中文字幕| 国产淫片久久久久久久久| 好男人在线观看高清免费视频| 在线观看午夜福利视频| 日韩亚洲欧美综合| 久久久国产成人免费| 久久久色成人| 久久午夜福利片| 嫩草影院入口| 99九九线精品视频在线观看视频| 久久人人精品亚洲av| 内地一区二区视频在线| 春色校园在线视频观看| 欧美性猛交╳xxx乱大交人| 变态另类丝袜制服| 久久精品影院6| 国产男人的电影天堂91| 国产精品一及| 国产黄a三级三级三级人| 国产黄a三级三级三级人| 欧美日韩中文字幕国产精品一区二区三区| 国产蜜桃级精品一区二区三区| 久久亚洲精品不卡| 久久久久久九九精品二区国产| 麻豆久久精品国产亚洲av| 又黄又爽又免费观看的视频| 国内毛片毛片毛片毛片毛片| 舔av片在线| 亚洲自偷自拍三级| 人人妻,人人澡人人爽秒播| 1000部很黄的大片| 亚洲自拍偷在线| 香蕉av资源在线| 搡女人真爽免费视频火全软件 | 色在线成人网| 欧美色欧美亚洲另类二区| 禁无遮挡网站| 麻豆国产av国片精品| 免费观看的影片在线观看| 欧美日韩中文字幕国产精品一区二区三区| 日本色播在线视频| 午夜影院日韩av| 91久久精品电影网| av天堂中文字幕网| 女同久久另类99精品国产91| 日韩在线高清观看一区二区三区 | 91久久精品国产一区二区三区| 亚洲精品成人久久久久久| 少妇丰满av| 婷婷精品国产亚洲av| 全区人妻精品视频| 国产高清不卡午夜福利| 亚洲经典国产精华液单| 国产极品精品免费视频能看的| 亚洲av二区三区四区| 午夜爱爱视频在线播放| 精品无人区乱码1区二区| 欧洲精品卡2卡3卡4卡5卡区| 欧美极品一区二区三区四区| 女的被弄到高潮叫床怎么办 | 久久久久久久久久成人| 18+在线观看网站| 淫妇啪啪啪对白视频| 日韩精品中文字幕看吧| 欧美人与善性xxx| 国产精品一区二区三区四区免费观看 | 婷婷色综合大香蕉| 国产精品,欧美在线| 丰满乱子伦码专区| 日韩欧美免费精品| 美女大奶头视频| 91久久精品电影网| 久久欧美精品欧美久久欧美| 日韩一本色道免费dvd| www.色视频.com| 免费电影在线观看免费观看| 久久婷婷人人爽人人干人人爱| 亚洲国产欧洲综合997久久,| 十八禁网站免费在线| 亚洲无线在线观看| 婷婷色综合大香蕉| 最近视频中文字幕2019在线8| 亚洲18禁久久av| 香蕉av资源在线| 精品久久久久久久久av| 一个人观看的视频www高清免费观看| 好男人在线观看高清免费视频| 三级毛片av免费| 99热这里只有是精品在线观看| 蜜桃亚洲精品一区二区三区| 午夜视频国产福利| 搞女人的毛片| 亚洲中文字幕日韩| 3wmmmm亚洲av在线观看| 久久久午夜欧美精品| 国产精品福利在线免费观看| 黄色女人牲交| 免费高清视频大片| 极品教师在线免费播放| 国产精品久久视频播放| 精品久久国产蜜桃| 哪里可以看免费的av片| 美女xxoo啪啪120秒动态图| 91在线精品国自产拍蜜月| a在线观看视频网站| 亚洲精华国产精华精| 国产精品久久视频播放| 韩国av在线不卡| 淫妇啪啪啪对白视频| 免费av不卡在线播放| 1000部很黄的大片| 黄色女人牲交| 免费看美女性在线毛片视频| www.www免费av| 精品一区二区三区视频在线| 日本欧美国产在线视频| 亚洲精品影视一区二区三区av| 午夜激情福利司机影院| 高清在线国产一区| 午夜视频国产福利| 欧美日韩综合久久久久久 | 国产精品野战在线观看| 在线播放无遮挡| 在线观看免费视频日本深夜| 久久久午夜欧美精品| 亚洲av免费在线观看| 亚洲性夜色夜夜综合| 欧美zozozo另类| 久99久视频精品免费| 亚洲av电影不卡..在线观看| 美女被艹到高潮喷水动态| 亚洲av五月六月丁香网| 精品福利观看| 夜夜看夜夜爽夜夜摸| av中文乱码字幕在线| 精品一区二区三区人妻视频| 美女免费视频网站| 男女边吃奶边做爰视频| 中亚洲国语对白在线视频| 国产视频内射| 国产精品一区二区免费欧美| 国产成人一区二区在线| 日韩人妻高清精品专区| 动漫黄色视频在线观看| 永久网站在线| 国产高清激情床上av| 亚洲无线在线观看| 国产亚洲精品久久久com| 中出人妻视频一区二区| 久9热在线精品视频| 一级av片app| 黄色欧美视频在线观看| 亚洲av二区三区四区| 亚洲成人精品中文字幕电影| 国产高清视频在线观看网站| 搡老妇女老女人老熟妇| 少妇人妻精品综合一区二区 | 女人被狂操c到高潮| 午夜久久久久精精品| 亚洲欧美精品综合久久99| 两人在一起打扑克的视频| 欧美一级a爱片免费观看看| 看免费成人av毛片| 久久人人精品亚洲av| 精品久久国产蜜桃| 国产欧美日韩精品亚洲av| 男女之事视频高清在线观看| 18禁在线播放成人免费| 禁无遮挡网站| 高清在线国产一区| 国产伦在线观看视频一区| 中文字幕人妻熟人妻熟丝袜美| 国产aⅴ精品一区二区三区波| 欧美激情国产日韩精品一区| 成人av一区二区三区在线看| 我的老师免费观看完整版| 欧美一区二区精品小视频在线| 日韩欧美三级三区| 午夜老司机福利剧场| 一进一出抽搐动态| 色尼玛亚洲综合影院| 国产欧美日韩精品一区二区| 亚洲av二区三区四区| 一a级毛片在线观看| 黄色视频,在线免费观看| 国产精华一区二区三区| 久久精品国产亚洲av天美| 91在线观看av| 男女那种视频在线观看| 看十八女毛片水多多多| 久久精品91蜜桃| 男人舔女人下体高潮全视频| 舔av片在线| 午夜爱爱视频在线播放| 又黄又爽又刺激的免费视频.| 亚洲欧美日韩东京热| 自拍偷自拍亚洲精品老妇| 有码 亚洲区| 日韩高清综合在线| 校园春色视频在线观看| 国产伦在线观看视频一区| 神马国产精品三级电影在线观看| 久久久久久久久大av| 99国产精品一区二区蜜桃av| 啪啪无遮挡十八禁网站| 一边摸一边抽搐一进一小说| 国产亚洲欧美98| 人人妻人人看人人澡| 欧美三级亚洲精品| 美女黄网站色视频| 人妻少妇偷人精品九色| 内射极品少妇av片p| 最好的美女福利视频网| 偷拍熟女少妇极品色| 久久精品国产99精品国产亚洲性色| 午夜日韩欧美国产| ponron亚洲| 欧美日韩中文字幕国产精品一区二区三区| 欧美色视频一区免费| 99久久精品一区二区三区| 一区二区三区激情视频| 国产成年人精品一区二区| 男人的好看免费观看在线视频| 又爽又黄a免费视频| 中文亚洲av片在线观看爽| 少妇的逼好多水| 久久午夜亚洲精品久久| 中文字幕久久专区| 国产午夜福利久久久久久| 亚洲av免费在线观看| 桃红色精品国产亚洲av| 中文字幕免费在线视频6| 国产亚洲精品久久久久久毛片| 免费av毛片视频| 99久久精品国产国产毛片| 欧美潮喷喷水| 99久久久亚洲精品蜜臀av| 亚洲av一区综合| 国产单亲对白刺激| 久久精品国产亚洲av涩爱 | 男插女下体视频免费在线播放| 精品久久久噜噜| 免费看a级黄色片| 午夜亚洲福利在线播放| 成人亚洲精品av一区二区| 波野结衣二区三区在线| 九九久久精品国产亚洲av麻豆| 国产真实乱freesex| 麻豆一二三区av精品| 99国产极品粉嫩在线观看| 俄罗斯特黄特色一大片| 女人被狂操c到高潮| .国产精品久久| 久久精品国产清高在天天线| 亚洲va在线va天堂va国产| 精品久久久噜噜| 国产精品综合久久久久久久免费| 色综合站精品国产| 亚洲欧美激情综合另类| 久久精品国产亚洲av涩爱 | 亚洲av不卡在线观看| 成人国产综合亚洲| 99热只有精品国产| 三级毛片av免费| netflix在线观看网站| 久久国产乱子免费精品| а√天堂www在线а√下载| 亚洲久久久久久中文字幕| 国产午夜福利久久久久久| 黄色视频,在线免费观看| 久久精品国产99精品国产亚洲性色| 欧美绝顶高潮抽搐喷水| 国产蜜桃级精品一区二区三区| 在线免费十八禁| 精品国内亚洲2022精品成人| 97碰自拍视频| 一夜夜www| 亚洲图色成人| 男女啪啪激烈高潮av片| 在线观看美女被高潮喷水网站| 精品99又大又爽又粗少妇毛片 | 久久久久久国产a免费观看| 国产成人aa在线观看| 内射极品少妇av片p| 国产精品免费一区二区三区在线| 很黄的视频免费| 少妇熟女aⅴ在线视频| 久久久久久伊人网av| 国产淫片久久久久久久久| 亚洲国产精品合色在线| 高清日韩中文字幕在线| 成人性生交大片免费视频hd| 老司机深夜福利视频在线观看| 天堂动漫精品| 亚洲国产色片| 91狼人影院| 久久久成人免费电影| 性欧美人与动物交配| 久久久精品大字幕| 村上凉子中文字幕在线| 久久亚洲真实| 最近最新中文字幕大全电影3| 12—13女人毛片做爰片一| 成人特级av手机在线观看| 又爽又黄a免费视频| av在线蜜桃| 亚洲av第一区精品v没综合| 精品人妻一区二区三区麻豆 | 亚洲欧美日韩高清在线视频| 一进一出抽搐gif免费好疼| 国产精品爽爽va在线观看网站| 国产精品一区www在线观看 | 露出奶头的视频| 啦啦啦啦在线视频资源| 小说图片视频综合网站| 淫秽高清视频在线观看| 99热这里只有精品一区| 久久精品91蜜桃| 嫩草影院精品99| 人妻丰满熟妇av一区二区三区| 久久人人精品亚洲av| 成人无遮挡网站| 久久国产精品人妻蜜桃| 2021天堂中文幕一二区在线观| 一区二区三区四区激情视频 | 国产毛片a区久久久久| 久久久久久久久久黄片| 久99久视频精品免费| 国产亚洲精品综合一区在线观看| 免费av毛片视频| 亚洲黑人精品在线| 国产成人a区在线观看| 亚洲成人久久性| 长腿黑丝高跟| 男人狂女人下面高潮的视频| 亚洲经典国产精华液单| 99九九线精品视频在线观看视频| 很黄的视频免费| av专区在线播放| 国产精品久久久久久亚洲av鲁大| 天美传媒精品一区二区| 欧美国产日韩亚洲一区| 国产高潮美女av| 天堂动漫精品| 一进一出抽搐动态| 日日夜夜操网爽| 午夜影院日韩av| 男女视频在线观看网站免费| 久久久精品大字幕| 亚洲自拍偷在线| 成人午夜高清在线视频| 欧美成人免费av一区二区三区| 国产精品亚洲一级av第二区| 国产一区二区三区在线臀色熟女| 99久国产av精品| 国产免费一级a男人的天堂| 天美传媒精品一区二区| 观看免费一级毛片| 在线观看av片永久免费下载| 日韩强制内射视频| 男人狂女人下面高潮的视频| 三级男女做爰猛烈吃奶摸视频| 成人国产麻豆网| 男插女下体视频免费在线播放| 亚洲在线观看片| 久久久午夜欧美精品| 欧美在线一区亚洲| 最近最新中文字幕大全电影3| 他把我摸到了高潮在线观看| 国产久久久一区二区三区| 黄色日韩在线| 五月玫瑰六月丁香| 特大巨黑吊av在线直播| 亚洲专区国产一区二区| 国产精品一区二区三区四区久久| 很黄的视频免费| 国产精品亚洲一级av第二区| 免费观看人在逋| 内地一区二区视频在线| 亚洲第一区二区三区不卡| 精品久久国产蜜桃| 免费黄网站久久成人精品| 国产乱人视频| 午夜精品在线福利| 美女被艹到高潮喷水动态| 成人三级黄色视频| 日韩欧美精品免费久久| 韩国av一区二区三区四区| 精品久久久久久久久av| av黄色大香蕉| 少妇人妻一区二区三区视频| 一边摸一边抽搐一进一小说| 国产精品亚洲美女久久久| 久久精品国产清高在天天线| 热99re8久久精品国产| 伦理电影大哥的女人| 丰满乱子伦码专区| 色尼玛亚洲综合影院| 国产精品久久久久久亚洲av鲁大| 国产精品久久视频播放| 搡老熟女国产l中国老女人| 亚洲国产欧洲综合997久久,| 国产色爽女视频免费观看| 狂野欧美激情性xxxx在线观看| 国产亚洲精品综合一区在线观看| 免费看美女性在线毛片视频| 亚洲人成网站在线播放欧美日韩| 春色校园在线视频观看| 1000部很黄的大片| 国产真实乱freesex| 国产高清视频在线播放一区| 久久精品国产亚洲网站| 99久久久亚洲精品蜜臀av| 综合色av麻豆| 十八禁网站免费在线| 88av欧美| 女人被狂操c到高潮| 久久久久国产精品人妻aⅴ院| ponron亚洲| 无遮挡黄片免费观看| 一个人免费在线观看电影| 免费黄网站久久成人精品| 精品一区二区三区人妻视频| 噜噜噜噜噜久久久久久91| 91精品国产九色| 亚洲性久久影院| 国产精品免费一区二区三区在线| 久久天躁狠狠躁夜夜2o2o| 国产亚洲精品综合一区在线观看| 麻豆成人午夜福利视频| 国内精品一区二区在线观看| 午夜亚洲福利在线播放| 啦啦啦观看免费观看视频高清| 久久久午夜欧美精品| 俺也久久电影网| 三级男女做爰猛烈吃奶摸视频| 99久久精品国产国产毛片| 免费观看在线日韩| 国产精品女同一区二区软件 | 欧美xxxx黑人xx丫x性爽| av黄色大香蕉| 久久久久久久久中文| 久久精品久久久久久噜噜老黄 | 国国产精品蜜臀av免费| 一级黄片播放器| 日本黄大片高清| 午夜亚洲福利在线播放| 亚洲av成人精品一区久久| 精品人妻偷拍中文字幕| 国产精品久久久久久av不卡| 人人妻人人看人人澡| 亚洲熟妇中文字幕五十中出| 成熟少妇高潮喷水视频| 无遮挡黄片免费观看| 国产真实乱freesex| 又黄又爽又免费观看的视频| 国产精品人妻久久久影院| 成人亚洲精品av一区二区| ponron亚洲| 嫁个100分男人电影在线观看| 白带黄色成豆腐渣| 一边摸一边抽搐一进一小说| 国产精品乱码一区二三区的特点| 国产欧美日韩一区二区精品| 此物有八面人人有两片| 国产69精品久久久久777片| 亚洲成人中文字幕在线播放| 色尼玛亚洲综合影院| 亚洲一级一片aⅴ在线观看| 99精品在免费线老司机午夜| 91麻豆精品激情在线观看国产| 91麻豆av在线| av天堂在线播放| 午夜福利视频1000在线观看| 日本熟妇午夜| 久久精品综合一区二区三区| 久久精品国产鲁丝片午夜精品 | 国产男人的电影天堂91| 国产大屁股一区二区在线视频| 国产精品综合久久久久久久免费| 毛片一级片免费看久久久久 | 欧美人与善性xxx| 亚洲熟妇中文字幕五十中出| 成熟少妇高潮喷水视频| 日本免费一区二区三区高清不卡| 国产成人影院久久av| 日韩中文字幕欧美一区二区| 免费人成视频x8x8入口观看| 国产精品人妻久久久久久| 日本精品一区二区三区蜜桃| 别揉我奶头 嗯啊视频| 无遮挡黄片免费观看| 99在线视频只有这里精品首页| h日本视频在线播放| 精品免费久久久久久久清纯| 国产亚洲精品av在线| 国产成年人精品一区二区| 国产精品伦人一区二区| 女同久久另类99精品国产91| 女生性感内裤真人,穿戴方法视频| 日韩欧美在线二视频| 国产精品久久久久久亚洲av鲁大| 久久6这里有精品| 免费av观看视频| 欧美精品国产亚洲| 亚洲人成网站在线播| 日韩精品青青久久久久久| 欧美zozozo另类| 一个人看视频在线观看www免费| 欧美一区二区国产精品久久精品| 久久精品国产亚洲av天美| 国产精品无大码| 中国美女看黄片| 国产精品亚洲美女久久久| 国产精品三级大全| 国产综合懂色| 一a级毛片在线观看| 97超级碰碰碰精品色视频在线观看| 亚洲精品在线观看二区| 国产成人aa在线观看| 午夜福利在线观看免费完整高清在 | 一级a爱片免费观看的视频| 亚洲欧美激情综合另类| 久久精品久久久久久噜噜老黄 | 国产综合懂色| 亚洲av第一区精品v没综合| 日韩中文字幕欧美一区二区|