The Generation of Soil Spectral Dynamic Feedback Using Landsat 8 Data for Digital Soil Mapping

: The soil spectral dynamic feedback captured from high temporal resolution remote sensing data, such as MODIS, during the soil drying process after a rainfall could assist with digital soil mapping. However, this method is ineffective in utilizing the images with a relatively high spatial resolution. There are an insufficient number of images in the soil drying process since those high spatial resolution images tend to have a low temporal resolution. This study is aimed at generating soil spectral dynamic feedback by integrating the feedback captured from the images with a high spatial resolution during the process of multiple drying after different rainfall events. The Landsat 8 data with a temporal resolution of 16 day was exemplified. Each single spectral feedback obtained from Landsat 8 was first adjusted to eliminate the impact of different rainfall magnitudes. Then, the soil spectral dynamic feedback was reorganized and generated based on the adjusted feedback. Finally, the soil spectral dynamic feedback generated based on Landsat 8 was used for mapping topsoil texture and compared with the mapping results based on the MODIS data and the fusion data of MODIS and Landsat 8. As revealed by the results, not only could the generated soil spectral dynamic feedback based on Landsat 8 data improve the details of the spatial distribution of soil texture, but it also enhances the accuracy of mapping. The mapping accuracy based on Landsat 8 data is higher than that based on the MODIS data and fusion data. The improvements of accuracy are more obvious in the areas with more complex surface conditions. This study widens the scope of application for soil spectral dynamic feedback and provides support for large ‐ scale and high ‐ precision digital soil mapping.


Introduction
The spatial distribution of soil property provides the basic information to guide crop planting, as well as the preservation and utilization of soil resources. Digital soil mapping (DSM) is an effective way to obtain soil spatial distribution information. It is defined as "computer-assisted production of digital representations of soil type or soil properties, which involve the creation and population of spatially-explicit information by the use of field and laboratory methods, coupled with spatial and non-spatial soil inference systems" [1][2][3]. McBratney et al. [4] and Arrouays [5] provided a comprehensive review of DSM and introduced numerous DSM approach. The traditional way to obtain the spatial distribution of soil properties is by field surveys or using interpolation method based on a large number of field samples, which is time-consuming and labor-consuming with huge costs. Recently, the most commonly used approach is constructing the quantitative relationships between soil and soil-forming or soil-altering factors based on statistic models to map soil properties or soil classes. However, it is difficult to apply this approach to the low-relief areas since the commonly-used soil environmental covariates such as topography and vegetation do not co-vary with soil conditions simultaneously [6][7][8][9]. In addition, the land use of low-relief areas is usually farmland and long-term human activities will weaken the relationships between soil and natural environmental covariates.
Remote sensing technology provides an effective way for digital soil mapping in low-relief area. Optical or radar images from sensors, such as Landsat TM, ASTER, Sentinel, etc. were commonly used to obtain the spatial distribution of soil surface characteristics [10][11][12][13][14][15]. For example, Bousbih et al. [16] mapped the topsoil clay content using radar and optical data from Sentinel-1 and Sentinel-2 in the Kairouan Plain of Tunisia. Zhai et al. analyzed the inversion ability of Landsat 8 and GF-1 images in soil organic matter and results show that the visible and near-infrared bands of the two remote sensing images is related to the soil organic matter content [17]. However, these studies used remote sensing image at one or a few specific times. Single images are static and represent only one moment in time. The mapping accuracy may be seriously affected by the atmospheric conditions and sensor noise [18]. In fact, the changes of soil reflectance from wet to dry are related to soil conditions. Just one or two images cannot capture such changes in soil reflectance.
A soil spectral dynamic feedback captured by a series of remote sensing data was proposed and proved to be effective in differentiating soil conditions over low relief areas [8,[19][20][21][22][23][24]. The soil spectral dynamic feedback is the changes of soil reflectance during the soil drying process after a rainfall event. The success of using soil spectral dynamic feedback for digital soil mapping relies on the high temporal resolution remote sensing data as it is necessary to employ remote sensors with high temporal resolutions to capture the dynamic changes in soil reflectance. The drying period usually lasts about 6-7 days or longer after a rainfall event and the temporal resolution of remote sensing data needs to be one day or higher. Therefore, Moderate-Resolution Imaging Spectroradiometer (MODIS) imagery with a temporal resolution of one day is taken as the main data to construct soil spectral feedback. However, the spatial resolution of remote sensing data with high temporal resolutions tends to be relatively low. The coarse cell size of the generated soil spectral dynamic feedback variables makes the predicted soil property maps less detailed and precise than required, especially in large-scale applications.
One solution to this problem is the fusion of remote sensing images with high temporal and spatial resolutions. Zeng et al. [25] applied an enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM) algorithm to process a series of MODIS data for the generation of fusion images with high spatial and temporal resolutions based on two pairs of MODIS images and Landsat images. The predicted soil texture maps using the fusion images showed a greater accuracy than those based on MODIS data. Nevertheless, the soil spectral dynamic feedback obtained using the fusion method involved only two remote sensing images with high spatial resolutions, while most of the spectral feedback information continued to be acquired from the high temporal remote sensing data. This prevents the obtained spectral feedback information from taking much advantage of remote sensing data with high spatial resolutions.
Despite the insufficient number of images with high spatial resolutions in a drying process for the generation of spectral dynamic feedback, multiple single images exist during multiple drying processes after different rainfall events at different times. It is possible to integrate those images after different rainfall events for generating soil spectral dynamic feedback with higher spatial resolutions. The key is rainfall magnitude that varies from one rainfall event to another, which renders the soil spectral feedback after different rainfalls incomparable and difficult to integrate for an effective spectral dynamic feedback. Therefore, this study aims at addressing two issues. The first one is whether an adjustment can be made to the single spectral feedback captured from the high spatial resolution image after different rainfall to integrate the adjusted feedback as effective soil spectral dynamic feedback. The other is whether the generated soil spectral dynamic feedback based on high spatial resolution images can be applied in the topsoil texture mapping and enhanced the accuracy of mapping result. The Landsat 8 data with temporal resolution of 16 day was exemplified in this paper. A case study of mapping soil sand and clay contents was conducted in the north of Xuancheng City (Anhui Province, China). Moreover, the prediction results were compared with the results based on the soil spectral dynamic feedback obtained from the MODIS data and the fusion data, respectively.

Study Area and Data
The study area is a farmland area located in Xuancheng City (Anhui Province, China). It is geographically located at 30°50′N-31°18′N, 118°40′E-119°18′E ( Figure 1). Most of the study areas have an elevation of less than 30 meters and a slope of lower than 2°. The climatic conditions in this region are mild and the mean annual precipitation is roughly 1429 mm. The precipitation is concentrated in the period between May and October annually. Paddy soil represents the main cultivated land soil in study area. A natural lake named Nanyi Lake covering an area of 189 square kilometer is located in the center of the study area. A total of fifty-six soil samples in the top of 0-20 cm were collected from the study area. By applying the laser diffraction technique, the percentage of soil sand and soil clay was chemically analyzed in the laboratory using a Mastersizer 2000 laser particle size analyzer (Malvern Instruments, Malvern, UK). The histograms of the soil sand and clay contents of the field samples are shown in Figure 2. The content of soil sand ranges from 0.3% to 60.2% with a standard deviation of 17.9%, while that of soil clay ranges from 3.2% to 16.9% with a standard deviation of 3.2%. The variance for the sand content is shown to be significantly greater compared to clay.   Due to the lack of daily interpolation climate data made available in this study area as well as the excessively low accuracy of precipitation products based on remote sensing data for use, the climate data including daily rainfall and cumulative potential evaporation was obtained based on the climate stations using an inverse distance weight interpolation method in this study. The climate stations were obtained from the China surface climate data daily data set (version 3.0) in the China Meteorological Data Sharing Service System (http://data.cma.cn/) and two climate stations were used in this study ( Figure 1). The reason for as few as two stations to be chosen is that there are only two stations located near the study area. Besides, the inverse distance weighted method was widely used to interpolate the spatial distribution of the climate data. It is suitable for the low-relief area and has no strict requirement on the number climate stations [26,27]. In consideration of this, this approach was adopted for our study. A weighted distance was used rather than an average of the two stations, which is based on the assumption that a closer station is assigned with a greater weight for interpolation.
A total of six rainfall events were chosen for this study to generate soil dynamic spectral feedback based on the Landsat data, which occurred on 24  Landsat 8 data was obtained from the USGS Global Visualization Viewer (https://glovis.usgs.gov/). Then, the data was orthographic and atmospheric adjusted using ENVI 5.3 (Exelis Visual Information Solutions, Boulder, CO, United States). The MODIS data was sourced from NASA's Earth Observing System Data and Information System (MYD09GA and MYD09GQ) (https://ladsweb.modaps.eosdis.nasa.gov/). Then, this data was re-projected using interactive data language (IDL). All of the images during the observation periods need to be with less cloud (less than 5%).

The Overview of the Soil Spectral Dynamic Feedback
Suppose that there is a rainfall input to the soil. Once the rainfall stops, the soil begins a drying process from wet to dry. This process is the feedback of the soil in response to the rainfall. The feedback pattern captured by remote sensors as changes in electromagnetic reflectance during the drying period is referred to as the soil spectral dynamic feedback [20,28]. The soil spectral dynamic feedback is dependent on the land surface conditions (such as topography and vegetation) as well as the soil conditions, assuming that the input rainfall is uniform. If the land surface conditions of two locations are very similar or the same, the spectral dynamic feedback will be largely related to their soils. Thus, the differences in soil conditions can be indicated by the differences in soil spectral dynamic feedback, which means that the feedbacks has a potential to be used for digital soil mapping. Under normal circumstances, the drying process is supposed to last 6-7 days or longer. It is possible to obtain around six or more remote sensing images in a drying period to comprise a soil spectral dynamic feedback for each pixel. Generally, there is no strict requirement on the number of the remote sensing images. However, the changes of the reflectance in the drying period cannot be captured if the number is too less. Thus, usually six or seven remote sensing images are recommended. If the drying period is long enough, more remote sensing data can be used to generate spectral dynamic feedback.
Soil spectral dynamic feedback is expressed as a 3D response surface ("evaporation-bands-reflectance") for each pixel ( Figure 4). The actual value of soil evaporation is difficult to be measured. Previous studies have demonstrated that there is a linear relationship between soil evaporation and the square root of cumulated potential evapotranspiration in the drying process [23,29,30]. Therefore, the easy-to-measure cumulated potential evapotranspiration of soil surface on the day of spectral acquisition was taken as the evaporation variable, which is indicated by the X-axis. The cumulated potential evapotranspiration is the summation of the daily potential evapotranspiration that starts on the first day to the desired day after rainfall [22]. The daily potential evapotranspiration is calculated using Penman-Monteith equation based on daily temperature, sunshine hours, water vapor pressure and other daily weather observations [31]. The Y-axis represents the bands of remote sensing images (usually is the first seven bands of the MODIS images), while the Z-axis denotes the soil reflectance [22,23].

The Generation of Soil Spectral Dynamic Feedback Using Landsat 8 Data
The temporal resolution of Landsat 8 data is 16 days. Therefore, the Landsat sensor usually can capture at most one image during the subsequent drying period to a rainfall, which means only one spectral feedback can be obtained. It is necessary to synthesize the remote sensing images obtained after multiple rainfalls to construct the soil spectral dynamic feedback. However, the magnitudes of precipitation for different rainfalls could be different. For a specific location, the different rainfall magnitude will cause soil moisture to vary, which renders the evapotranspiration that indicates the content of soil moisture during the drying process incomparable. Therefore, it is impractical to compare the feedback captured from the Landsat 8 data after different rainfall. The idea of this paper is to adjust the evapotranspiration of feedback after each rainfall event against a unified reference for them to be comparable before organizing all the adjusted feedbacks to generate soil spectral dynamic feedback.
Therefore, the generation of soil spectral dynamic feedback using Landsat 8 data includes the following three procedures ( Figure 5). The first step is to determine a unified reference for evapotranspiration variable at each pixel. The spectral dynamic feedback of MODIS band 7 after saturated rainfall was treated as a unified reference in this study. The second step is to obtain a series of evapotranspiration adjustments after different rainfall based on MODIS band 7 and construct the relationships between evapotranspiration adjustment and rainfall magnitudes. At a given pixel, we denote that the cumulative potential evapotranspiration as CET' and its corresponding spectral reflectance as r in the unified reference. For a non-saturating rainfall with a magnitude of mx at the same pixel, we denoted that the cumulative potential evapotranspiration as CET when the corresponding spectral reflectance is also r. CET' is larger than CET because the soil moisture loss during the drying period after the saturated rainfall is more than that after the non-saturating rainfall. This difference between CET' and CET will be treated as the evapotranspiration adjustment (ΔCET). The relationships between evapotranspiration adjustment and rainfall magnitudes can be constructed when a series of evapotranspiration adjustments after multiple rainfall (mn) were obtained. The final step is to adjust the evapotranspiration of the feedback captured from Landsat 8 according to the relationship with the assumption that this relationship for every single Landsat pixel located in one MODIS pixel is the same. Then the adjusted feedbacks can be reorganized to generate soil spectral dynamic feedback for each Landsat pixel. The detailed explanations for each step are provided in the following Sections 2.3.1 to 2.3.3.

Determination of a Unified Reference Based on the MODIS Data
This unified reference refers to a spectral dynamic feedback with reduced dimensionality, which relies solely on the reflectance of MODIS band 7 and the corresponding evapotranspiration value after a rainfall with a specific magnitude. MODIS band 7 is chosen for its shortwave infrared band with a spectrum ranging from 2105 nm to 2155 nm, which is most sensitive to soil moisture and is superior in reflecting the change in soil moisture during the drying process [22,32,33].
The key to determining a unified reference is to choose a specific rainfall. As we know, if rainfall magnitude increases from 0 to a certain amount, the moisture in the soil tends to be saturated and will remain unchanged due to the limited water holding capacity of soil. The daily cumulative potential evapotranspiration after each saturated rainfall can be used to indicate the state of soil moisture during the drying process, which is basically consistent for specific soils. This makes the soil spectral dynamic feedback comparable. In this case, the evapotranspiration and the corresponding MODIS 7 spectral feedback after saturated rainfall can be treated as unified reference. A rainfall threshold (20 millimeter) was set according to the results of Zeng et al (2017) [24], for which the rainfall that exceeds the threshold is referred to as saturated rainfall. The dynamic feedback following saturated rainfall can be defined as follows: where , refers to the unified reference, which is the relationship between the reflectance of MODIS band 7 ( , ) and daily cumulated potential evapotranspiration ( ′ , ) based on a statistics model (f) during the drying process after the saturated rainfall. , is calculated by summing up the daily potential evapotranspiration that starts from the first day to the day of remote sensing image acquisition after rainfall. The daily potential evapotranspiration is calculated using Penman-Monteith equation based on daily temperature, sunshine hours, water vapor pressure and other daily weather observations [31]. f can be the exponential model, which was deduced from the equations used in soil evaporation studies and the relationship between soil water content and soil surface reflectance found in previous studies [32,34]. The details of deduction have been described in Guo et al. (2015) [23]. Thus, the relationship between the soil reflectance of a shortwave infrared band and the daily cumulated potential evapotranspiration can be expressed as: where ai,j and bi,j denote the fitting coefficients at the pixel i, j.
In this study, the evapotranspiration variable of feedback will be adjusted. We need to obtain the expression formula of evapotranspiration as the reference for further adjusting. Thus, we converted equation (2) to the expression of evapotranspiration using the reflectance ( , ): where ′ , is deemed as the unified reference for evapotranspiration variable based on the saturation rainfall.

Construction of the Relationship between Evapotranspiration Adjustment and Rainfall Magnitude
In order to construct the relationship between evapotranspiration adjustment and rainfall magnitude, the first step is to calculate the evapotranspiration adjustment during the multiple drying periods under different rainfalls. Besides, the calculation of evapotranspiration adjustment requires the involvement of not only the MODIS data with a high temporal resolution, but also the reflectance of MODIS band 7. Before each rainfall, the soil is required to endure a persistent period (longer than 15 days) of low or even no precipitation to ensure that the soil is in a consistently dry state. The drying period is better to be fallow period with no vegetation or sparse vegetation, such as the winter or early spring of each year.
Generally, the spectral value will decrease as the soil moisture increases and the spectral value tends to be stabilized until the moisture reaches a certain level [23,35,36]. For a specific location, the spectral reflectance is assumed to be r when the cumulated evapotranspiration is CET'r in the reference. Assuming the cumulated evapotranspiration is CETmr when the corresponding spectral reflectance is r as well after an unsaturated rainfall m. Because the initial soil moisture after a saturating rainfall is more than a non-saturating rainfall, an equal spectral reflectance implies that the moisture in soil loss after saturated rainfall is more than that after non-saturating rainfall. The evapotranspiration CETmr is smaller than CET'r. The difference between CETmr and CET'r can be treated as evapotranspiration adjustment (ΔCET) for the rainfall m.
Therefore, for each location (i,j), the n spectral feedback of short wave infrared band ( , , , ) of MODIS can be obtained after the rainfall (mi,j), and the evapotranspiration adjustment (Δ , , , ) for each feedback can be calculated as follows: where 20 mm is taken as the rainfall threshold for this study area, and , , , is the cumulated potential evapotranspiration on the nth day during the drying process after a rainfall event with a magnitude of m for pixel (i, j), which is the summation of the daily potential evapotranspiration that starts on the first day to n day after rainfall [22]. The daily potential evapotranspiration is calculated using Penman-Monteith equation [31].
′ , , , is the cumulated potential evapotranspiration in the unified reference when the spectral feedback is , , , as well, which is calculated based on Equation (3) Once a series of evapotranspiration adjustments under different rainfalls were made for each pixel, the relationship between rainfall and evapotranspiration adjustment can be constructed by applying a statistical model for each pixel. The logarithmic model, exponential model, and polynomial regression model were used as the candidate statistical models and the best one with highest R 2 and lowest AIC will be selected. In this study, the cubic polynomial model was finally used because of its best performance in fitting the relationship between rainfall and evapotranspiration adjustment according to the findings of Zeng et al. [37]: where mi,j is the rainfall magnitude for pixel (i, j).

Organization Soil Spectral Dynamic Feedback Using Landsat 8 Images after Multiple Rainfall Events
Every single spectral feedback obtained from Landsat 8 data following each rainfall event is subject to adjustment based on the relationship between the evapotranspiration adjustment and rainfall magnitude. Then, all the adjusted feedbacks can be organized as spectral dynamic feedback for each pixel.
As the relationship between the evapotranspiration adjustment and rainfall magnitude is constructed using MODIS data, it is assumed that the relationship between the evapotranspiration adjustment and rainfall magnitude for each location is the same within a MODIS pixel. Further with this, the model can be applied to calculate the evapotranspiration adjustments according to the magnitude of rainfall for each Landsat pixel. Firstly, the spatial resolution of the coefficients between the evapotranspiration adjustment and rainfall magnitude is resampled to 30 meters, which is the same as Landsat 8 data. Secondly, the evapotranspiration adjustment after each rainfall for every pixel is calculated for the adjustment of the evapotranspiration variables in the spectral feedback captured using Landsat 8 data. If the magnitude (k) of rainfall is no less than 20 mm, there is no need to make adjustment to the evapotranspiration of the feedback. Otherwise, the evapotranspiration adjustment (ΔCET'i,j,k) could be calculated based on the rainfall magnitude (k) using the cubic polynomial model obtained in Section 2.3.2.
Then, the feedback obtained under the rainfall k is adjusted according to ΔCET'i,j,k: , , , , where Fi,j,k is the feedback at location (i, j) after rainfall k, which contains the reflectance of the seven bands (λ) of Landsat 8 and the adjusted evapotranspiration ( , , Δ ′ , , ).
, , denotes the potential evapotranspiration accumulated on a daily basis for the location on the day of Landsat image acquisition after rainfall.
A series of adjusted feedbacks can be obtained when all of the feedbacks captured by Landsat 8 data after different rainfall are adjusted. These adjusted feedbacks can be reorganized as soil spectral dynamic feedback ( , , , , , , … , , , ) by sorting the volume of evapotranspiration with ascending order for each pixel. The evapotranspiration in , , is smaller than , , , the evapotranspiration in the , , is smaller than , , , and so on. Then, an exponential algorithm (equation 2) is applied to obtain the relationship between the soil reflectance of each band and evapotranspiration to determine the changes in soil reflectance as the soil dries [23].

Methods for comparisons
In order to evaluate our proposed method, the method of fusion Landsat 8 data and MODIS data for generating soil spectral dynamic feedback variables was also applied for soil mapping and the predicted soil maps were compared [25]. The enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM) was applied to generate the fusion data with both high spatial and temporal resolutions on the desired date [38]. The entire procedure is detailed as follows: Firstly, the similar neighbouring pixels around the central pixel are identified. Based on Landsat 8 images, the difference of spectral reflectance between the central pixel and the neighbouring pixels in a local window is calculated. If the difference is less than the standard deviation of the spectral reflectance in each local window, this pixel is identified as a similar neighbouring pixel.
Secondly, each similar neighbouring pixel is assigned a weight for determining their contribution to predicting the value of central pixel of the fused data. The weight of similar neighbouring pixel is calculated depending on their spatial distance to the central pixel as well as the spectral similarity between the similar pixel and its corresponding MODIS pixel [38,39]. The neighbouring pixel with a shorter distance and a greater similarity will be assigned with a larger weight.
Finally, the conversion coefficient between the Landsat 8 and MODIS reflectance is calculated in the local window through the construction of a linear regression model with all of the similar neighbouring pixels and their weight used. The coefficients are applied to predict the fusion data on the desired date.
In addition, the soil spectral dynamic generated using the MODIS data with a high temporal resolution is also applied for soil mapping and the result was compared with the result based on the proposed method. There are three cases of the soil spectral dynamic feedback constructed based on the MODIS data in this study.

Extracting the Characteristics from Soil Spectral Dynamic Feedback
The soil spectral dynamic feedbacks for each pixel were organized as a 3D surface (Figure 4). The value of spectral reflectance changes with cumulated potential evapotranspiration and the bands of remote sensing data at the same time, which shows the fluctuation of the surface and reflects the internal structure of the spectral signals. This fluctuation and the internal structure are the important features to characterizing the soil spectral dynamic feedback. Because wavelet analysis has advantages in extracting the features of a surface in multiple directions (such as horizontal, vertical and diagonal), this study uses two-dimensional discrete wavelet analysis to extract the features of the soil spectral dynamic feedback, and takes the obtained feature variables as environmental variables [20,40,41].
The key features extracted from the dynamic feedback using wavelet analysis at each location were described by four coefficient matrices, which including the approximation coefficients matrix "cA", vertical detail coefficients matrix "cV", diagonal detail coefficients matrix "cD", and horizontal detail coefficients matrix "cH" [20]. These coefficients matrixes represented the approximation characteristics ("cA") and detailed characteristics ("cV", "cD", and "cH") from different direction of the response surface of the soil spectral dynamic feedbacks. The average and standard deviations of each coefficient matrix were calculated to represent the key characteristics of the "evaporation-bands-reflectance" surface at each location. As a result, eight variables (cAmean, cAstd, cHmean, cHstd, cVmean, cVstd, cDmean, cDstd) characterizing the spatial variation of soil spectral dynamic feedback were obtained as the environmental covariates for predicting the spatial distribution of soil texture.

Digital Soil Mapping Method
The covariates as obtained from soil spectral dynamic feedback using the proposed method and the other two methods (Section 2.4) were applied to predict the percentage of soil sand and soil clay in the study area. As one of the efficient digital soil mapping methods, the spatial prediction method based on the Third Law of Geography was applied in this study [42], which is because it has no strict requirement on the representativeness of samples and is suitable for the study area with limited samples. The basic idea of the Third Law of Geography for digital soil mapping is that similar environmental conditions will generate similar soils. The similarity in geographic configuration between each sample point and a prediction (unsampled) point is first computed. The similarity is then used as the weight in the prediction of the soil texture at the certain prediction point [42].
In this study, the covariates obtained from the soil spectral dynamic feedback were used to characterize environmental conditions. The similarities of environmental conditions , between the prediction point i and the sampled point j are calculated as follows: where coi and coj represent the values of the oth covariate at locations i and j, t indicates the number of covariates. E refers to the function to calculate the similarity between coi and coj. This similarity was calculated based on the Gower distance [43]. P refers to a minimum operator in this study and is used to integrate the similarities of the t environmental covariates into a single value. The minimum operator indicates that the lowest similarity is taken as the ultimate environmental similarity. Therefore, the environmental similarities between the predicted points i and all soil samples can be organized as a vector ⃗ : where n indicates the number of soil samples, while Si,j denotes the integrated environmental similarity for the predicted location i to soil sample j. The soil property Vi at each location i is predicted using the similarity vector and the soil property at each sampled location: where , represents the highest similarity among the environmental similarity of the location i to all the soil samples. indicates the soil property at the location m with the highest environmental similarity.
The leave-one-cross-validation was conducted to validate the accuracy of the prediction result based on the mean absolute estimation error (MAE) and root mean square error (RMSE). In addition, the predicted value was used to fit the observed value based on linear regression, and the R-square (R 2 ) was taken as an evaluation index. A higher value of R 2 suggests a higher level of accuracy.

The Extracted Characteristics of the Soil Spectral Dynamic Feedback
The characteristics extracted from the soil spectral dynamic feedback based on the Landsat 8 data, fusion data and only MODIS data are presented in Figures 6-8, respectively. In these figures, cAmean and cAstd are approximation characteristics or trend characteristics of the soil spectral dynamic feedback. cHmean and cHstd, cVmean and cVstd, cDmean and cDstd are the detailed characteristics at horizontal direction, vertical direction, diagonal direction of the soil spectral dynamic feedback, respectively. Those characteristics are used as the environmental covariates to predict soil texture.
The spatial distribution of those covariates exhibits similarity and spatial aggregation across several regions, despite some obvious differences in those maps. The spatial details of soil covariates derived from Landsat 8 data and fusion data have been significantly improved. For example, the field road and the spatial heterogeneity between fields can be displayed prominently. The mean values of approximation coefficients (cAmean) and horizontal detail coefficients (cHmean) show more effectiveness in the presentation of road details. This information is smoothed in the covariates based solely on MODIS data due to the coarse spatial resolution (250 meters).
The value of soil covariates based on fusion data shows the widest range and plenty of outliers, especially for the standard deviations of four coefficient matrix, which is significantly higher compared to those covariates based on Landsat 8 and MODIS data, for example, the range of "cAstd" varies from 0.2286 to 1.7486 with the standard deviation of 0.3466 based on the fusion data. In comparison, the range of "cAstd" varies from 0.0176 to 0.4788 with the standard deviation of 0.0319 and 0 to 0.1266 with the standard deviation of 0.0158 based on Landsat 8 and MODIS, respectively. This is because the errors are inevitably involved in the fusion process. In this case, there is a possibility that noise is caused for the fusion data. The spatial heterogeneity of each soil covariates based on MODIS data is the minimum, and the transition of spatial variation is relatively smooth due to the low spatial resolution of MODIS. The spatial variation of soil covariates obtained from Landsat 8 data is relative reasonable, which not only contain specific information, but also show gradual transition.

Compared the mapping result based on the Landsat 8 data and MODIS
Figures 9 and 10 present the predicted maps of soil sand and soil clay using the soil spectral dynamic feedback obtained from Landsat 8 data based on the proposed method and the MODIS data of case 1, respectively. The north of the study area features relative higher soil sand content than the south area and the higher content of soil clay is concentrated to the east and west of Nanyi Lake.  The details of spatial distribution of soil sand and soil clay based on Landsat 8 images show an obvious improvement. The spatial heterogeneity of soil in different fields with different crops has been reflected as well to some extent. Additionally, more detailed information, such as road, is also visible in the maps based on Landsat 8.
Regardless of soil sand or soil clay, the predicted spatial distribution maps based on MODIS data are rougher than based on Landsat 8. In the study area, farmland is divided by ridges or roads. The soil properties of each unit tend to be different due to different varieties of crops, farming measures and irrigation conditions. The area of each unit is considerably less than 250*250 square meters, which suggests that MODIS data is unsuited to the prediction of soil properties in this area. The prediction maps using Landsat 8 data with 30 m resolution can be more capable to reflect the heterogeneity of soil properties for farmland.

Result Comparison Based on the Landsat 8 Data and Fusion Data
The spatial distribution of soil sand and soil clay predicted from soil spectral dynamic feedback based on the proposed method and fusion data is shown in Figures 11 and 12, respectively. Overall, the trend of the two results is similar, with both of them being capable to reflect the spatial heterogeneity of soil. Nevertheless, there remain differences in many of the locations whether for the prediction results of soil sand or soil clay when different data is used. For soil sand, the difference is observed across these regions.  For soil clay, the difference is primarily concentrated in the east and west of Nanyi Lake, which is possibly attributed to the clay content to the east and west of Nanyi Lake being relatively high. Besides, the predicted difference in these locations based on different remote sensing images is found obvious. Table 1 shows the prediction errors of soil sand and soil clay based on soil spectral dynamic feedback constructed using different remote sensing data. The prediction accuracy using the soil spectral dynamic feedback based on the integration of soil spectral feedback of Landsat 8 after different rainfall events is relative higher than using the fusion data of MODIS and Landsat 8 data, which is because the spectral feedback information obtained from fusion data is solely based on two Landsat 8 data during the drying process after one rainfall event. The other information is acquired from the MODIS data as well in the soil drying process. For the method proposed in this study, all of the spectral feedbacks are obtained from Landsat 8 data after different rainfalls. Thus, it can provide more specific spectral information. In addition, it is unavoidable for the acquisition of fusion data from the MODIS and Landsat 8 data by applying spatiotemporal fusion algorithms to introduce error and some noise, which will affect the prediction of soil as well.

Mapping Accuracies
For soil sand content, only RMSE and MAE in case 1 based on MODIS data are slightly lower than those based on the Landsat 8. For soil clay, the RMSEs and MAEs based on MODIS data for the three cases are invariably higher compared to the RMSE (2.32%) and MAE (1.84%) based on Landsat 8 data. In addition, the R 2 based on the Landsat 8 data is generally the highest of all. In this study, it can be concluded that the prediction accuracy based on Landsat data is higher compared to that based on MODIS data or the fusion of Landsat and MODIS data.

Discussion
In order to compare the prediction error between the results based on the Landsat data and those based on the MODIS data, the difference of the absolute error of the two was examined for each sample. Firstly, the absolute prediction error of the result based on Landsat data (ErrorL) and MODIS data (ErrorM) for each sample was calculated, respectively. The result of using MODIS data is the selection of the most accurate prediction result. Then, the difference of the absolute error between ErrorL and ErrorM was calculated using ErrorL-ErrorM. A difference that is less than zero suggests that the accuracy has been improved based on Landsat for this sample. Finally, the standard deviation of the reflectance of Landsat band 7 within a MODIS pixel at each sample location was calculated. A more significant standard deviation indicates a greater heterogeneity at a specific location.
The scatter diagrams of difference between the absolute error and the standard deviation for all field samples for the soil sand and soil clay are shown in Figure 13. When the standard deviation is relatively low (< 0.035), the difference of the absolute error is irregular, as some are lower than zero but some are not. It is thus suggested that if the heterogeneous of land surface is simple, the prediction accuracy based on MODIS or Landsat does not conform to regular distribution. Due to the Landsat 8 image without cloud cover is limited, the dates of Landsat 8 data used in this paper span from November to the next February. The long time span causes instability for surface coverage, which results in noise or uncertainties. Thus, the prediction accuracy in the locations with simple land coverage is not necessarily improved based on Landsat compared to MODIS data.
However, when the standard deviation becomes larger (> 0.035), the difference of the absolute error of most samples was negative, which indicates that if the heterogeneous of the location is large, the prediction accuracy will improve based on the data with a higher spatial resolution. The location with large heterogeneity means complex land cover, which may contain farmland, vegetable fields, houses, roads and so on. The MODIS data with spatial resolution is merely 250 meter, which is inevitable to cause failure to record much of this spatial detail, thus having a negative impact on prediction accuracy. Therefore, the integration of feedback from Landsat 8 is deemed more necessary in the location with a large heterogeneity. The rainfall magnitude threshold of unified reference is set to 20 mm in this study, which is determined based on the results of Zeng et al [24] regarding the impact on the digital soil mapping of soil spectral dynamic feedback obtained after rainfalls of different magnitude. The threshold plays an important role in the proposed method and could vary from region to region. The different soil conditions will have impact on the threshold, which can also be determined according to specialist experience or knowledge about the soil moisture capacity. In addition, the threshold is not fixed. Actually, it can be slightly larger but cannot be excessively small. For example, it has little effect on the relationship between the evapotranspiration adjustment and the rainfall magnitude if the threshold was set to 25 mm in this study. The rainfall magnitude exceeding 20 mm will have little impact on the feedback. Therefore, the feedback obtained after rainfall with a magnitude ranging from 20 to 25 mm is similar to that obtained after a rainfall with a magnitude of 25 mm. The adjusted value of evapotranspiration will be similar to that as obtained when the threshold is set to 20 mm. If the rainfall threshold is set to 15 mm, it suggests that the adjustment value following the rainfall with a magnitude ranging from 15 to 20 mm is zero. However, the actual value of adjustment is higher than 0.
Generally, the soil-forming or soil-altering factors, such as climate, biology, topography, parent material, and other factors contribute to the distribution of soil. However, for our study area with relative small area and gentle terrain, the spatial distribution of climate, biology and parent material are too homogeneous to indicate the spatial variability of soils. The topography, such as elevation and slope is too gentle, which also do not co-vary with soil conditions over space. Human activities may as the important factors impact on the formation of soil in this study area. However, the spatial distribution of human activity factors, such as irrigation and fertilization, is difficult to be obtained. Under such condition, using remote sensing is an effective alternative way to predict the soil spatial distribution. The soil spectral dynamic feedback used rainfall as the input to land surface. The land surface begins a drying process from wet to dry when the rain stopped. The reflectance of different bands of the remote sensing images during the drying process can be considered as the feedback of the land surface in response to the rainfall. The soil texture is sensitive to soil moisture. Finer-textured soils, such as soils higher in clay content retain more moisture and hold it more tightly due to their finer pore-size. The change of reflectance of remote sensing data caused by the change of soil moisture with time during the drying process can indicate the variation of soil texture. Locations with different soil texture would exhibit different dynamic feedback to the rainfall if other surface conditions are the same. Thus, the soil spectral dynamic feedback can be used for digital soil mapping.
In addition, the usability of the precipitation products based on satellite, such as CMORPH and PERSIANN, was evaluated as well using observation stations in this study. The results suggested that the satellite precipitation products achieved low accuracy in this study area. Thus, the spatial distribution of daily rainfall and cumulative potential evaporation were determined by means of interpolation through climate stations, since there was no daily spatial distribution interpolation data made available. There is few interpolation methods can be chosen in this study because of the limited climate stations. The inverse distance weighted method was adopted since it has no requirement on the station number and is suitable for the low-relief area. Compared with previous studies, that one climate station is taken to represent the magnitude of an entire study area, our research is an improvement in the acquisition of meteorological data. If more climate stations available, more interpolation methods can be used to obtain the spatial distribution of climate data and a best one can be chosen. Representing the spatial variability of the climate data will cause uncertainty, which is possible to make the proposed method less effective. How the uncertainty will impact on the soil spectral dynamic feedback calculation and the effectiveness of digital soil mapping would be an important question worthy of further study.

Conclusions
The construction of soil spectral dynamic feedback for digital soil mapping in gentle areas requires remote sensing data with high temporal resolutions, such as MODIS data with temporal resolution of one day. It is difficulty to capture soil spectral dynamic feedback from Landsat 8 data due to low temporal resolution. This study is purposed to integrate the soil spectral feedback of Landsat 8 after multiple rainfall events. The major obstacle to this is to eliminate the impact of different rainfalls on the feedback. The different rainfall magnitudes will lead to the evapotranspiration, which indicates that the state of soil moisture in soil spectral dynamic feedback is incomparable. The idea of this study is to adjust the feedback captured by Landsat 8 after each rainfall to a unified reference before the integration of them. The soil spectral dynamic feedback with a higher spatial resolution was applied to digital soil mapping in the gentle area of North of Xuancheng to validate its effectiveness.
As revealed by the results, the mapping accuracy based on Landsat 8 data is higher than that based on the fusion of Landsat 8 and MODIS. It is also higher than the mapping accuracy achievable when only the MODIS data is used. Besides, the improvement of accuracy is made more obvious in the areas with more complex surface conditions. Additionally, the spatial distribution of soil sand and soil clay as predicted based on Landsat 8 data shows obvious improvement compared to that based on the MODIS data. Not only can it provide more specific information on the characteristics of land resources, it can also reflect the spatial differences of soil properties between different farmland units, which are smoothed in the prediction results based on MODIS data. Therefore, the method proposed in this study is effective in integrating the soil spectral feedback captured from Landsat 8 data in the drying period after different rainfall events. The mapping accuracy will be improved and the predicted spatial distribution can be made more specific when the integrated feedback was applied to digital soil mapping in gentle areas.