Improving Soil Moisture Estimation by Identiﬁcation of NDVI Thresholds Optimization: An Application to the Chinese Loess Plateau

: Accuracy soil moisture estimation at a relevant spatiotemporal scale is scarce but beneﬁcial for understanding ecohydrological processes and improving weather forecasting and climate models, particularly in arid and semi-arid regions like the Chinese Loess Plateau (CLP). This study proposed Criterion 2, a new method to improve relative soil moisture (RSM) estimation by identiﬁcation of normalized difference vegetation index (NDVI) thresholds optimization based on our previously proposed iteration procedure of Criterion 1. Apparent thermal inertia (ATI) and temperature vegetation dryness index (TVDI) were applied to subregional RSM retrieval for the CLP throughout 2017. Three optimal NDVI thresholds (NDVI 0 was used for computing TVDI, and both NDVI ATI and NDVI TVDI for dividing the entire CLP) were ﬁrstly identiﬁed with the best validation results (R) of subregions for 8-day periods. Then, we compared the selected optimal NDVI thresholds and estimated RSM with each criterion. Results show that NDVI thresholds were optimized to robust RSM estimation with Criterion 2, which characterized RSM variability better. The estimated RSM with Criterion 2 showed increased accuracy (maximum R of 0.82 ± 0.007 for Criterion 2 and of 0.75 ± 0.008 for Criterion 1) and spatiotemporal coverage (45 and 38 periods (8-day) of RSM maps and the total RSM area of 939.52 × 10 4 km 2 and 667.44 × 10 4 km 2 with Criterion 2 and Criterion 1, respectively) than with Criterion 1. Moreover, the additional NDVI thresholds we applied was another strategy to acquire wider coverage of RSM estimation. The improved RSM estimation with Criterion 2 could provide a basis for forecasting drought and precision irrigation management.


Introduction
Accurate and timely soil moisture (SM) information has essential applications in different fields, such as flood/drought forecasting, climate and weather modeling, water resources, and agriculture management [1]. Proper water resource management is crucial in declining vulnerability to drought and other extreme events that may occur with increasing frequency because of climate change. This has been widely recognized in the arid and

Satellite Data and Image Pre-processing
The MODIS data used in this study are composed of MODIS/Terra 500-m resolution 8-day surface reflectance (MOD09A1), and MODIS/Terra 1-km resolution 8-day LST products (MOD11A2) in 2017 for calculating ATI, NDVI, and TVDI. Moreover, the MOD09GA 500-m daily products were used to extract the acquisition time, serving for collecting the corresponding time of in situ RSM observations. We also used the 1 km for 2017 MCD12Q1. Type2 land cover product from the University of Maryland (UMD) land cover

Satellite Data and Image Pre-Processing
The MODIS data used in this study are composed of MODIS/Terra 500-m resolution 8-day surface reflectance (MOD09A1), and MODIS/Terra 1-km resolution 8-day LST products (MOD11A2) in 2017 for calculating ATI, NDVI, and TVDI. Moreover, the MOD09GA 500-m daily products were used to extract the acquisition time, serving for collecting the corresponding time of in situ RSM observations. We also used the 1 km for 2017 MCD12Q1. Type2 land cover product from the University of Maryland (UMD) land cover classification scheme (15 different land cover types over the CLP in Figure 1b) masked water bodies. Here, the mask of water bodies we used covered the water bodies' extent derived from the normalized difference water index (NDWI) for each 8-day period. The NDWI was first proposed by McFeeters in 1996 to monitor changes related to water content in water bodies [39]. To cover the CLP, five granules of MODIS data were mosaicked, re-projected, Remote Sens. 2021, 13, 589 4 of 25 and re-sampled at the resolution of 1/224 • (~500 m) using the MODIS Reprojection Tool and were clipped in ArcGIS 10.2 (ESRI Inc., Redlands, CA, USA).

In Situ Observation Data
Hourly RSM-relative soil moisture (%) of the 20-cm soil layer provided by the 213 automatic soil moisture observation stations were used in this study. The number of in situ observations used for each 8-day period in this study is given in Section 4.1. To narrow the temporal gaps between the in situ observation data and the 8-day composite products, the 8-day average in situ RSM value at each station was calculated by averaging the corresponding daily RSM. In detail, the daily granule acquisition time of the MOD09GA products (from the beginning to ending date-time) were collected first to serve as the reference for selecting corresponding in situ RSM observations. Precipitation (in mm/h) and elevation data of these stations were provided by the China Meteorological Data Service Center (http://data.cma.cn/en).

Subregional RSM Estimation
Subregional RSM estimation was applied with the two criteria. Here, we took Criterion 2 as an example ( Figure 2) and a detailed comparison of the two criteria is shown in Section 3.2.2. During the early stages of crop growth, the monitoring accuracy of ATIapparent thermal inertia (defined in Section 3.1.1) is better than that of TVDI-temperature vegetation dryness index (defined in Section 3.1.2). However, as crop growth progresses, the advantages of TVDI become evident. The average of the ATI and TVDI value, ranging from 0 to 1, was calculated, where NDVI varied from NDVI ATI and NDVI TVDI in the ATI/TVDI subregion. The idea of averaging ATI and TVDI as an assigned value in the ATI/TVDI subregion was inspired by previous studies [33,34]. A model-level integrated approach was used to effectively retrieve regional-scale daily SM. The average value of SM from the ATI-based model and the TVDI-based model was regarded as the SM when NDVI ranged from 0.10 to 0.18 [40]. Similarly, SM was obtained by averaging the PDI-based model and TVDI-based model [33].
To estimate RSM, the whole CLP was divided into three subregions (the ATI subregion, the TVDI subregion, and the ATI/TVDI subregion) according to the NDVI of individual pixels. While the ATI-based models (ATI) were merely applied to the ATI subregion (NDVI ≤ NDVI ATI ), the TVDI-based models (TVDI) were merely applied to the TVDI subregion (NDVI > NDVI TVDI ). Then, the ATI-based models and the TVDI-based models together (the ATI/TVDI joint models-the average of ATI and TVDI) were applied to the ATI/TVDI subregion (NDVI ATI < NDVI ≤ NDVI TVDI ). Here, as mentioned before, the ATIbased model is routinely used to estimate the RSM of bare soil and sparsely vegetated areas with low NDVI and the TVDI-based model is more suitable for dense vegetation coverage with high NDVI. Thus, the case of NDVI ATI greater than NDVI TVDI is not considered in an individual subregion.
The RSM estimation using the ATI-based, the TVDI-based, and the ATI/TVDI joint models should contain two key procedures. One is the iteration procedure applying three NDVI thresholds to calculate R. The other is the identification of optimal NDVI thresholds for subregional RSM estimation. In this study, we proposed a new method to select optimal NDVI thresholds for improving RSM estimation involving increased accuracy and spatiotemporal coverage, which was defined as Criterion 2. Both Criterion 1 and Criterion 2 shared similar iteration procedures (with different maximum iterations) to obtain R. and spatiotemporal coverage, which was defined as Criterion 2. Both Criterion 1 and Criterion 2 shared similar iteration procedures (with different maximum iterations) to obtain R ̅ .
In addition, the NDVI threshold (NDVI0), the lower limit of NDVI, below which the data are excluded when we derive dry/wet edges, would be used for generating TVDI [35]. To obtain optimal NDVI thresholds, NDVI0 and NDVIATI, both ranging from 0 to 0.5, and NDVITVDI, ranging from 0 to 0.7 with an interval of 0.01, were successively tested in the iterative process ( Figure 2). We carried out NDVI ranges of the iteration in the programing design used for loops statement and it was easy to run programming using the Figure 2. Flowchart of relative soil moisture (RSM) estimation by the apparent thermal inertia (ATI)-based models, the ATI/TVDI-temperature vegetation dryness index joint models, and the TVDI-based models with Criterion 2 (adapted from [35]). NDVI 0 was used for computing TVDI, and both NDVI ATI and NDVI TVDI were applied for dividing the entire Chinese Loess Plateau (CLP). The three subregions, namely the ATI subregion (NDVI ≤ NDVI ATI ), the ATI/TVDI subregion (NDVI ATI < NDVI ≤ NDVI TVDI ), and the TVDI subregion (NDVI > NDVI TVDI ), were assigned by calculated ATI, the average of ATI and TVDI, and TVDI, respectively.
In addition, the NDVI threshold (NDVI 0 ), the lower limit of NDVI, below which the data are excluded when we derive dry/wet edges, would be used for generating TVDI [35]. To obtain optimal NDVI thresholds, NDVI 0 and NDVI ATI , both ranging from 0 to 0.5, and NDVI TVDI , ranging from 0 to 0.7 with an interval of 0.01, were successively tested in the iterative process ( Figure 2). We carried out NDVI ranges of the iteration in the programing design used for loops statement and it was easy to run programming using the fixed large range for each 8-day period in 2017. For one iteration, linear relation analysis between an assigned value (e.g., ATI) and the corresponding in situ RSM observations was performed through the 10-fold cross-calibration in the calibration process. Detailed calibration and Remote Sens. 2021, 13, 589 6 of 25 validation processes are presented in Section 3.2.2. After the completion of all iterations, three groups of optimal NDVI thresholds with a maximum R in validation were selected with Criterion 2 for one 8-day period rather than one group with Criterion 1. The overall 8-day RSM map was ultimately produced with the generated subregional RSM by applying the selected optimal thresholds. The comparison between the two criteria was conducted as follows.
3.1.1. Apparent Thermal Inertia (ATI) Soil TI-thermal inertia is a thermal property of soil and describes the resistance of soil to temperature variations. TI proportionally increases as SM increases because moist soil has higher water thermal conductivity and heat capacity, thereby exhibiting a lower diurnal temperature fluctuation [27,29]. ATI-apparent thermal inertia is a simplified calculation for TI [28,41]. As one of the SM estimation methods, ATI can be easily calculated by measuring surface albedo and the diurnal temperature range as follows: where ATI is the apparent thermal inertia [K −1 ], A is the broadband surface albedo, and ∆LST corresponds to the diurnal temperature range between day and night [K]. Theoretically, the MODIS-derived ATI was computed as the ratio of the daily surface albedo and the diurnal temperature range [42]. However, the availability of the MODIS-derived ATI on some days was mainly limited by the availability of satellite LST observations because LST was often absent due to the presence of cloud during the satellite overpass times. Thus, it was difficult to obtain a continuous time series of ATI. In detail, the merged Terra 8-day surface reflectance and temperature products (MOD09A1/MOD11A2) are distributed on 8-day synthesis periods of clear sky data accumulation and each 8-day composite pixel contains the best possible observation according to specified criteria [43]. In addition, previous studies applied MOD11A2 LST products and MOD09A1 surface reflectance data (time resolution of 8 days) to Equation (1) [15,34,44]. In this case, our study was also carried out for an 8-day temporal resolution to calculate ATI. ATI varies between 0 and 1. The broadband albedo can be computed in shortwave spectral ranges from Terra MODIS surface reflectance [45]: where b 1 -b 5 and b 7 are the reflectance of bands 1, 2, 3, 4, 5, and 7 of MODIS, respectively [15,42,46,47]. MODIS's six bands are excellent in making the broadband albedo conversions under the general atmospheric conditions [45].

Temperature Vegetation Dryness Index (TVDI)
As the scatter plot of NDVI against LST forms a triangle or trapezium (hereinafter called triangle), Sandholt et al. [48] proposed the concept of TVDI, which can represent RSM, which is formulated as follows: LST min = a wet × NDVI + b wet (5) where LST represents the MODIS-derived LST in each of the pixels, LST min and LST max refer to the minimum/maximum LST in the triangle space defining the wet/dry edge at a given NDVI, respectively. a wet , a dry , b wet , and b dry are the linear regression parameters (slope and intercept) of dry/wet edges, respectively. Based on TVDI, RSM can be related to LST min and LST max with the following equation [49,50]: From Equation (6) RSM can be found as: where RSM is the relative soil moisture at any given pixel, RSM w is the maximum RSM according to wet edge, and RSM d is the minimum RSM corresponding to dry edge. The trend line of LST max gives the dry edge and that of the LST min represents the wet edge. The NDVI-LST triangle space defining the dry/wet edges on DOY-day of the year 113 over the CLP is shown in Figure 3. where LST represents the MODIS-derived LST in each of the pixels, LSTmin and LSTmax refer to the minimum/maximum LST in the triangle space defining the wet/dry edge at a given NDVI, respectively. a , a , b , and b are the linear regression parameters (slope and intercept) of dry/wet edges, respectively. Based on TVDI, RSM can be related to LSTmin and LSTmax with the following equation [49,50]: From Equation (6) RSM can be found as: where RSM is the relative soil moisture at any given pixel, RSM is the maximum RSM according to wet edge, and RSM is the minimum RSM corresponding to dry edge. The trend line of LST gives the dry edge and that of the LST represents the wet edge. The NDVI-LST triangle space defining the dry/wet edges on DOY-day of the year 113 over the CLP is shown in Figure 3.  [35,48,51]). Theoretically, in the triangular figure (pink region area), the base edge of the triangle with maximum evapotranspiration pixels, and the top edge of the triangle with zero evapotranspiration pixels are displayed. As the NDVI increases, the maximum LST decreases and can be fitted to a negative slope using the least square method, which is defined as the dry edge in red color lines (LSTmax). The base line of the triangle represents the wet edge in blue color lines (LSTmin), which is calculated by averaging a group of points in the lower limits of the scatterplot. The TVDI increases from 0 to 1 (a black arrow going from TVDI = 0 to TVDI = 1), indicating a land surface change from extreme wetness to extreme drought. Linear equations were generated when NDVI0 equals 0, 0.1, and 0.2, respectively. The linear regression coefficients (slopes, intercepts, and R 2 ) of dry/wet edges varied with different NDVI0.
Generally, the TVDI value of the dry edge, referring to the driest region of the study area, equals 1, while that of the wet edge (being the most humid region) is close to 0 [52].   [35,48,51]). Theoretically, in the triangular figure (pink region area), the base edge of the triangle with maximum evapotranspiration pixels, and the top edge of the triangle with zero evapotranspiration pixels are displayed. As the NDVI increases, the maximum LST decreases and can be fitted to a negative slope using the least square method, which is defined as the dry edge in red color lines (LST max ). The base line of the triangle represents the wet edge in blue color lines (LST min ), which is calculated by averaging a group of points in the lower limits of the scatterplot. The TVDI increases from 0 to 1 (a black arrow going from TVDI = 0 to TVDI = 1), indicating a land surface change from extreme wetness to extreme drought. Linear equations were generated when NDVI 0 equals 0, 0.1, and 0.2, respectively. The linear regression coefficients (slopes, intercepts, and R 2 ) of dry/wet edges varied with different NDVI 0 .
Generally, the TVDI value of the dry edge, referring to the driest region of the study area, equals 1, while that of the wet edge (being the most humid region) is close to 0 [52]. LST linearly changes with NDVI in the conditions of the same RSM. Between two edges (dry/wet edges), all intermediate conditions can occur, and all RSM conditions can consequently be represented within the NDVI-LST triangle space [31,53]. However, the maximum and minimum LST at lower NDVI in the scatterplot (two red circles in Figure 3) do Remote Sens. 2021, 13, 589 8 of 25 not seem to contribute to the formation of the dry/wet edges, which has been noted by previous studies [54][55][56]. Thus, in our study, the NDVI 0 is the low limit of NDVI, below which the data are excluded when we derive dry/wet edges. The TVDI (NDVI < NDVI 0 ) with Criterion 2 was calculated using the linear regression parameters of derived dry/wet edges, and the TVDI (NDVI < NDVI 0 ) was not computed because of setting NDVI 0 smaller than NDVI ATI for Criterion 1. In other words, TVDI would only be used both in the ATI/TVDI subregion and the TVDI subregion, where NDVI was higher than NDVI ATI , not to mention NDVI 0 .

RSM Estimation with Criterion 2
The overall RSM was produced with three groups of selected optimal NDVI thresholds (Criterion 2) using MODIS-derived ATI, TVDI, and the mean of ATI and TVDI against in situ RSM observations for one 8-day period. The equations were used as follows: where RSM overall represents the overall RSM and it is combined by three subregional RSM (RSM ATI , RSM ATI/TVDI , and RSM TVDI ). RSM ATI and RSM TVDI are the RSM estimated by the ATI-based and TVDI-based models, respectively, and RSM ATI/TVDI is the RSM estimated by the ATI/TVDI joint model. a ATI and b ATI are coefficients from fitting the ATI values and in situ RSM observations in the ATI subregion. a TVDI and b TVDI are coefficients from fitting the TVDI values and in situ RSM observations in the TVDI subregion. a ATI/TVDI and b ATI/TVDI are coefficients from fitting the average value of ATI and TVDI and in situ RSM observations in the ATI/TVDI subregion. NDVI ATI and NDVI TVDI are the selected optimal thresholds for generating subregions.

Calibration and Validation Processes for the Two Criteria
The two criteria differ in how the optimal NDVI thresholds are identified, namely NDVI 0 ≤ NDVI ATI ≤ NDVI TVDI for Criterion 1 and NDVI ATI < NDVI TVDI for Criterion 2. In other words, the limit condition of Criterion 1 was stricter than Criterion 2 and the value of NDVI 0 might be greater/lower than NDVI ATI for Criterion 2. As we mentioned in Section 3.1.2., dry/wet edges were calculated when NDVI was greater than NDVI 0 . In this case, the value of TVDI (NDVI < NDVI 0 ) was calculated based on the dry/wet edge derived from NDVI greater than NDVI 0 . The size relationship between NDVI 0 and NDVI ATI is the main difference of the two criteria, which directly lead to a different assigned value in the ATI/TVDI subregion, thereby affecting the calibration process when using the ATI/TVDI joint models. The detailed differences between subregional RSM estimation with each criterion are represented in Figures 4 and 5. For Criterion 2, the average of ATI and TVDI was assigned to the pixels in the ATI/TVDI subregion, in which TVDI includes two intervals (NDVI ATI < NDVI < NDVI 0 and NDVI 0 ≤ NDVI ≤ NDVI TVDI ) with Criterion 2 ( Figure 5). TVDI with Criterion 1 is from only one interval, namely NDVI ATI < NDVI ≤ NDVI TVDI in the ATI/TVDI subregion from Figure 4.
Importantly, because optimal NDVI thresholds were identified based on the R values in the validation, both criteria performed the iteration procedure (calibration and validation processes) in three subregions (the ATI subregion, the TVDI subregion, and the ATI/TVDI subregion) to calculate R for 8-day periods (procedures with a purple background in Figure 2). In order to not limit the size of NDVI 0 and NDVI ATI in the iteration procedure, more iterations would be implemented with Criterion 2 (maximum of 97,546 iterations) than that of Criterion 1 (maximum of 48,620 iterations) [35]. The number of iterations was calculated by all combinations of the three thresholds.  Figure 4. Schematic diagram of the calibration and validation processes with Criterion 1 (NDVI0 ≤ NDVIATI ≤ NDVITVDI). NDVIATI and NDVITVDI were applied to divide the whole CLP. NDVI0 was lower than NDVIATI and merely used for calculating TVDI (not shown in the figure). ATI, the average of ATI and TVDI, and TVDI were assigned to the ATI subregion (NDVI ≤ NDVIATI), the ATI/TVDI subregion (NDVIATI < NDVI ≤ NDVITVDI), and the TVDI subregion (NDVI > NDVITVDI), respectively. RSMATI and RSMTVDI were the RSM estimated by the ATI-based and TVDI-based models, respectively, and RSMATI/TVDI was the RSM estimated by the ATI/TVDI joint model. aATI, bATI, aTVDI, and bTVDI in the equations were coefficients from fitting the ATI values and TVDI values with in situ RSM observations in the ATI subregion and the TVDI subregion, respectively. aATI/TVDI and bATI/TVDI were coefficients from fitting the average value of ATI and TVDI and in situ RSM observations in the ATI/TVDI subregion. After 10 rounds 10-fold cross-calibration, R , R / , and R in validation were calculated by averaging RATI, RATI/TVDI, and RTVDI, respectively.  Figure 5. Schematic diagram of the calibration and validation processes with Criterion 2 (NDVIATI < NDVITVDI). NDVIATI and NDVITVDI were applied to divide the whole CLP. NDVI0 might be lower/greater than NDVIATI and appeared when NDVI0 was greater than NDVIATI, indicating the calculated TVDI from two NDVI intervals (blue and light orange color regions) in the ATI/TVDI subregion. The value of TVDI with Criterion 2 in the blue color regions (NDVI < NDVI0) was calculated based on the dry/wet edge derived from NDVI greater than NDVI0. The computed TVDI with Criterion 2 in the light orange color regions (NDVI ≥ NDVI0) was the same as that of Criterion 1 (orange color regions in Figure 4). For the meaning of the other parameters (RSMATI, RSMTVDI, RSMATI/TVDI, RATI, RATI/TVDI, RTVDI, R , R / , and R ), please refer to the caption of Figure 4. . Schematic diagram of the calibration and validation processes with Criterion 1 (NDVI 0 ≤ NDVI ATI ≤ NDVI TVDI ). NDVI ATI and NDVI TVDI were applied to divide the whole CLP. NDVI 0 was lower than NDVI ATI and merely used for calculating TVDI (not shown in the figure). ATI, the average of ATI and TVDI, and TVDI were assigned to the ATI subregion (NDVI ≤ NDVI ATI ), the ATI/TVDI subregion (NDVI ATI < NDVI ≤ NDVI TVDI ), and the TVDI subregion (NDVI > NDVI TVDI ), respectively. RSM ATI and RSM TVDI were the RSM estimated by the ATI-based and TVDI-based models, respectively, and RSM ATI/TVDI was the RSM estimated by the ATI/TVDI joint model. a ATI , b ATI , a TVDI , and b TVDI in the equations were coefficients from fitting the ATI values and TVDI values with in situ RSM observations in the ATI subregion and the TVDI subregion, respectively. a ATI/TVDI and b ATI/TVDI were coefficients from fitting the average value of ATI and TVDI and in situ RSM observations in the ATI/TVDI subregion. After 10 rounds 10-fold cross-calibration, R ATI , R ATI/TVDI , and R TVDI in validation were calculated by averaging R ATI , R ATI/TVDI, and R TVDI , respectively.  Figure 4. Schematic diagram of the calibration and validation processes with Criterion 1 (NDVI0 ≤ NDVIATI ≤ NDVITVDI). NDVIATI and NDVITVDI were applied to divide the whole CLP. NDVI0 was lower than NDVIATI and merely used for calculating TVDI (not shown in the figure). ATI, the average of ATI and TVDI, and TVDI were assigned to the ATI subregion (NDVI ≤ NDVIATI), the ATI/TVDI subregion (NDVIATI < NDVI ≤ NDVITVDI), and the TVDI subregion (NDVI > NDVITVDI), respectively. RSMATI and RSMTVDI were the RSM estimated by the ATI-based and TVDI-based models, respectively, and RSMATI/TVDI was the RSM estimated by the ATI/TVDI joint model. aATI, bATI, aTVDI, and bTVDI in the equations were coefficients from fitting the ATI values and TVDI values with in situ RSM observations in the ATI subregion and the TVDI subregion, respectively. aATI/TVDI and bATI/TVDI were coefficients from fitting the average value of ATI and TVDI and in situ RSM observations in the ATI/TVDI subregion. After 10 rounds 10-fold cross-calibration, R , R / , and R in validation were calculated by averaging RATI, RATI/TVDI, and RTVDI, respectively.  Figure 5. Schematic diagram of the calibration and validation processes with Criterion 2 (NDVIATI < NDVITVDI). NDVIATI and NDVITVDI were applied to divide the whole CLP. NDVI0 might be lower/greater than NDVIATI and appeared when NDVI0 was greater than NDVIATI, indicating the calculated TVDI from two NDVI intervals (blue and light orange color regions) in the ATI/TVDI subregion. The value of TVDI with Criterion 2 in the blue color regions (NDVI < NDVI0) was calculated based on the dry/wet edge derived from NDVI greater than NDVI0. The computed TVDI with Criterion 2 in the light orange color regions (NDVI ≥ NDVI0) was the same as that of Criterion 1 (orange color regions in Figure 4). For the meaning of the other parameters (RSMATI, RSMTVDI, RSMATI/TVDI, RATI, RATI/TVDI, RTVDI, R , R / , and R ), please refer to the caption of Figure 4. . Schematic diagram of the calibration and validation processes with Criterion 2 (NDVI ATI < NDVI TVDI ). NDVI ATI and NDVI TVDI were applied to divide the whole CLP. NDVI 0 might be lower/greater than NDVI ATI and appeared when NDVI 0 was greater than NDVI ATI , indicating the calculated TVDI from two NDVI intervals (blue and light orange color regions) in the ATI/TVDI subregion. The value of TVDI with Criterion 2 in the blue color regions (NDVI < NDVI 0 ) was calculated based on the dry/wet edge derived from NDVI greater than NDVI 0 . The computed TVDI with Criterion 2 in the light orange color regions (NDVI ≥ NDVI 0 ) was the same as that of Criterion 1 (orange color regions in Figure 4). For the meaning of the other parameters (RSM ATI , RSM TVDI , RSM ATI/TVDI , R ATI , R ATI/TVDI , R TVDI , R ATI , R ATI/TVDI , and R TVDI ), please refer to the caption of Figure 4.
For each subregion, to decrease the variability of the calibration, the linear fit between the assigned value (ATI, the average of ATI and TVDI, and TVDI, respectively) and observed RSM was performed using 10 rounds of 10-fold cross-calibration. To be more specific, after the random splitting of the paired assigned values and in situ RSM observations (when the number of available RSM observation stations was greater than 20) into 10 subsamples with nine as training data and one as testing data, acquired regression parameters (slopes and intercepts) in training data were applied to testing data. Then, one group of estimated RSM and in situ RSM observations was generated. After the 10th fold was accomplished, the validation data of 10 groups of estimated RSM and in situ RSM observations were computed, including R, as well as the p-value for a significance test (p-value < 0.05), which was just called "one-round". After 10 rounds of iterations, we averaged the R (R) from the 10 rounds of validation as the reference to identify corresponding optimal NDVI thresholds.

Evaluation of Estimated RSM for the Two Criteria
The threshold qualification when identifying optimal NDVI thresholds for the two criteria was different, but both were based on the averaged validation results (R), which directly represents the accuracy of RSM estimation. For an individual subregion, only the maximum R greater than 0.23 for Criterion 2 and slightly lower (0.17) for Criterion 1 should be selected; otherwise, no optimal NDVI thresholds were chosen for that subregion. For Criterion 1, the NDVI thresholds corresponding to maximum R among the three subregions were chosen as the optimal thresholds when the NDVI 0 , NDVI ATI , and NDVI TVDI thresholds existed simultaneously for one 8-day period [35]. In this way, they could not guarantee that maximum R would be reached for all the three subregions and the generated three subregions according to NDVI ATI and NDVI TVDI always excluding each other. Thus, only one group of optimal NDVI thresholds were chosen for subregional RSM estimation and then the overall RSM with Criterion 1 was combined by subregional RSM for each 8-day period [35].
However, for Criterion 2, more selected processes for optimal NDVI thresholds were carried out. The NDVI thresholds corresponding to maximum R in each subregion, namely the highest R ATI , R ATI/TVDI , and R TVDI , were chosen as the optimal thresholds. The optimal NDVI thresholds of the three subregions were not influenced by each other. Thus, three groups of optimal NDVI thresholds regarding three subregions were selected for subregional RSM estimation and then the overall RSM with Criterion 2 was combined by that subregional RSM for each 8-day period. In this case, the three subregions for one 8-day period might be overlapped due to the three groups of optimal NDVI thresholds individually selected instead of always excluding each other for Criterion 1. Based on the relationships of the selected optimal NDVI thresholds of the subregions, we list all cases in terms of different combinations of models used to produce an overall RSM map for one 8-day period (Table 1). Overlaps might exist when combining different models for RSM estimation for one 8-day period (see cases 4, 5, 6, and 7 in Table 1).
We combined the subregional RSM generated using the corresponding models with the selected optimal NDVI thresholds to obtain an overall RSM map for one 8-day period. Then, we evaluated and compared the selected optimal NDVI thresholds with each criterion in terms of the validation outcomes ( Figure 6). Furthermore, the performance of the selected optimal NDVI thresholds for RSM retrieval, involving accuracy and spatiotemporal coverage using the two criteria, was compared. Monthly, seasonal, and yearly RSM maps were generated by combining 8-day RSM maps under Criterion 2 and were eventually examined. One type model is considered to be used once in one case. 2 The subscripts a in the NDVI ATI_a and j in the NDVI ATI_j represent the selected optimal NDVI ATI thresholds in the ATI subregion and the ATI/TVDI subregion, respectively. The subscripts j in NDVI 0_j and t in NDVI 0_t mean the selected optimal NDVI 0 in the ATI/TVDI subregion and the TVDI subregion, respectively. The subscripts j in NDVI TVDI_j and t in NDVI TVDI_t refer to the selected optimal NDVI TVDI in the ATI/TVDI subregion and the TVDI subregion, respectively. 3 NDVI ATI < NDVI TVDI for an individual subregion because the ATI-based model was suitable for regions with low NVDI and the TVDI-based model was applicable for regions with high NDVI [15,57]. NDVIATI_a ≥ NDVITVDI_t Yes 1 One type model is considered to be used once in one case. 2 The subscripts a in the NDVIATI_a and j in the NDVIATI_j represent the selected optimal NDVIATI thresholds in the ATI subregion and the ATI/TVDI subregion, respectively. The subscripts j in NDVI0_j and t in NDVI0_t mean the selected optimal NDVI0 in the ATI/TVDI subregion and the TVDI subregion, respectively. The subscripts j in NDVITVDI_j and t in NDVITVDI_t refer to the selected optimal NDVITVDI in the ATI/TVDI subregion and the TVDI subregion, respectively. 3 NDVIATI < NDVITVDI for an individual subregion because the ATI-based model was suitable for regions with low NVDI and the TVDI-based model was applicable for regions with high NDVI [15,57].

Comparison of Validation Results
The measure R in validation reflects the accuracy of the estimated RSM compared to the observed RSM. Lower R represents poorer estimated RSM. Optimal NDVI thresholds were identified corresponding to the highest R in subregions. Then, RSM maps could be generated using the selected optimal NDVI thresholds for each 8-day period. Eventually, we selected 45 8-day periods (out of 46 except for DOY 73 in 2017) with Criterion 2 and only 38 8-day periods with Criterion 1 (Figure 7) [35]. Thus, 45 8-day periods of estimated RSM maps (improved 7 8-day periods compared to with Criterion 1) would be generated with Criterion 2. The R values with Criterion 2 were slightly higher than those with Criterion 1 but roughly showed a similar trend throughout the year. These  The measure R in validation reflects the accuracy of the estimated RSM compared to the observed RSM. Lower R represents poorer estimated RSM. Optimal NDVI thresholds were identified corresponding to the highest R in subregions. Then, RSM maps could be generated using the selected optimal NDVI thresholds for each 8-day period. Eventually, we selected 45 8-day periods (out of 46 except for DOY 73 in 2017) with Criterion 2 and only 38 8-day periods with Criterion 1 (Figure 7) [35]. Thus, 45 8-day periods of estimated RSM maps (improved 7 8-day periods compared to with Criterion 1) would be generated with Criterion 2. The R values with Criterion 2 were slightly higher than those with Criterion 1 but roughly showed a similar trend throughout the year. These values fluctuated irregularly with peak values (0.82 ± 0.007 with Criterion 2 and 0.75 ± 0.008 with Criterion 1, respectively) in winter [35]. In this case, the accuracy of RSM estimation with Criterion 2 was improved.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 26 with Criterion 1, respectively) in winter [35]. In this case, the accuracy of RSM estimation with Criterion 2 was improved. The applied models corresponding to the highest R varied with the 8-day periods. Table 2 displays the validation results as well as the models used for each 8-day period with Criterion 2. The seven rows in blue (DOYs 1, 9, 33, 65, 201, 225, and 241) indicate the new validation results with Criterion 2 for RSM retrieval compared to Criterion 1. The ATI/TVDI joint models were used in almost all 8-day periods (34/45 8-day periods), which was well in line with Criterion 1 [35]. The highest R (0.82 ± 0.007 on DOY 361) with Criterion 2 was higher than that of Criterion 1 (0.75 ± 0.008 on DOY 313). The validation score of R was better than previous results in the study area [35].
For DOYs 145 and 305, the ATI/TVDI joint models were applied twice in the two complementary ATI/TVDI subregions with different NDVI ranges (0.00 < NDVI ≤ 0.20 with R of 0.61 ± 0.044 and 0.25 < NVDI ≤ 0.27 with R of 0.63±0.017 on DOY 145, and 0.13 < NDVI ≤ 0.17 with R of 0.47 ± 0.021 and 0.30 < NDVI ≤ 0.34 with R of 0.44 ± 0.020 on DOY 305, respectively). In the ATI/TVDI subregions, different NDVI thresholds (NDVIATI and NDVITVDI) means the generated RSM maps with these NDVI intervals. Here, higher R , not merely the highest R , could be selected as well and their corresponding thresholds were regarded as additional thresholds for RSM retrieval. In this context, the two complementary ATI/TVDI subregional RSM could be merged to produce the overall RSM map. Thus, the more additional thresholds we selected, the more completed the obtained RSM map would be. The additional NDVI thresholds we applied, here, was another improved strategy to acquire wider coverage of RSM estimation. The number of stations used with each criterion was the same for each 8-day period. The applied models corresponding to the highest R varied with the 8-day periods. Table 2 displays the validation results as well as the models used for each 8-day period with Criterion 2. The seven rows in blue (DOYs 1, 9, 33, 65, 201, 225, and 241) indicate the new validation results with Criterion 2 for RSM retrieval compared to Criterion 1. The ATI/TVDI joint models were used in almost all 8-day periods (34/45 8-day periods), which was well in line with Criterion 1 [35]. The highest R (0.82 ± 0.007 on DOY 361) with Criterion 2 was higher than that of Criterion 1 (0.75 ± 0.008 on DOY 313). The validation score of R was better than previous results in the study area [35].
For DOYs 145 and 305, the ATI/TVDI joint models were applied twice in the two complementary ATI/TVDI subregions with different NDVI ranges (0.00 < NDVI ≤ 0.20 with R of 0.61 ± 0.044 and 0.25 < NVDI ≤ 0.27 with R of 0.63±0.017 on DOY 145, and 0.13 < NDVI ≤ 0.17 with R of 0.47 ± 0.021 and 0.30 < NDVI ≤ 0.34 with R of 0.44 ± 0.020 on DOY 305, respectively). In the ATI/TVDI subregions, different NDVI thresholds (NDVI ATI and NDVI TVDI ) means the generated RSM maps with these NDVI intervals. Here, higher R, not merely the highest R, could be selected as well and their corresponding thresholds were regarded as additional thresholds for RSM retrieval. In this context, the two complementary ATI/TVDI subregional RSM could be merged to produce the overall RSM map.
Thus, the more additional thresholds we selected, the more completed the obtained RSM map would be. The additional NDVI thresholds we applied, here, was another improved strategy to acquire wider coverage of RSM estimation. The number of stations used with each criterion was the same for each 8-day period.

The Optimal NDVI Thresholds
The thresholds NDVI ATI and NDVI TVDI were used to divide the entire study area into three subregions for subregional RSM estimation. In terms of the fixed ranges of NDVI 0 , NDVI ATI , and NDVI TVDI we used during iteration, to cover more combinations of NDVI thresholds, the variation of NDVI for all seasons was taken into consideration and maximum possible ranges were selected throughout the year. Indeed, the NDVI ranges we used were a little bit wide compared to other studies, to some extent, resulting in many redundant iterations [33,40]. There was no significant trend in these optimal thresholds throughout the year. It would be explained by taking into account seasonal change (or phenological variation) effects on prediction. To better explain the variation of the NDVI threshold, the construction of seasonal models should be considered over long periods for further studies. Moreover, the NDVI threshold could improve the application of other phenology-based RSM estimation methods meant basically for the growing season [58][59][60].
The selected minimum and maximum NDVI ATI (0.00 and 0.50) and NDVI TVDI (0.09 and 0.68) are shown in Table 2. The obvious difference between the two criteria was the identification of NDVI 0 , which not only determines the value of TVDI but also affects the value of pixels in the ATI/TVDI subregion. The pixels with NDVI ATI < NDVI ≤ NDVI TVDI were assigned the average of ATI and TVDI in the ATI/TVDI subregion for two criteria, but some of the TVDI (NDVI ATI < NDVI < NDVI 0 ) were calculated by the wet/dry edges produced from those NDVI higher than NDVI 0 for Criterion 2. The identification of NDVI ATI was not affected by the other two thresholds with Criterion 2 and the biggest NDVI ATI was 0.50 solely using the ATI-based models on DOY 201. The selected NDVI TVDI with Criterion 2 was up to 0.68, approaching the upper limited value (0.70) on DOY 81 using the ATI/TVDI joint models.
In terms of NDVI 0 for the NDVI-LST scatter plots, the value of selected NDVI 0 with Criterion 2 fluctuated over 8-day periods but remained relatively low and stable in winter (DOYs 1-49 and 337-361). Most of R 2 dry were higher than R 2 wet , particularly in summer (DOYs 153-233), and R 2 wet peaked in spring and autumn (Figure 8b). The slopes (a dry and a wet ) and intercepts (b dry and b wet ) of the regression equations changed symmetrically ( Figure 8a) and showed similar trends regardless of the criterion [35].

Evaluation of Estimated RSM at the Regional Scale
The estimated RSM area and average RSM obtained with each of the two criteria varied throughout the year (Figure 9). There were 45 8-day periods of the RSM maps produced in 2017 with Criterion 2, compared to 38 with Criterion 1 [35]. This clearly shows that using Criterion 2 improved the temporal coverage of RSM estimation.
In terms of spatial coverage, the two criteria also led to different RSM estimation results. The total area of estimated RSM with Criterion 2 was 939.52 × 10 4 km 2 in 2017, which was 40.76% higher than with Criterion 1 (667.44 × 10 4 km 2 ). Among the 45 8-day

Evaluation of Estimated RSM at the Regional Scale
The estimated RSM area and average RSM obtained with each of the two criteria varied throughout the year (Figure 9). There were 45 8-day periods of the RSM maps produced in 2017 with Criterion 2, compared to 38 with Criterion 1 [35]. This clearly shows that using Criterion 2 improved the temporal coverage of RSM estimation.
In terms of spatial coverage, the two criteria also led to different RSM estimation results. The total area of estimated RSM with Criterion 2 was 939.52 × 10 4 km 2 in 2017, which was 40.76% higher than with Criterion 1 (667.44 × 10 4 km 2 ). Among the 45 8-day periods, there were 24 8-day periods for which Criterion 2 produced a larger geographical area of estimated RSM than Criterion 1. The increased area of estimated RSM ranged from 0.30 × 10 4 km 2 on DOY 81 to 52.58 × 10 4 km 2 on DOY 321 and the average increased area was 19.77 × 10 4 km 2 [35]. Therefore, RSM estimation with Criterion 2 also improved the spatial coverage compared with Criterion 1.
Despite such differences in spatiotemporal coverage, estimated RSM with each criterion also shared some similarities. Both average RSM were highest (~22%) in autumn (DOYs 241-329), which agrees with the study by Jiao et al. (2016) that reported autumn had higher average SM than other seasons in 1998-2000 over the CLP [61]. In addition, we also found that spring (DOYs 57-145) had higher RSM than winter (DOYs 1-49 and 337-361) and summer (DOYs 153-233) (Figure 9).
It should be noted that the area produced by Criterion 1 could occasionally be greater than that of Criterion 2 (e.g., DOYs 49,105,153,161,249,and 289). The used optimal NDVI ATI , NDVI TVDI , R in validation, and estimated RSM area with each criterion for these 8-day periods are shown in Table 3. Models performed better with Criterion 2 (higher R in validation) despite a greater estimated RSM area with Criterion 1. Thus, we could not conclude that Criterion 1 was better than Criterion 2 by merely considering the spatial coverage. In addition, as we mentioned in Section 4.1.1, to acquire wider coverage of RSM estimation, we could choose the additional NDVI thresholds to estimate RSM, which corresponds to higher R rather than the highest R. Therefore, we might test more additional NDVI thresholds in subregions to combine the overall RSM map with the desired accuracy for further study.

Evaluation of Estimated RSM at the Station Scale
At the station scale, in order to reflect the overall result, stations should be uniformly distributed as much as possible with diverse geographical characterization (e.g., precipitation, elevation, and land cover). In addition, stations with more than 23 8-day periods (half of the 46 8-day periods) of estimated RSM could be considered to better reveal the variation of the estimated RSM and the observed RSM. In this case, six stations (station 53,553, station 53,771, station 53,845, station 53,857, station 57,031, and station 57,048) were selected from those stations as samples. The estimated RSM with each criterion and observed RSM with the total 8-day precipitation of six stations (described in Table 4) are presented in Figure 10. A detailed comparison reveals that despite missing some of the estimated RSM, especially for Criterion 1, the values of estimated and observed RSM have the same tendency. Importantly, the estimated RSM (closer to the observation) with Criterion 2 kept a better trend with the observed RSM at the station scale. Moreover, more estimated RSM with Criterion 2 was observed than with Criterion 1 throughout the period among six stations (the blue lines look more consistent than the red lines in Figure 10). The maximum estimated RSM values (near DOY 281) were observed in autumn among the six stations, consistent with heavy rainfall at that time.

Estimated Monthly RSM
The maps of monthly estimated RSM via Criterion 2 are demonstrated in Figure 11. The change of color from orange, red, green, and blue represents the gradual increase of RSM. RSM was higher in the southern part of the CLP in winter (December, January, and February), and in the western area in spring (March, April, and May) based on the clustering of the blue color. The soil in autumn (September, October, and November) was generally wetter than other seasons because most of the area was covered by green. The wetter soil might be mainly due to concentrated precipitation in autumn. Obviously, the estimated RSM maps in winter were more complete. Some researchers applied other dryness indices (e.g., PDI) or modified TVDI to estimate RSM [16,33,63,64]. To obtain more complete (monthly) RSM maps, ATI and TVDI could be replaced by other dryness indices to estimate RSM with Criterion 2 in our future research.  [62] and RSM might be increased dramatically after rainfall due to a delayed response of RSM changes to rainfall [13]. Moreover, high estimated RSM was observed in the periods when there were records of rainfall events. Meanwhile, one limitation of this study was that there were some temporal differences between in situ observations and the remote sensing imagery used.

Evaluation of Estimated RSM with Criterion 2 4.3.1. Estimated Monthly RSM
The maps of monthly estimated RSM via Criterion 2 are demonstrated in Figure 11. The change of color from orange, red, green, and blue represents the gradual increase of RSM. RSM was higher in the southern part of the CLP in winter (December, January, and February), and in the western area in spring (March, April, and May) based on the clustering of the blue color. The soil in autumn (September, October, and November) was generally wetter than other seasons because most of the area was covered by green. The wetter soil might be mainly due to concentrated precipitation in autumn. Obviously, the estimated RSM maps in winter were more complete. Some researchers applied other dryness indices (e.g., PDI) or modified TVDI to estimate RSM [16,33,63,64]. To obtain more complete (monthly) RSM maps, ATI and TVDI could be replaced by other dryness indices to estimate RSM with Criterion 2 in our future research. Importantly, the number of RSM stations for RSM verification varied among months. The available RSM stations are related not only to observation throughout the month but also to estimation in that month. In this case, there are fewer RSM stations for validation in January and February because of no RSM observations in the frozen soil and in April and October because of relatively incomplete RSM maps, respectively. The maximum station for validation was 180 in November and the minimum station was 59 in January and April ( Figure 12). There is a rule of thumb for interpreting the size of the correlation coefficient using the absolute value of the Pearson's r: 0.00-0.30, very weak; 0.30-0.50, weak; 0.50-0.70, moderate; 0.70-0.90, strong; and 0.90-1.00, very strong [65,66]. Estimated RSM had a weak correlation with the observed RSM in most months and the highest Pearson's r was 0.68 in January. The root mean square error (RMSE) varied from 3.77% in April to 6.10% in October. The highest mean absolute error (MAE) also appeared in October (4.97%) and the lowest MAE was 3.02% in March. Importantly, the number of RSM stations for RSM verification varied among months. The available RSM stations are related not only to observation throughout the month but also to estimation in that month. In this case, there are fewer RSM stations for validation in January and February because of no RSM observations in the frozen soil and in April and October because of relatively incomplete RSM maps, respectively. The maximum station for validation was 180 in November and the minimum station was 59 in January and April ( Figure 12). There is a rule of thumb for interpreting the size of the correlation coefficient using the absolute value of the Pearson's r: 0.00-0.30, very weak; 0.30-0.50, weak; 0.50-0.70, moderate; 0.70-0.90, strong; and 0.90-1.00, very strong [65,66]. Estimated RSM had a weak correlation with the observed RSM in most months and the highest Pearson's r was 0.68 in January. The root mean square error (RMSE) varied from 3.77% in April to 6.10% in October. The highest mean absolute error (MAE) also appeared in October (4.97%) and the lowest MAE was 3.02% in March. Remote Sens. 2021, 13, x FOR PEER REVIEW 20 of 26 The estimated monthly RSM and their area is illustrated in Figure 13. The average RSM varied from 8.64% in February to 16.42% in August and the highest and lowest maximum RSM (60.96% and 20.93%) appeared in September and January, respectively. The minimum RSM value was 0.00% for most months. However, there were two peaks in spring and the average RSM was clustered in the peak valley (11.42% in March, 14.18% in April, and 10.50% in May, respectively). The total area of the generated RSM area in February was the greatest (62.39 × 10 4 km 2 ). The area of generated RSM in April, July, and October was quite small (less than half the area of the CLP). Through comparison of the area distribution of the estimated RSM to the frequency of the observed RSM at available stations, similar tendencies were demonstrated for most months except in January and February. This might be attributed to the scarce RSM observations at that time, thus failing to reflect the spatial distribution of the whole CLP.   The estimated monthly RSM and their area is illustrated in Figure 13. The average RSM varied from 8.64% in February to 16.42% in August and the highest and lowest maximum RSM (60.96% and 20.93%) appeared in September and January, respectively. The minimum RSM value was 0.00% for most months. However, there were two peaks in spring and the average RSM was clustered in the peak valley (11.42% in March,14.18% in April, and 10.50% in May, respectively). The total area of the generated RSM area in February was the greatest (62.39 × 10 4 km 2 ). The area of generated RSM in April, July, and October was quite small (less than half the area of the CLP). Through comparison of the area distribution of the estimated RSM to the frequency of the observed RSM at available stations, similar tendencies were demonstrated for most months except in January and February. This might be attributed to the scarce RSM observations at that time, thus failing to reflect the spatial distribution of the whole CLP. Remote Sens. 2021, 13, x FOR PEER REVIEW 21 of 26

Estimated Seasonal and Yearly RSM
The seasonal and yearly RSM maps with the statistics of RSM are shown in Figure  14. Although estimated RSM had a weak correlation with the observed RSM among seasons (Pearson's r ranges from 0.53 in spring to 0.67 in winter and autumn), there was a moderate correlation (0.73) for annual RSM. RMSE varied from 3.74% in winter to 4.41% in autumn. The highest MAE also appeared in autumn (3.64%) and the lowest MAE was 3.00% for annual RSM. We found that autumn had the greatest error among the four seasons, which was highly related to the corresponding months (September, October, and November) with the higher error of RSM estimation ( Figure 12). In addition, the generated RSM area of these months in autumn was small and the merged seasonal RSM map would also have a larger error because the seasonal RSM of a certain pixel might be computed by the estimated RSM of a specific month instead of all months in autumn.
The available RSM stations for validation were 49 in the year 2017, which means only 49 observation stations have continuous RSM observations throughout the year. Meanwhile, the scarce RSM observations in winter and 2017 could not reflect the spatial distribution of the estimated RSM for the whole CLP (Figure 14a3,e3). The seasonal difference of RSM over the CLP is distinct in 1998-2000 and 2008-2010 [61], which agrees with our results especially for summer (dryer in the southeastern area and wetter in the remaining regions). The southeastern regions were affected by drought and extreme drought after analyzing the drought variation trends in different subregions of the CLP over four decades [67]. However, mean RSM was highest in autumn and lowest in spring [61], which

Estimated Seasonal and Yearly RSM
The seasonal and yearly RSM maps with the statistics of RSM are shown in Figure 14. Although estimated RSM had a weak correlation with the observed RSM among seasons (Pearson's r ranges from 0.53 in spring to 0.67 in winter and autumn), there was a moderate correlation (0.73) for annual RSM. RMSE varied from 3.74% in winter to 4.41% in autumn. The highest MAE also appeared in autumn (3.64%) and the lowest MAE was 3.00% for annual RSM. We found that autumn had the greatest error among the four seasons, which was highly related to the corresponding months (September, October, and November) with the higher error of RSM estimation ( Figure 12). In addition, the generated RSM area of these months in autumn was small and the merged seasonal RSM map would also have a larger error because the seasonal RSM of a certain pixel might be computed by the estimated RSM of a specific month instead of all months in autumn.
The available RSM stations for validation were 49 in the year 2017, which means only 49 observation stations have continuous RSM observations throughout the year. Meanwhile, the scarce RSM observations in winter and 2017 could not reflect the spatial distribution of the estimated RSM for the whole CLP (Figure 14(a3,e3)). The seasonal difference of RSM over the CLP is distinct in 1998-2000 and 2008-2010 [61], which agrees with our results especially for summer (dryer in the southeastern area and wetter in the remaining regions). The southeastern regions were affected by drought and extreme drought after analyzing the drought variation trends in different subregions of the CLP over four decades [67]. However, mean RSM was highest in autumn and lowest in spring [61], which is slightly different from our results (highest (13.91%) in autumn but lowest (9.08%) in winter). The reason for such a difference might be that they just compared SM in three seasons (spring, summer, and autumn) while we evaluated all four seasons [61]. is slightly different from our results (highest (13.91%) in autumn but lowest (9.08%) in winter). The reason for such a difference might be that they just compared SM in three seasons (spring, summer, and autumn) while we evaluated all four seasons [61].

Conclusions
This study aimed to improve the relative soil moisture (RSM) estimation that is based on the apparent thermal inertia (ATI) and temperature vegetation dryness index (TVDI). By optimizing the identification of NDVI thresholds, Criterion 2 (NDVI ATI <NDVI TVDI ) improved both the accuracy and spatiotemporal coverage of estimation compared with Criterion 1 (NDVI 0 ≤ NDVI ATI ≤ NDVI TVDI ). In addition, monthly, seasonal, and yearly RSM maps of the Chinese Loess Plateau (CLP) in 2017 were produced via the 8-day RSM maps and then examined. From our results, we conclude that:

•
The ATI/TVDI joint models not only have higher applicability than the ATI-based and TVDI-based models for all 8-day periods but also for simultaneous use within different NDVI ranges in the ATI/TVDI subregions for one 8-day period. Thus, in addition to the optimal NDVI thresholds, the additional NDVI thresholds we applied was another improved strategy to acquire wider spatial coverage of RSM estimation. • NDVI thresholds were optimized for robust RSM estimation with Criterion 2 for each 8-day period over the CLP and the selected optimal thresholds constantly changed throughout the study period. The applicability of Criterion 2, involving spatiotemporal coverage (45 and 38 8-day periods of RSM maps and the total RSM area of 939.52 × 10 4 km 2 and 667.44 × 10 4 km 2 with Criterion 2 and Criterion 1, respectively) and the accuracy (maximum R of 0.82 ± 0.007 for Criterion 2 and of 0.75 ± 0.008 for Criterion 1) of estimated RSM, was better than that of Criterion 1.

•
The estimated RSM (closer to the observation) with Criterion 2 kept a better trend with the observed RSM at the station scale. Moreover, more estimated RSM with Criterion 2 was observed than with Criterion 1 throughout the period.

•
High estimated RSM was observed in the periods when there were records of rainfall events, especially in autumn (mean RSM of 13.91 ± 2.65%)-wetter than other seasons. With a mean annual RSM of 10.16 ± 2.21%, the annual RSM map shows dryer areas in the southeastern part of the CLP.
This study focused on improving RSM estimation from MODIS imagery with Criterion 2, and such effort is still fundamental to the general research of SM remote sensing. We improved the mapping of RSM using Criterion 2. The next steps following this research would extend the retrieval to other sensors like Landsat and Sentinel satellites, to perform reliable estimates of SM at a high spatiotemporal resolution.  Informed Consent Statement: Not applicable for this study.

Data Availability Statement:
The input data used in this research can be accessed freely from online sources.