Soil Moisture Estimation for the Chinese Loess Plateau Using MODIS-derived ATI and TVDI

Abstract

Timely and effective estimation and monitoring of soil moisture (SM) provides not only an understanding of regional SM status for agricultural management or potential drought but also a basis for characterizing water and energy exchange. The apparent thermal inertia (ATI) and Temperature Vegetation Dryness Index (TVDI) are two widely used indices to reflect SM from remote sensing data. While the ATI-based model is routinely used to estimate the SM of bare soil and sparsely vegetated areas, the TVDI-based model is more suitable for areas with dense vegetation coverage. In this study, we present an iteration procedure that allows us to identify optimal Normalized Difference Vegetation Index (NDVI) thresholds for subregions and estimate their relative soil moisture (RSM) using three models (the ATI-based model, the TVDI-based model, and the ATI/TVDI joint model) from 1 January to 31 December 2017, in the Chinese Loess Plateau. The initial NDVI (NDVI₀) was first introduced to obtain TVDI value and two other thresholds of NDVI_ATI and NDVI_TVDI were designed for dividing the whole area into three subregions (the ATI subregion, the TVDI subregion, and the ATI/TVDI subregion). The NDVI values corresponding to maximum R-values (correlation coefficient) between estimated RSM and in situ RSM measurements were chosen as optimal NDVI thresholds after performing as high as 48,620 iterations with 10 rounds of 10-fold cross-calibration and validation for each period. An RSM map of the whole study area was produced by merging the RSM of each of the three subregions. The spatiotemporal and comparative analysis further indicated that the ATI/TVDI joint model has higher applicability (accounting for 36/38 periods) and accuracy than the ATI-based and TVDI-based models. The highest average R-value between the estimated RSM and in situ RSM measurements was 0.73 ± 0.011 (RMSE—root mean square error, 3.43 ± 0.071% and MAE—mean absolute error, 0.05 ± 0.025) on the 137th day of 2017 (DOY—day of the year, 137). Although there is potential for improved mapping of RSM for the entire Chinese Loess Plateau, the iteration procedure of identifying optimal thresholds determination offers a promising method for achieving finer-resolution and robust RSM estimation in large heterogeneous areas.

1. Introduction

Soil moisture (SM) is a key hydrological variable influencing water availability for vegetation and plays a fundamental role in land-atmosphere interactions [1]. Various remote sensing data, spanning the electromagnetic spectrum from microwave to optical range, have been tested for SM estimation since the 1970s [2,3,4]. Particularly, the use of optical and thermal infrared satellite images, with their fine spatial resolutions and wide coverage, has offered the potential for evaluating surface/sub-surface SM [5,6,7].

In the optical and thermal infrared regions of the electromagnetic spectrum, a range of SM indicators (e.g., TVDI—Temperature Vegetation Dryness Index and ATI—apparent thermal inertia) has successfully been applied for SM retrieval. TVDI is derived from vegetation indices (e.g., NDVI—normalized difference vegetation index) and land surface temperature (LST) [8,9,10], whereas ATI is obtained from LST and albedo [11,12,13,14,15,16]. Soil thermal inertia (TI) is a factor that reflects the temperature change of soil surface and is highly correlated with SM [17,18]. Based on this principle, simplified TI, namely ATI, is routinely used to estimate the SM of bare and less-vegetated soil [19,20,21,22,23,24,25,26]. There are no strict limitations on the number and time of day of LST observations using the ATI-based model comparing to other methods [27]. In addition, based on the relationship between remotely sensed LST and NDVI (generally called NDVI-LST space) [28], Sandholt et al. defined TVDI [29], which can represent the degree of SM in the area with vegetation coverage. Since then, this index has been adopted by many researchers to estimate SM from Moderate Resolution Imaging Spectroradiometer (MODIS) data on global and regional scales [18,30,31,32,33,34,35]. The potential of the MODIS-derived land surface dryness index for SM estimation over large geographical areas has been well demonstrated in specific seasons [20,36]. During the early stages of crop growth, the monitoring accuracy of ATI is better than that of TVDI. However, as crop growth progresses, the advantages of TVDI become evident. Therefore, it is inappropriate to employ a single model throughout the year because a single model (e.g., ATI or TVDI) generally is suitable for a similar land cover class rather than seasonal variation in surface coverage [27,37,38]. Wang et al. provided a viable method for time series SM monitoring and overcome the limitations of a single method (PDI—perpendicular drought index or TVDI) [10]. To date, few attempts have been made to test the joint use of the ATI-based and TVDI-based models and whether such joint can yield a more accurate SM estimation should be further investigated [17,39,40].

The scatter plots of remotely sensed LST and the NDVI often exhibit triangular or trapezoidal shapes and are called the NDVI-LST space. As an effective method based on the NDVI-LST feature space, the TVDI-based model considers vegetation coverage in SM estimation and has been widely applied to vegetated areas. Previous studies related the NDVI-LST triangle or trapezoid space with the linear fitting of the dry edge (the boundary derived based on observed maximum LST for the same NDVI representing limited water availability) and wet edge (the boundary derived based on observed minimum LST for the same NDVI representing unlimited water access) to represent SM variations [41,42,43]. Several researchers proposed parabolic or logarithmic fitting of the dry/wet edge in the space [44,45] or applied other types of vegetation index (VI) for an improved space [31,46,47]. However, the initial NDVI (NDVI₀) in the NDVI-LST space is rarely considered in constructing dry/wet edges, applying different NDVI₀ will derive different dry/wet edges, thereby further affecting SM estimation. The NDVI₀ is the low limit of NDVI, below which the data are excluded and should be tested for SM retrieval from the TVDI-based model.

K-fold cross-validation is a sampling technique for model validation when observation data are limited and can also be used for evaluating the performance of a machine-learning model [48]. It is known that cross-validation provides an unbiased estimation result because every observation has an opportunity to appear in training and testing datasets [49,50]. Cross-validation on models demonstrated the potential of radar and LANDSAT data for the estimation of surface moisture [51,52]. Baldwin et al. run a 20-fold cross-validation procedure to test the predictive skill of Soil Moisture Analytical Relationship infiltration model-Kalman filter (SMAR-EnKF) for the parameter regressions [53]. Relative soil moisture (RSM) represents the percentage of SM that accounts for the moisture storage capacity and is used to describe the soil moisture levels in the present study. It would be interesting to adopt a cross-validation method in the process of calibration (e.g., a 10-fold cross-calibration [54]) to examine the relationship between dryness indices (ATI or TVDI) and in situ RSM measurements, since they are random variables. After the 10th fold calibration, the correlation coefficient (R) between estimated RSM and in situ RSM measurements is calculated to evaluate the performance of the model (e.g., the ATI-based model).

Deep and loose loess soil in the Chinese Loess Plateau can retain a high amount of water, while high soil porosity causes strong evaporation and scarcity of surface water. Large-scale dryness detecting and monitoring in dry areas are essential for improving agricultural management and formulating policies that support drought risk reduction. As mentioned above, a new method is needed to not only estimate soil moisture accurately but also help to address the above problems. The main goal of the presented study is to improve soil moisture estimation for the Chinese Loess Plateau using the ATI-based model, the TVDI-based model, and the ATI/TVDI joint model. Three optimal NDVI thresholds are identified for subregional RSM retrieval. Then, the mapping the retrieved RSM for each period from 1 January to 31 December 2017, and spatial distribution of RSM will be also presented by season.

2. Study Area and Data

2.1. Study Area

The study area is the Chinese Loess Plateau (100°54′–114°33′E, 33°43′N-41°16′N, Figure 1) in northwestern China, covering a total area of approximately 6.4 × 10⁵ km² and ranging from 200 m to 3000 m in elevation. This area is a typical water-limited terrestrial ecosystem with a temperate continental monsoon climate. The long-term mean annual temperature of the CLP is 4 °C in the northwest and 14 °C in the southeast. The annual precipitation ranges from 200 mm in the northwest to 750 mm in the southeast [55] and occurs in the form of thunderstorms during summer months from July to September (accounting for 60% of the annual precipitation). The dominant porous loess soil and heavy rainfall make the CLP prone to erosion. The main land use types include native grassland, farmland, and introduced vegetated land. There are 213 Chinese automatic soil moisture observation stations (CASMOS) on the CLP (Figure 1).

2.2. Satellite Data and Image Pre-processing

A series of 6 Collection Terra Moderate Resolution Imaging Spectroradiometer (MODIS) products were freely downloaded from the website of Level-1 and Atmosphere Archive and Distribution System (LAADS) Distributed Archive Center (DAAC) [56]. The Terra MODIS image data, imaged in the optical and thermal infrared region of the electromagnetic spectrum in the year of 2017 (Julian day 1 to 365), were used to estimate RSM in this study. The satellite products are summarized in

Table 1. MODIS data used in the study.

Data Product	Variables Used	Use	Spatial/Temporal Resolution
MOD09GA	The temporal range for granule acquisition	Collecting the corresponding times of in situ RSM measurements	500 m, 1 day
MOD09A1	Surface reflectance	Calculating NDVI, ATI, and TVDI	500 m, 8 day
MOD11A2	Daytime/nighttime LST	Calculating TVDI and ATI	1 km, 8 day
MCD12Q1. Type2	Land covers	Producing a water body mask	1 km, 1 year

. 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, which includes quality flags and information on influence from clouds. Five granules of MODIS data (h25v04, h26v04, h25v05, h26v05, and h27v05) were used to cover the CLP for every single period and each 8-day composite pixel contains the best possible observation according to specified criteria [57]. Five granules were then mosaicked and re-sampled at the resolution of 1/224° (~500m) and re-projected into Geographic Coordinate System on the WGS 84 geoid using the MODIS Reprojection Tool (MRT) (https://lpdaac.usgs.gov/tools/modis_reprojection_tool/), and were clipped using the vector borders of the CLP in ArcGIS 10.2 (ESRI Inc., Redlands, CA, USA). In the following experiments, MODIS-derived NDVI was calculated and after the quality flagging procedure, the daytime LST was used to compute TVDI [58]. In addition, the variable ATI was produced with surface reflectance and the difference between LST of daytime and nighttime.

2.3. In Situ Measured RSM Data

Hourly RSM at 20 cm depth measurements at 213 automatic SM observation stations (Figure 1; see

Table A1. Locations of the automatic RSM observation stations within the Chinese Loess Plateau.

Serial Number	Station	Latitude (°E)	Longitude (°N)	Serial Number	Station	Latitude (°E)	Longitude (°N)
1	52,765	37.38	101.61	44	53,488	40.04	113.59
2	52,797	37.19	104.06	45	53,490	40.44	114.06
3	52,855	36.69	101.25	46	53,512	39.79	106.80
4	52,862	36.97	101.66	47	53,513	40.72	107.37
5	52,863	36.82	101.95	48	53,522	40.05	107.84
6	52,866	36.73	101.75	49	53,529	39.09	107.96
7	52,868	36.02	101.37	50	53,533	39.81	108.71
8	52,869	36.49	101.58	51	53,543	39.82	110.01
9	52,874	36.49	102.41	52	53,545	39.56	109.71
10	52,875	36.5	102.10	53	53,547	39.10	109.03
11	52,876	36.33	102.84	54	53,553	39.85	111.22
12	52,877	36.10	102.26	55	53,562	39.92	111.66
13	52,884	36.35	103.93	56	53,564	39.37	111.21
14	52,885	36.75	103.25	57	53,565	39.44	111.50
15	52,895	36.57	104.69	58	53,574	39.52	112.27
16	52,896	36.55	104.15	59	53,576	39.51	112.82
17	52,963	35.94	102.02	60	53,578	39.37	112.43
18	52,972	35.85	102.46	61	53,579	39.05	112.92
19	52,974	35.54	102.03	62	53,580	39.83	113.10
20	52,980	35.97	103.30	63	53,582	39.71	113.67
21	52,982	35.48	103.56	64	53,584	39.56	113.16
22	52,983	35.87	104.14	65	53,585	39.17	113.26
23	52,984	35.58	103.18	66	53,590	39.74	114.26
24	52,985	35.41	103.34	67	53,594	39.44	114.22
25	52,986	35.36	103.86	68	53,644	38.59	108.83
26	52,993	35.68	105.06	69	53,646	38.27	109.78
27	52,995	35.58	104.60	70	53,651	38.83	110.46
28	52,998	35.13	104.20	71	53,659	37.97	111.01
29	53,337	41.05	108.28	72	53,662	38.71	111.57
30	53,348	41.02	109.13	73	53,663	38.93	111.82
31	53,419	40.33	106.99	74	53,665	38.28	111.63
32	53,420	40.87	107.13	75	53,666	38.36	111.94
33	53,433	40.73	108.65	76	53,669	38.07	111.81
34	53,446	40.53	109.88	77	53,673	38.74	112.71
35	53,455	40.55	110.55	78	53,674	38.39	112.69
36	53,457	40.39	110.03	79	53,676	38.5	112.98
37	53,463	40.86	111.57	80	53,677	37.94	112.48
38	53,464	40.73	111.17	81	53,678	38.07	112.65
39	53,466	40.76	111.71	82	53,679	37.75	112.54
40	53,467	40.25	111.25	83	53,681	38.84	113.36
41	53,469	40.4	111.82	84	53,685	38.08	113.42
42	53,478	39.99	112.46	85	53,687	37.79	113.63
43	53,486	40.37	113.77	86	53,725	37.6	107.60
87	53,730	38.19	107.47	130	53,882	36.06	113.03
88	53,732	37.86	108.72	131	53,884	36.51	113.03
89	53,735	37.6	108.80	132	53,888	36.20	113.44
90	53,738	36.92	108.18	133	53,906	35.52	105.71
91	53,740	37.96	109.29	134	53,908	35.21	105.23
92	53,748	37.15	109.69	135	53,915	35.53	106.66
93	53,754	37.49	110.26	136	53,917	35.22	106.06
94	53,759	36.99	110.83	137	53,923	35.73	107.63
95	53,763	37.90	112.17	138	53,924	35.07	107.62
96	53,764	37.51	111.11	139	53,925	35.68	107.19
97	53,767	37.33	111.18	140	53,926	35.34	107.35
98	53,768	37.15	111.75	141	53,927	35.20	106.62
99	53,769	37.24	111.78	142	53,928	35.30	107.02
100	53,770	37.36	112.35	143	53,929	35.21	107.77
101	53,771	37.41	112.05	144	53,930	36.45	107.99
102	53,774	37.58	112.37	145	53,931	36.00	109.38
103	53,775	37.42	112.59	146	53,934	35.78	107.98
104	53,776	37.70	112.79	147	53,937	35.53	107.89
105	53,777	37.51	112.14	148	53,938	35.17	108.28
106	53,778	37.17	112.18	149	53,941	35.19	109.58
107	53,780	37.91	113.15	150	53,942	35.79	109.36
108	53,783	37.60	113.72	151	53,947	35.06	109.08
109	53,788	37.33	113.57	152	53,948	34.89	109.63
110	53,821	36.57	107.30	153	53,950	35.23	110.15
111	53,845	36.60	109.50	154	53,954	35.62	110.97
112	53,853	36.70	110.95	155	53,955	35.52	110.46
113	53,856	36.46	110.75	156	53,956	35.42	110.84
114	53,857	36.06	110.18	157	53,957	35.61	110.71
115	53,859	36.09	110.66	158	53,958	35.17	110.78
116	53,861	35.89	111.39	159	53,959	35.11	111.07
117	53,863	37.06	111.94	160	53,961	35.65	111.49
118	53,865	36.66	111.56	161	53,962	35.73	111.70
119	53,866	36.23	111.66	162	53,963	35.65	111.36
120	53,868	36.06	111.49	163	53,964	35.62	111.21
121	53,869	36.58	111.70	164	53,965	35.50	111.57
122	53,871	36.86	112.87	165	53,966	35.98	111.83
123	53,872	36.76	112.69	166	53,967	35.33	111.20
124	53,873	36.10	112.87	167	53,968	35.28	111.66
125	53,874	36.25	111.90	168	53,970	35.69	112.19
126	53,877	36.16	112.25	169	53,973	35.78	112.94
127	53,878	36.51	113.36	170	53,975	35.49	112.41
128	53,879	36.32	112.89	171	53,976	35.50	112.86
129	53,880	36.34	113.23	172	53,978	35.09	112.63
173	56,092	34.99	104.65	194	57,042	34.78	109.19
174	57,001	34.75	105.33	195	57,043	34.80	109.97
175	57,003	34.89	106.84	196	57,044	34.40	109.23
176	57,004	34.73	104.88	197	57,045	34.40	109.49
177	57,014	34.56	105.87	198	57,047	34.16	109.32
178	57,020	34.36	107.39	199	57,048	34.40	108.72
179	57,022	34.68	107.80	200	57,051	34.79	111.19
180	57,024	34.44	107.65	201	57,052	34.88	110.45
181	57,025	34.51	107.38	202	57,053	34.70	110.71
182	57,026	34.37	107.88	203	57,055	34.59	110.13
183	57,027	34.30	107.73	204	57,060	35.16	111.22
184	57,029	34.49	108.45	205	57,061	34.84	111.20
185	57,030	34.69	108.14	206	57,066	34.40	111.67
186	57,031	34.83	108.55	207	57,071	34.80	112.47
187	57,032	34.14	108.21	208	57,074	34.42	112.40
188	57,033	34.56	108.53	209	57,076	34.74	112.79
189	57,034	34.31	108.24	210	57,080	34.74	112.97
190	57,035	34.54	108.25	211	57,123	34.29	108.07
191	57,037	34.92	108.98	212	57,131	34.45	108.97
192	57,038	34.28	108.48	213	57,132	34.14	108.58
193	57,041	34.63	108.93

for their locations) in 2017 were obtained from the Chinese Meteorological Data Service Center (CMDC) [59]. It is worth mentioning that the number of observation stations varied temporarily with maximal 213 stations available among the 46 periods in 2017 (DOY—day of the year, 1-DOY 361), which includes 45 periods from the first day to the 360th day of 2017 with an 8-day interval and one period (DOY 361) composited with the last five days of the year. For an accurate temporal match between the in situ measurement data and the 8-day composite products, the daily granule acquisition time of the MOD09GA products (from the beginning to ending date-time) should be collected first to serve as the reference for selecting corresponding in situ RSM measurements. In particular, the revisit period of the Terra satellite is 16 days, and the varying acquisition time daily within 16 days is observed according to viewing information from the website (https://ladsweb.modaps.eosdis.nasa.gov/). The acquisition times of the MODIS granules and their corresponding in situ RSM measurements on DOY 113 and at Station 52,765 are shown in

Table 2. The acquisition times of the MODIS granules and their corresponding in situ RSM measurements on DOY 113 and at Station 52,765.

Day	The Acquisition Times of the MODIS Granules						In Situ RSM Measurements
	Time Range	h25v04 (UTC)	h25v05 (UTC)	h26v04 (UTC)	h26v05 (UTC)	h27v05 (UTC)	Time Range (UTC)	Daily RSM (%)
23 Apr.	Start time	02:05:00	03:45:00	02:05:00	03:45:00	02:05:00	03:00:00	25.00
	End time	05:30:00	05:30:00	03:50:00	03:55:00	03:55:00	06:00:00
24 Apr.	Start time	02:45:00	-	01:10:00	02:50:00	02:50:00	02:00:00	25.00
	End time	02:55:00	-	02:55:00	03:00:00	03:00:00	03:00:00
25 Apr.	Start time	01:50:00	03:35:00	01:50:00	03:35:00	01:55:00	02:00:00	25.40
	End time	05:15:00	05:20:00	03:40:00	05:20:00	03:40:00	06:00:00
26 Apr.	Start time	02:35:00	04:15:00	00:55:00	02:40:00	02:35:00	01:00:00	26.00
	End time	04:20:00	06:05:00	04:20:00	04:25:00	04:25:00	07:00:00
27 Apr.	Start time	03:20:00	03:20:00	01:40:00	03:20:00	01:40:00	02:00:00	25.00
	End time	05:05:00	05:10:00	03:25:00	05:05:00	03:30:00	06:00:00
28 Apr.	Start time	02:25:00	04:05:00	00:45:00	02:25:00	02:25:00	01:00:00	25.00
	End time	05:45:00	05:50:00	04:10:00	04:10:00	04:10:00	06:00:00
29 Apr.	Start time	03:05:00	03:10:00	01:25:00	03:10:00	01:30:00	02:00:00	24.00
	End time	04:55:00	04:55:00	04:50:00	04:55:00	03:15:00	05:00:00
30 Apr.	Start time	02:10:00	03:50:00	02:10:00	02:15:00	02:15:00	03:00:00	26.00
	End time	05:35:00	05:40:00	03:55:00	04:00:00	04:00:00	06:00:00
8-day average in situ RSM measurement (DOY: 113; Station: 52,765)								25.18

, for example. The 8-day average in situ RSM value on DOY 113 at Station 52,765 was 25.18%, which was calculated by averaging the daily RSM (the figures in the last column in Table 2).

3. Methods

3.1. ATI and TVDI

3.1.1. Apparent Thermal Inertia (ATI)

Soil thermal inertia (TI) is described as a thermal property of soil that characterizes its resistance to temperature change and has been used for near-surface SM retrieval [60,61]. However, the calculation of soil TI is relatively complex because of the difficulty in acquiring the values of the required parameters. Price proposed that solar radiation incident in the calculation should be treated as a constant under some conditions [11]. ATI, to simplify the TI, is calculated by spectral surface albedo and the diurnal LST range for RSM retrieval [62]. ATI normally ranges from 1 to 0, estimated by the following equation [63]:

$$\mathrm{ATI}=\frac{1-\mathrm{A}}{\Delta \mathrm{LST}}$$

(1)

$$\Delta \mathrm{LST}=LS{T}_{daytime}-LS{T}_{nighttime}$$

(2)

where ATI is apparent thermal inertia [K⁻¹], A is full band albedo, $\Delta \mathrm{LST}$ is the difference between $LS{T}_{daytime}$ and $LS{T}_{nighttime}$ [K], $LS{T}_{daytime}$ is the 8-day average daytime LST [K], and $LS{T}_{nighttime}$ is the 8-day average nighttime LST [K].

To estimate the full band albedo, Liang et al. proposed an empirical formula for MODIS image data [64]:

$$\mathrm{A}=0.16b_1+0.291b_2+0.243b_3+0.11b_4+0.112b_5+0.081b_7-0.0015$$

(3)

where $b_1$–$b_5$ and $b_7$ refer to the reflectance of MODIS bands.

3.1.2. Interpretation of the NDVI-LST Triangle Space

Previous researchers have noted a negative relationship between NDVI and LST, which usually forms a triangular shape [29,65,66,67] in the scatter plot of NDVI against LST. It can be related to vegetation coverage and SM, thus is widely used to estimate SM [31,34,68,69,70,71,72].

The conceptual NDVI-LST triangle space is shown in Figure 2. In the triangle, the edges in red or blue can be defined as representing extreme conditions of SM and evapotranspiration. While the decreasing dry edges in red (CE, C’E’) represent limiting conditions (maximum LST and no evapotranspiration) of SM and evapotranspiration for different land cover classes, the horizontal wet edges in blue (DE, D’E’) represent the LST (LST_min) under unlimited water access conditions and potential evapotranspiration. The pixels closest to the dry edge reflect extremely stressed surfaces with low SM. As green vegetation increases along the x-axis, LST_max declines. Theoretically, in the triangular space, the temperatures of full vegetation cover and bare soil are the same when soil moisture is saturated [46]. The left and right edges (CD, C’D’) represent bare soil, ranging from dry to wet (LST_max to LST_min, respectively). Thus, the green color lines in the NDVI-LST space can be considered as a set of SM contours and the location of a pixel in the NDVI-LST triangular space is influenced mainly by SM availability (Figure 2).

3.1.3. Temperature Vegetation Dryness Index (TVDI)

To obtain information on RSM, the spatial distribution of TVDI shows clearly that TVDI is higher (maximum 1) and soil moisture is lower (water-stressed) at the dry edge, whereas TVDI is lower (minimum 0) and soil moisture is higher (saturated) at the wet edge [44]. The dryness index TVDI acquired from the NDVI-LST triangle space is defined by Sandholt et al. [29] as follows:

$$\mathrm{TVDI}=\frac{\mathrm{LST}-{\mathrm{LST}}_{\min }}{{\mathrm{LST}}_{\max }-{\mathrm{LST}}_{\min }}$$

(4)

where LST is the MODIS-derived LST in a given pixel, LST_min refers to the minimum LST in the triangle space defining the wet edge at a given NDVI, LST_max is the maximum LST in the triangle space defining the dry edge at a given NDVI.

As demonstrated in Figure 2, both dry and wet edges in the NDVI-LST triangle space are linearly regressed for an area large enough and the corresponding formula are shown as:

$${\mathrm{LST}}_{\max }={\mathrm{a}}_{\mathrm{dry}} \times \mathrm{NDVI}+{\mathrm{b}}_{\mathrm{dry}}$$

(5)

$${\mathrm{LST}}_{\min }={\mathrm{a}}_{\mathrm{wet}} \times \mathrm{NDVI}+{\mathrm{b}}_{\mathrm{wet}}$$

(6)

with Equations (5) to (6), Equation (4) for TVDI can be transformed as follows [8,10,75,76]:

$$\mathrm{TVDI}=\frac{\mathrm{LST}-\left({\mathrm{a}}_{\mathrm{wet}} \times \mathrm{NDVI}+{\mathrm{b}}_{\mathrm{wet}}\right)}{\left({\mathrm{a}}_{\mathrm{dry}} \times \mathrm{NDVI}+{\mathrm{b}}_{\mathrm{dry}}\right)-\left({\mathrm{a}}_{\mathrm{wet}} \times \mathrm{NDVI}+{\mathrm{b}}_{\mathrm{wet}}\right)}$$

(7)

where ${\mathrm{a}}_{\mathrm{wet}}$, ${\mathrm{a}}_{\mathrm{dry}}$, ${\mathrm{b}}_{\mathrm{wet}},$ and ${\mathrm{b}}_{\mathrm{dry}}$ are the linear fitting coefficients (slope and intercept) defining the dry and wet edges.

3.1.4. Applying Initial NDVI₀ for Determination of TVDI

The scatterplot for the relationship between NDVI and LST, exemplified by DOY 113, is demonstrated in Figure 3. As the initial NDVI₀ increases, both dry edge and wet edges change slightly. Such relationships (LST_max or LST_min and NDVI) are indicated when the initial NDVI₀ equals 0, 0.1, and 0.2, respectively, with the linear equations in the scatterplot (Figure 3), positive slope values for the dry edge equations, and negative slope values for the wet edge equations. Thus, the TVDI value for a given pixel, to some extent, is influenced by the fluctuated wet or dry edge. Additionally, the TVDI-based model is more applicable to vegetated land with relatively high NDVI. In order to improve TVDI-based RSM retrieval, it is necessary to test NDVI₀ from 0 to 0.5 with a small interval of 0.01 in the study area. The new TVDI value was calculated based on a given NDVI₀ in an iteration. It is worth mentioning that the pixels of NDVI less than 0 were ignored in this study since these pixels are mainly water body, cloud, or snow,which are regarded as 100% RSM.

3.2. Establishment of the Retrieval Model

The linear relationship between ATI and in situ RSM measurements, and between TVDI and in situ RSM measurements is simple, but which subregion using the ATI-based model or TVDI-based model for the whole study area seems rather complex. In that case, it is challenging to get corresponding reliable thresholds for subregions using different models. We proposed an iteration procedure to obtain optimal NDVI thresholds for subregions. In detail, as can be seen in Figure 4, the initial NDVI (NDVI₀) should be used first to determine the TVDI-based model. Once the NDVI₀ is determined, the wet and dry edges and their linear regression equations (Equations (5) and (6)) will be formed with the specific value of a_dry, a_wet, b_dry, b_wet, R²_wet, and R²_dry. The second and third NDVI threshold (NDVI_ATI and NDVI_TVDI, respectively) are then used to divide the entire study area into three subregions, namely the ATI subregion (ATI), the TVDI subregion (TVDI), and the ATI/TVDI subregion (M_ATI/TVDI). While the area with NDVI less than NDVI_ATI was regarded as the ATI subregion, where only the ATI-based model was applied to the retrieval of RSM, the area with NDVI greater than NDVI_TVDI was regarded as the TVDI subregion where only the TVDI-based model was applied. For the ATI/TVDI subregion (NDVI_ATI ≤ NDVI ≤ NDVI_TVDI), the average of ATI and TVDI (M_ATI/TVDI) was used to estimate RSM. The procedure of calibration, that the linear regression between ATI and in situ RSM measurements, between M_ATI/TVDI and in situ RSM measurements, and between TVDI and in situ RSM measurements, was performed in each iteration for the three subregions.

Here, we took the ATI subregion as an example. Both NDVI₀ and NDVI_ATI ranged from 0 to 0.5 with a small interval of 0.01 and NDVI_TVDI ranged from 0 to 0.7 with an interval of 0.01 in the iterative optimization procedure. For one iteration, linear regression analysis was performed to examine the relation between ATI and in situ RSM measurements through the 10-fold crossing-calibration in the calibration process. In the validation procedure, ten sets of statistics between in situ RSM measurements and estimated RSM, including R-values, RMSE (root mean square error), MAE (mean absolute error), and the p-values for significance tests were computed with corresponding NDVI₀, NDVI_ATI, and NDVI_TVDI thresholds. The mathematical expressions of the goodness validation between in situ RSM measurements and estimated RSM in the study are presented in

Table A2. Mathematical expressions of the goodness of validation.

Abbreviation	Formula
R2	${\left(\frac{{{\displaystyle \sum }}_{i=0}^n\left(O_i-\overline{O_i}\right)\left(P_i-\overline{P_i}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-\overline{O_i}\right)}^2{{\displaystyle \sum }}_{i=1}^n{\left(P_i-\overline{P_i}\right)}^2}}\right)}^2$
R	$\frac{{{\displaystyle \sum }}_{i=0}^n\left(O_i-\overline{O_i}\right)\left(P_i-\overline{P_i}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-\overline{O_i}\right)}^2{{\displaystyle \sum }}_{i=1}^n{\left(P_i-\overline{P_i}\right)}^2}}$
RMSE	$\sqrt{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(O_i-P_i\right)}^2}$
MAE	$\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n\left|O_i-P_i\right|$
STD	$\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(P_i-\overline{P_i}\right)}^2}{n}}$

where $O_i$ and $P_i$ are the in situ RSM measurements and estimated RSM, $\overline{O}$ is the mean of in situ RSM measurements and $\overline{P}$ is the mean of estimated RSM.

. Then the average R values of ten rounds were calculated as measures of the goodness-of-fit between in situ RSM measurements and estimated RSM. The calibration and validation steps in the ATI subregion were performed simultaneously in the TVDI subregion and the ATI/TVDI subregion with corresponding NDVI₀, NDVI_ATI, and NDVI_TVDI thresholds.

After completing all iterations, the NDVI thresholds corresponding to maximum R-value were chosen as the optimal thresholds when NDVI₀, NDVI_ATI, and NDVI_TVDI thresholds exist simultaneously for RSM estimation in the whole CLP. The subregional RSM was estimated by applying selected optimal thresholds (NDVI₀, NDVI_ATI, and NDVI_TVDI) for each period. Generally, the whole region was divided into three subregions by selected NDVI_ATI and NDVI_TVDI, and the overall RSM map, ultimately produced using a combination of the subregional RSM for each period. The completely iterative optimization procedure was repeated for all 46 periods in 2017. The spatiotemporal distribution of RSM was characterized by derived RSM imagery.

In terms of establishing NDVI thresholds, it is worth mentioning that NDVI₀ should be smaller than or equal to NDVI_ATI, which is smaller than or equal to NDVI_TVDI (NDVI₀ ≤ NDVI_ATI ≤ NDVI_TVDI). In this way, the TVDI value for each pixel could be calculated in the TVDI subregion. Similarly, the mean of ATI and TVDI for each pixel could be assigned in the ATI/TVDI subregion.

3.2.1. RSM Estimation for the ATI Subregion

In the following experiment, the focus was on identifying an NDVI_ATI threshold which allows the CLP to be divided into two subregions first. The area with NDVI less than or equal to NDVI_ATI (hereinafter ATI subregion) was assigned ATI value. The NDVI_ATI threshold was selected from 0 to 0.5 with a small interval of 0.01. In general, ATI has a positive correlation with RSM, and the ATI-based model is mainly applied to barely or sparsely vegetate ground [77]. Thus, we chose the area where NDVI was less than the NDVI_ATI threshold as the ATI subregion for the application of the ATI-based model in the calibration process. The regression parameters (${\mathrm{a}}_{\mathrm{ATI}}$ and ${\mathrm{b}}_{\mathrm{ATI}}$) of the ATI subregion can be calculated using the following empirical linear equation:

$$\mathrm{RSM}={\mathrm{a}}_{\mathrm{ATI}} \times \mathrm{ATI}+{\mathrm{b}}_{\mathrm{ATI}}$$

(8)

where $\mathrm{RSM}$ is the relative soil moisture value (unit: %) measured at depth of 20 cm by automatic soil moisture observation stations, $\mathrm{ATI}$ is the $\mathrm{ATI}$ value for a given pixel (NDVI ≤ NDVI_ATI), ${\mathrm{a}}_{\mathrm{ATI}}$ and ${\mathrm{b}}_{\mathrm{ATI}}$ are the regression parameters using $\mathrm{ATI}$ value and in situ $\mathrm{RSM}$ measurements.

3.2.2. RSM Estimation for the TVDI Subregion

In this study, the TVDI-based model was also used for RSM estimation. Since ATI generally results in weaker correlations with in situ RSM measurements in vegetated areas than TVDI, it is better to consider the application of TVDI for RSM estimation to the TVDI subregion. To outline the TVDI subregion, an additional threshold of NDVI_TVDI was considered, which was greater than or equal to NDVI_ATI. The NDVI_TVDI threshold was selected from 0 to 0.7 with an interval of 0.01. The area for which NDVI was larger than or equal to NDVI_TVDI was regarded as the TVDI subregion and the TVDI value was applied to such an area for constructing a linear regression model with in situ RSM measurements in the calibration process. Therefore, the calibration of the model is presented as follows:

$$\mathrm{RSM}={\mathrm{a}}_{\mathrm{TVDI}} \times \mathrm{TVDI}+{\mathrm{b}}_{\mathrm{TVDI}}$$

(9)

where $\mathrm{TVDI}$ refers to the $\mathrm{TVDI}$ for a given pixel (NDVI ≥ NDVI_TVDI)$;{ \mathrm{a}}_{\mathrm{TVDI}}$ and ${\mathrm{b}}_{\mathrm{TVDI}}$ are the regression parameters using $\mathrm{TVDI}$ value and in situ $\mathrm{RSM}$ measurements

3.2.3. RSM Estimation for the ATI/TVDI Subregion

To achieve higher RSM estimation accuracy, both ATI and TVDI are utilized to retrieve the RSM of the remaining area of the CLP (i.e., the CLP except for the ATI subregion and TVDI subregion; hereinafter as the ATI/TVDI subregion), which overcomes the limitations of a single method (ATI-based or TVDI-based models). Zhao et al. (2012) proposed a model-level integrated approach to effectively retrieve large-scale daily soil moisture. They regarded the average value of soil moisture from the ATI-based model and the TVDI-based model as the soil moisture when NDVI ranged from 0.1 to 0.18 [78]. Wang et al. (2020) applied the soil moisture obtained by averaging the PDI (perpendicular drought index)-based model and TVDI-based model [10]. In our study, the average of ATI and TVDI value was calculated where NDVI varied from NDVI_ATI and NDVI_TVDI in the ATI/TVDI subregion, as a joint index M_ATI/TVDI ranging from 0 to 1. However, the nature of the correlation (positive or negative) between M_ATI/TVDI and in situ RSM measurements varied. The calibration equations for calculating M_ATI/TVDI and in situ RSM measurements in the ATI/TVDI subregion are given below:

$${\mathrm{M}}_{\mathrm{ATI}/\mathrm{TVDI}}=\frac{\mathrm{ATI}+\mathrm{TVDI}}{2}$$

(10)

$$\mathrm{RSM}={\mathrm{a}}_{\mathrm{ATI}/\mathrm{TVDI}} \times { \mathrm{M}}_{\mathrm{ATI}/\mathrm{TVDI}}+{\mathrm{b}}_{\mathrm{ATI}/\mathrm{TVDI}}$$

(11)

where ${\mathrm{M}}_{\mathrm{ATI}/\mathrm{TVDI}}$ is the average of ATI and TVDI for a given pixel (NDVI_ATI < NDVI < NDVI_TVDI) and ${\mathrm{a}}_{\mathrm{ATI}/\mathrm{TVDI}}$ and ${\mathrm{b}}_{\mathrm{ATI}/\mathrm{TVDI}}$ are the regression parameters using ${\mathrm{M}}_{\mathrm{ATI}/\mathrm{TVDI}}$ and in situ RSM measurements.

3.2.4. Calibration and Validation

The 10-fold cross-calibration method (traditionally called 10-fold cross-validation in statistics [79]) would be used to improve the accuracy of RSM estimation. In general, this method is used primarily for assessing how accurately a predictive model performs [80,81]. Moreover, to reduce the variability of the calibration, 10 rounds of 10-fold cross-calibration are performed between generated dryness indices (ATI, TVDI, M_ATI/TVDI) and in situ RSM measurements in those three subregions, respectively.

The ATI subregion was taken as an example to illustrate the calibration and validation procedures (Figure 5). The paired ATI and in situ RSM measurements were regarded as original samples and the RSM estimation was allowed only when the number of available RSM observation stations was greater than 20 in the ATI subregion. In detail, one round of cross-calibration involved random splitting of the original samples of paired ATI and in situ RSM measurements into complementary 10 subsamples with 9 as training data and one as testing data. The detailed calibration is shown in Figure 6. Training data in our analysis were calibrated using Equation (8) one-fold. The estimated RSM for one-fold was calculated using ATI value in testing data and obtained regression parameters (a_ATI and b_ATI) in the process of one-fold calibration. In this way, one group of estimated RSM and in situ RSM measurements was generated. We repeated this process until every 10-fold served as the testing data in practice. After the 10th fold was accomplished, we collected 10 groups of estimated RSM and in situ RSM measurements as the validation data. The validation results, including R_ATI, RMSE, and MAE, as well as the p-value for a significance test, were computed, which was just called “one-round”.

The original samples of paired ATI and in situ RSM measurements would be randomly split into 10 subsamples and the same calibration and validation were performed again. After ten rounds of iterations, we averaged the ten rounds of validation outcomes of the R ($\overline{R_{ATI}}$) as the reference to evaluate corresponding NDVI thresholds, since the ATI subregion was generated from NDVI_ATI. It should be noted that the same procedures of calibration and validation were conducted in the TVDI subregion as well as in the ATI/TVDI subregion. Similarly, the $\overline{R_{TVDI}}$ value was determined by the thresholds NDVI₀ and NDVI_TVDI. The three thresholds (NDVI₀, NDVI_ATI, and NDVI_TVDI) together determined the value of $\overline{R_{ATI/TVDI}}$.

3.2.5. Identifying and Applying the Optimal NDVI Thresholds for Subregional RSM Retrieval

The optimal NDVI thresholds were identified dependent on the average validation results of the correlation coefficient ($\overline{R_{ATI}}$, $\overline{R_{ATI/TVDI}}$ and $\overline{R_{TVDI}}$), which directly impact the accuracy of RSM retrieval. It is worth mentioning that the NDVI thresholds corresponding to maximum R-value among three subregions were chosen as the optimal thresholds when NDVI₀, NDVI_ATI, and NDVI_TVDI thresholds exist simultaneously for one period. Here, we could not guarantee that all three subregions reach the highest averaged R-value in their subregions at the same time with selected optimal NDVI thresholds. To some extent, subregional RSM with the lower R-value could not be produced.

In addition, it was preferable to set NDVI_ATI larger than or equal to NDVI₀ and NDVI_TVDI larger than or equal to NDVI_ATI, namely NDVI₀ ≤ NDVI_ATI ≤ NDVI_TVDI. In that case, the TVDI value was given to all pixels in the TVDI subregion when NDVI₀ was less than or equal to NDVI_TVDI. NDVI_ATI less than or equal to NDVI_TVDI could prevent certain overlapping between the ATI subregion and the TVDI subregion. If NDVI₀ was bigger than NDVI_ATI, there would be some pixels filled with ATI value rather than the average of ATI and TVDI value in the ATI/TVDI subregion. As mentioned before, both NDVI₀ and NDVI_ATI range from 0 to 0.5, and [0, 0.7] is for NDVI_TVDI with an interval of 0.01. The total number of NDVI thresholds collocation for one period is 48,620, which means a maximum of 48,620 iterations performed in three subregions when the number of available RSM observation stations was greater than 20.

3.2.6. Combining Subregional RSM Maps to Generate overall RSM Maps

The subregional RSM was separately estimated using the corresponding model with selected optimal NDVI thresholds. Then, we combined them to obtain an overall RSM map for one period. The combined equations were used as follows:

$$RS{M}_{\mathit{overall}}=\left\{\begin{array}{c}RS{M}_{ATI}=a_{ATI} \times ATI+b_{ATI} NDVI\in \left[0,NDV{I}_{ATI}\right]\\ RS{M}_{ATI/TVDI}=a_{ATI/TVDI} \times \frac{ATI+TVDI}{2}+b_{ATI/TVDI} NDVI\in \left(NDV{I}_{ATI},NDV{I}_{TVDI}\right)\\ RS{M}_{TVDI}=a_{TVDI} \times TVDI+b_{TVDI} NDVI\in \left[NDV{I}_{TVDI}\mathrm{,1}\right]\end{array}\right.$$

(12)

where $RS{M}_{overall}$ represents RSM obtained with the optimal NDVI thresholds using the combined method, $RS{M}_{ATI}$ and $RS{M}_{TVDI}$ are the RSM estimated by the ATI-based and TVDI-based models, respectively, $RS{M}_{ATI/TDVI}$ represents 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 measurements in the ATI subregion, $a_{TVDI}$ and $b_{TVDI}$ are coefficients from fitting the TVDI values and in situ RSM measurements 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 measurements in the ATI/TVDI subregion, and $NDV{I}_{ATI}$ and $NDV{I}_{TVDI}$ are the selected optimal thresholds for generating three subregions.

To illustrate the improvement in the proposed method, ATI-based and TVDI-based models were solely applied for RSM estimation. Quantitative comparison between observed and estimated RSM was performed at some stations (even distributed over the CLP) on DOYs 17 (in winter), 97 (spring), 137 (spring), 161 (summer), 209 (summer), 297 (autumn), 329 (autumn) and 353 (winter).

4. Results

4.1. Evaluation of the Calibration and Validation Results

4.1.1. Comparing Validation Results

While 46 periods of data were analyzed in 2017 (DOYs 1-361 with an 8-day interval), only 38 RSM maps (except for DOYs 1, 9, 33, 65, 73, 201, 225, and 241) were obtained eventually using the above-mentioned models over the study area. For the CLP, winter (December-February) was the coldest season of 2017 (Julian day/DOY337-361 and DOY 1-49). While DOY 57 to 145 and DOY 153 to 233 were generally considered as spring and summer, respectively, the rest of the days of the year belonged to autumn (DOY 241-329). It should be noted that 74-213 stations in the CLP were used and the number of stations (~100) used in winter was less than the other periods (usually exceeding 200) (

Table 3. Calibration and validation results including selected optimal NDVI thresholds and wet/dry edge fitting R².

DOY	Month/ Season	Calibration							Validation
		No. of Stations	Model	$\mathbf{N}\mathbf{D}\mathbf{V}{\mathbf{I}}_0$	$\mathbf{N}\mathbf{D}\mathbf{V}{\mathbf{I}}_{\mathbf{A}\mathbf{T}\mathbf{I}}$	$\mathbf{N}\mathbf{D}\mathbf{V}{\mathbf{I}}_{\mathbf{T}\mathbf{V}\mathbf{D}\mathbf{I}}$	${\mathbf{R}}_{\mathbf{d}\mathbf{r}\mathbf{y}}^2$	${\mathbf{R}}_{\mathbf{w}\mathbf{e}\mathbf{t}}^2$	$\mathbf{R}\mathbf{M}\mathbf{S}\mathbf{E} \pm \mathbf{S}\mathbf{T}\mathbf{D} \left(\mathbf{\%}\right)$	$\mathbf{M}\mathbf{A}\mathbf{E} \pm \mathbf{S}\mathbf{T}\mathbf{D}$	$\mathbf{R} \pm \mathbf{S}\mathbf{T}\mathbf{D}$
17	Jan/Winter	75	ATI/TVDI	0.03	0.15	0.21	0.51	0.67	5.06 ± 0.066	0.11 ± 0.038	0.57 ± 0.011
25	Jan/Winter	74	ATI	0.05	0.14	0.55	0.70	0.71	4.63 ± 0.067	0.03 ± 0.018	0.50 ± 0.018
41	Feb/Winter	86	ATI	0.00	0.15	0.22	0.17	0.59	4.47 ± 0.090	0.12 ± 0.068	0.36 ± 0.030
49	Feb/Winter	101	ATI/TVDI	0.06	0.14	0.64	0.62	0.45	5.23 ± 0.071	0.05 ± 0.041	0.38 ± 0.023
57	Mar/Spring	157	ATI/TVDI	0.09	0.19	0.26	0.92	0.57	5.46 ± 0.091	0.07 ± 0.036	0.36 ± 0.030
81	Mar/Spring	209	ATI/TVDI	0.10	0.18	0.21	0.94	0.92	3.07 ± 0.067	0.05 ± 0.031	0.52 ± 0.024
89	Apr/Spring	209	ATI/TVDI	0.14	0.17	0.22	0.81	0.80	3.63 ± 0.037	0.02 ± 0.014	0.23 ± 0.030
97	Apr/Spring	210	ATI/TVDI	0.11	0.21	0.31	0.90	0.92	3.75 ± 0.068	0.06 ± 0.032	0.20 ± 0.051
105	Apr/Spring	210	ATI/TVDI	0.12	0.21	0.35	0.92	0.86	4.59 ± 0.044	0.03 ± 0.018	0.25 ± 0.028
113	Apr/Spring	210	ATI/TVDI	0.15	0.18	0..21	0.94	0.92	3.70 ± 0.068	0.04 ± 0.022	0.48 ± 0.025
121	May/Spring	210	ATI/TVDI	0.17	0.22	0.27	0.97	0.83	4.43 ± 0.058	0.03 ± 0.022	0.44 ± 0.018
129	May/Spring	211	ATI/TVDI	0.33	0.50	0.70	0.81	0.73	4.52 ± 0.063	0.09 ± 0.047	0.62 ± 0.013
137	May/Spring	211	ATI/TVDI	0.26	0.36	0.52	0.87	0.04	3.43 ± 0.071	0.05 ± 0.025	0.73 ± 0.011
145	May/Spring	211	ATI/TVDI	0.20	0.38	0.57	0.85	0.83	4.62 ± 0.074	0.11 ± 0.060	0.44 ± 0.022
153	Jun/Summer	211	ATI	0.23	0.34	0.41	0.82	0.79	5.56 ± 0.042	0.02 ± 0.018	0.35 ± 0.016
			ATI/TVDI						4.42 ± 0.070	0.12 ± 0.056	0.50 ± 0.020
161	Jun/Summer	208	TVDI	0.29	0.49	0.69	0.76	0.55	5.59 ± 0.187	0.16 ± 0.085	0.34 ± 0.044
			ATI/TVDI						5.39 ± 0.077	0.14 ± 0.048	0.23 ± 0.044
			ATI						5.22 ± 0.035	0.01 ± 0.009	0.54 ± 0.008
169	Jun/Summer	211	ATI/TVDI	0.09	0.36	0.43	0.95	0.59	4.08 ± 0.060	0.07 ± 0.034	0.61 ± 0.015
177	Jun/Summer	211	ATI/TVDI	0.07	0.33	0.46	0.93	0.66	4.44 ± 0.039	0.02 ± 0.016	0.50 ± 0.013
185	Jul/Summer	213	ATI/TVDI	0.33	0.45	0.57	0.90	0.38	4.68 ± 0.038	0.03 ± 0.017	0.63 ± 0.006
193	Jul/Summer	212	ATI/TVDI	0.43	0.44	0.60	0.87	0.28	4.92 ± 0.065	0.08 ± 0.042	0.18 ± 0.041
209	Jul/Summer	212	ATI/TVDI	0.31	0.45	0.56	0.89	0.68	6.01 ± 0.127	0.05 ± 0.038	0.58 ± 0.020
217	Aug/Summer	213	ATI/TVDI	0.00	0.32	0.51	0.72	0.68	5.10 ± 0.094	0.04 ± 0.028	0.57 ± 0.020
233	Aug/Summer	213	ATI/TVDI	0.29	0.29	0.56	0.80	0.60	5.57 ± 0.101	0.08 ± 0.034	0.21 ± 0.036
249	Sep/Autumn	213	ATI/TVDI	0.24	0.36	0.50	0.85	0.82	6.02 ± 0.105	0.18 ± 0.062	0.37 ± 0.040
			ATI						5.06 ± 0.135	0.12 ± 0.098	0.46 ± 0.049
257	Sep/Autumn	213	ATI/TVDI	0.25	0.50	0.55	0.84	0.82	4.41 ± 0.053	0.16 ± 0.035	0.21 ± 0.021
			ATI						5.59 ± 0.050	0.03 ± 0.030	0.32 ± 0.022
265	Sep/Autumn	212	ATI/TVDI	0.11	0.20	0.31	0.81	0.86	4.99 ± 0.130	0.23 ± 0.077	0.50 ± 0.041
273	Oct/Autumn	212	ATI/TVDI	0.11	0.25	0.33	0.93	0.76	6.79 ± 0.212	0.17 ± 0.130	0.17 ± 0.060
			ATI						4.00 ± 0.034	0.03 ± 0.020	0.43 ± 0.013
281	Oct/Autumn	212	ATI/TVDI	0.28	0.31	0.38	0.58	0.75	4.79 ± 0.130	0.38 ± 0.071	0.47 ± 0.046
289	Oct/Autumn	212	ATI/TVDI	0.00	0.26	0.38	0.74	0.88	5.42 ± 0.024	0.04 ± 0.036	0.32 ± 0.018
297	Oct/Autumn	212	ATI/TVDI	0.25	0.25	0..33	0.82	0.75	3.31 ± 0.066	0.03 ± 0.027	0.67 ± 0.014
305	Nov/Autumn	212	ATI/TVDI	0.02	0.17	0.26	0.51	0.84	3.95 ± 0.036	0.02 ± 0.025	0.45 ± 0.014
313	Nov/Autumn	213	ATI/TVDI	0.00	0.22	0.25	0.51	0.77	2.83 ± 0.039	0.05 ± 0.031	0.75 ± 0.008
321	Nov/Autumn	189	ATI/TVDI	0.20	0.20	0.23	0.72	0.53	3.84 ± 0.046	0.05 ± 0.034	0.47 ± 0.014
329	Nov/Autumn	185	ATI/TVDI	0.17	0.23	0.31	0.70	0.60	4.35 ± 0.054	0.19 ± 0.051	0.67 ± 0.008
337	Dec/Winter	185	ATI/TVDI	0.00	0.23	0.31	0.40	0.36	4.72 ± 0.066	0.25 ± 0.056	0.64 ± 0.010
345	Dec/Winter	185	ATI/TVDI	0.08	0.20	0.32	0.68	0.65	4.60 ± 0.069	0.04 ± 0.017	0.54 ± 0.017
353	Dec/Winter	185	ATI/TVDI	0.03	0.19	0.22	0.45	0.56	4.54 ± 0.045	0.03 ± 0.021	0.60 ± 0.010
			ATI/TVDI	0.03	0.25	0.35			3.86 ± 0.026	0.04 ± 0.034	0.61 ± 0.007
361	Dec/Winter	182	TVDI	0.04	0.16	0.18	0.82	0.84	5.09 ± 0.041	0.02 ± 0.012	0.32 ± 0.018
			ATI/TVDI						4.89 ± 0.166	0.50 ± 0.066	0.52 ± 0.025

). Fewer measurements were available for the study area due to RSM not being measured when the soil was frozen. The RSM maps on DOYs 1, 9, 33, 65, 73 were uncalculated in the cold days since the average R-value was relatively low (less than 0.18).

The last column of Table 3 has the highest R values in the validation corresponding to optimal NDVI thresholds—after 10 rounds of 10-fold cross-calibration for each 8-day data. As the DOY increased, the value of R fluctuated over time with no clear trends, while the R of DOY 313 imagery was the largest, up to 0.75 ± 0.008, using the ATI/TVDI joint model, and the minimum R was 0.17 ± 0.060, using the ATI/TVDI joint model, on DOY 273. In addition, the error of the model was calculated with the largest and smallest RMSE values being 6.79 ± 0.212% (DOY 273) and 2.83 ± 0.039% (DOY 313), both using the ATI/TVDI joint model. In terms of MAE, the maximal MAE was 0.50 ± 0.066 using the ATI/TVDI joint model on DOY 361 and the minimum MAE 0.01 ± 0.009 using the ATI-based model on DOY 161.

4.1.2. The NDVI_ATI and NDVI_TVDI Thresholds for Generating Subregions

As shown in Table 3, no fixed NDVI_ATI and NDVI_TVDI thresholds were applying to data throughout the year. In this study, it is observed that most of the identified NDVI_ATI thresholds were in the vicinity of 0.17 in the first four months and in November and December. However, NDVI_ATI in almost the entire summer (DOYs 153-217) was larger (more than 0.3) than in the other seasons. It is noted that only a few periods, e.g., DOYs 129, 161, and 257, saw their NDVI_ATI value approximated the upper limit (0.5) that was set. Table 3 also shows that the smallest NDVI_ATI (0.14) was used for DOY 25 and 49.

Regarding the NDVI_TVDI thresholds over the study area, the months from May to September witnessed larger NDVI_TVDI values (~0.5). Additionally, the minimum NDVI_TVDI appeared in the later months of 2017 when NDVI_ATI was quite small as well. As for modeling, the joint index M_ATI/TVDI was used for nearly all periods (36/38). In other words, most of the RSM was estimated using the ATI/TVDI joint model in the ATI/TVDI subregion. The value of R was relatively low (0.34 ± 0.044 on DOY 161 and 0.32 ± 0.018 on DOY 361, respectively) solely using the TVDI-based model, which is in agreement with the results from Lu et al. [17]. The subregional RSM maps were merged for DOY 257 and DOY 361 (Figure 7 and Figure 8). There are some discrepancies at the boundary of the subregions for the imagery of two periods because they are constructed according to different NDVI ranges—namely the selected optimal NDVI thresholds.

4.1.3. The Threshold NDVI₀ for TVDI

Our results demonstrate that NDVI₀ was generally small, particularly in the first four months and the last month of the year. It was 0, at its lowest, on DOYs 41, 217, 289, 313, and 337, rising to 0.43 on DOY 193. The NDVI and LST on the dry edge and wet edge generally showed significant correlations with ${\mathrm{R}}_{dry}^2$ and ${\mathrm{R}}_{wet}^2$ up to 0.97 on DOY 121 and 0.92 on DOYs 81 and 113, respectively (Table 3).

The parameters of dry and wet edge equations (slope a, intercept b, and correlation coefficient R) are illustrated in Figure 9 for a better understanding of the fitting results. These observations were consistent with the TVDI concept—i.e., the dry edge with a negative slope and the wet edge with a positive slope formed the NDVI-LST triangular space for all the periods of 2017 (Figure 9). Despite similar patterns, the slope of dry edge (a_dry) had higher values in cold days than on warm days, which was the opposite in the case of the slope of wet edge (a_wet). The intercepts b_dry and b_wet also had similar trends (values were higher on warm days than on cold days) but the trend of b_dry was more prominent. Regarding the correlation coefficients, ${ \mathrm{R}}_{dry}^2$ was generally higher than ${ \mathrm{R}}_{wet}^2$.

4.2. Analysis of RSM Estimation Using an Optimal Model and Optimal NDVI Thresholds

4.2.1. Calibration and Validation Analysis of RSM Estimation

The subregional RSM was estimated by applying the optimal thresholds (NDVI₀, NDVI_ATI, and NDVI_TVDI) for each period. In practice, no RSM estimation was performed for subregions with low R-values. According to

Table A3. Data description and accuracy of these models for RSM retrieval.

DOY	Calibration			Validation
	Model	Fitting Equation	R2	R2	RMSE	MAE
17	ATI/TVDI	y = 70.61x − 6.39	0.36 ***	1.00	0.91	0.90
25	ATI	y = −239.74x + 19.98	0.28 *	1.00	1.89	1.36
41	ATI	y = 259.54x + 2.732	0.23 **	1.00	0.99	0.74
49	ATI/TVDI	y = 32.06x + 1.64	0.20 **	0.57	3.50	3.39
57	ATI/TVDI	y = 56.46x − 2.16	0.19 **	0.97	1.00	0.97
81	ATI/TVDI	y = 35.11x + 6.41	0.38 **	1.00	1.95	1.82
89	ATI/TVDI	y = 30.24x + 7.22	0.12 **	0.43	2.96	2.49
97	ATI/TVDI	y = 20.34x + 8.26	0.30 ***	0.48	3.57	3.56
105	ATI/TVDI	y = 36.52x + 8.20	0.10 **	0.68	1.37	1.18
113	ATI/TVDI	y = 75.22x + 3.21	0.33 **	0.67	2.11	1.38
121	ATI/TVDI	y = 71.01x − 0.785	0.25 **	0.55	2.09	1.93
129	ATI/TVDI	y = 49.47x − 2.46	0.43 **	1.00	3.10	2.44
137	ATI/TVDI	y = −75.55x + 22.47	0.66 ***	0.72	3.63	2.71
145	ATI/TVDI	y = −92.69x + 27.00	0.30 **	0.93	2.81	2.48
153	ATI	y = −303.76x + 25.47	0.12 ***	0.86	2.21	2.08
	ATI/TVDI	y = 70.05x + 2.73	0.22 **	1.00	1.54	1.30
	TVDI	y = −26.78x + 21.12	0.24 **	0.48	1.57	1.21
161	ATI/TVDI	y = −237.90x + 24.39	0.28 ***	0.73	4.82	4.06
	ATI	y = 37.57x + 3.66	0.19 *	0.92	5.25	3.82
169	ATI/TVDI	y = 53.34x + 5.72	0.44 **	1.00	0.30	0.27
177	ATI/TVDI	y = 57.34x + 1.22	0.29 **	0.48	2.12	1.60
185	ATI/TVDI	y = 67.79x − 1.69	0.42 ***	0.99	3.47	2.88
193	ATI/TVDI	y = 38.35x + 1.11	0.11 *	0.33	5.78	4.75
209	ATI/TVDI	y = 92.55x − 8.24	0.40 ***	0.99	2.48	1.95
217	ATI/TVDI	y = 88.95x − 1.72	0.37 ***	0.63	3.78	3.75
233	ATI/TVDI	y = 40.50x + 4.53	0.17 *	0.60	3.39	2.61
249	ATI/TVDI	y = −235.66x + 27.89	0.32 ***	0.39	2.77	2.36
	ATI	y = 54.38x + 5.53	0.23 **	0.97	1.89	1.56
257	ATI/TVDI	y = −159.38x + 24.87	0.15 ***	0.46	4.44	3.61
	ATI	y = 64.42x + 1.02	0.39 ***	0.40	7.49	6.40
265	ATI/TVDI	y = 58.62x + 5.45	0.42 ***	0.58	2.97	2.73
273	ATI/TVDI	y = −38.13x + 22.98	0.23 ***	0.90	2.58	2.37
	ATI	y = 63.82x + 2.86	0.27 **	0.28	6.75	5.73
281	ATI/TVDI	y = −57.46x + 32.93	0.46 ***	0.66	1.85	1.81
289	ATI/TVDI	y = −35.11x + 25.12	0.15 **	0.78	3.85	3.47
297	ATI/TVDI	y = 85.81x-8.40	0.47 ***	0.98	0.89	0.81
305	ATI/TVDI	y = 53.85x + 4.76	0.23 ***	0.70	2.79	2.40
313	ATI/TVDI	y = 21.85x + 7.78	0.15 *	0.84	2.68	2.19
321	ATI/TVDI	y = 55.92x − 1.92	0.26 ***	0.97	1.12	1.03
329	ATI/TVDI	y = 70.01x − 7.01	0.48 ***	1.00	1.88	1.31
337	ATI/TVDI	y = 57.54x − 4.17	0.47 ***	0.86	1.94	1.55
345	ATI/TVDI	y = 56.01x − 5.73	0.32 ***	0.56	2.36	1.87
353	ATI/TVDI	y = 49.16x − 4.13	0.43 ***	0.90	2.05	2.04
	ATI/TVDI	y = 62.46x − 5.29	0.41 ***	0.81	1.89	1.80
361	TVDI	y = 22.66x − 1.25	0.14 ***	0.53	2.91	2.43
	ATI/TVDI	y = 91.96x − 14.32	0.36 ***	0.95	2.41	2.02

*, **, *** significant at p < 0.1, p < 0.05 and p < 0.01, respectively.

, in the process of calibration, the linear regression between three indices (ATI, M_ATI/TVDI, and TVDI) and in situ RSM measurements was various with both positive and negative relationships in the certain subregion.

For example, M_ATI/TVDI showed a positive correlation with in situ RSM measurements on most of the periods (~80%) except for DOY 137, 145, 161, 249, 257, 273, 281, and 289. In our research, two negative slope values occurred in the fitting equation between ATI and in situ RSM measurements on DOY 25 and 153. TVDI was positively correlated with in situ RSM measurements on DOYs 161 and 361 with the regression coefficient 0.24 and 0.14, respectively, in the calibration procedure. The largest R² was 0.66 and 1 in the calibration and validation process, respectively. It is worth mentioning that the fitting outcomes for each period passed the hypothesis test with confidence of 0.1 (Table A3).

4.2.2. Estimated RSM Maps over the CLP

The spatial distribution of RSM is illustrated by season (Figure 10, Figure 11, Figure 12 and Figure 13). It is clear that only the RSM in the pixels whose NDVI was within the specified NDVI range could be estimated for each imagery. The used subregions could be marked in the upper left corner of each subfigure. While NDVI_ATI < NDVI < NDVI_TVDI means that the joint index M_ATI/TVDI was used for the ATI/TVDI subregion, NDVI ≤ NDVI_ATI, and NDVI_TVDI ≥ NDVI represent that indices ATI and TVDI were used for the ATI subregion and TVDI subregion, respectively.

The estimated RSM tended to cluster in the northwestern of the study area when the ATI-based model was applied (Figure 10b,c, Figure 12a and Figure 13a,b,d). These regions are mainly covered by sparse vegetation or bare soil because of lower NDVI. However, the TVDI-based model could estimate RSM over the southern and eastern parts (Figure 10h).

From the perspective of RSM seasonal variation, RSM in red color was low, particularly in the southern part (RSM less than 5%) (Figure 10b,d–h) and RSM in green color was relatively high in northern China in winter (Figure 10a,e,f,h). The overall RSM in green and light blue color gradually increased with a rise to 40% in the west area in spring (Figure 11d and Figure 10e,f,g). Nearly 50% of the RSM in ultra-blue color was observed in the west area in summer (Figure 12c and Figure 10d–h). RSM with the yellow and light green color in most areas increased to its peak in August and September (Figure 12h,i and Figure 13a–d). Then RSM in orange and red color dropped gradually in the following days (Figure 13g,i–k).

5. Discussion

5.1. Applicability of the ATI/TVDI Joint Model

In this paper, the indices ATI and TVDI were used for modeling and estimating RSM in the CLP during the whole year of 2017. NDVI is a vital factor influencing not only the wet and dry edges for obtaining TVDI but also the application of models in the subregions. Compared to the ATI-based and TVDI-based models, the ATI/TVDI joint model achieved higher accuracy in calibration and the ATI-based and TVDI-based models’ advantages had a good performance in RSM estimation. It is observed from Table A3 that both ATI-based and TVDI–based models are seldom used for RSM estimation, which was consistent with the observations reported by Lu [17]. Du et al. performed the original TVDI and the modified TVDI_m for drought monitoring in semi-arid regions of China and found that the improved TVDI_m exhibits better performance in both meteorological and agricultural drought monitoring compared to the original TVDI. The R² between TVDI_m and in situ RSM was 0.21, which is lower than most of the R² in calibration (Table A3) [73]. Zhang et al. calculated the correlation coefficients between TVDI and SM at depths of 10 cm and 20 cm in Southwest China and the maximum R² was 0.29 on DOY 281, which is lower than our maximum R² (0.66) on DOY 137 and 0.46 on DOY 281 [47]. NDVI₀ that we applied played an important role in RSM estimation by the TVDI-based model according to the selected optimal NDVI thresholds. Comparative analysis indicated that the ATI/TVDI joint model had higher applicability (accounting for 36/38 periods) than the ATI-based and TVDI-based models. As shown in Table A3, there is no significant regularity for the parameters in the fitting equations since the studied period is short.

Figure 14 and

Table A4. Comparison of the estimated RSM and R-value in calibration using different models at some stations on DOYs 17, 353 (in winter), 97, 137 (spring), 161, 209 (summer), 297, and 329 (autumn).

Month/Season	DOY	Station	Observed RSM (%)	Combined Model		ATI-Based Model		TVDI-Based Model
				Estimated RSM (%)	R-Value in Calibration	Estimated RSM (%)	R-Value in Calibration	Estimated RSM (%)	R-Value in Calibration
Jan/Winter	17	53,735	0.80	5.99	0.60	9.51	0.37	6.96	0.45
		53,931	5.85	6.89		10.99		7.90
		57,055	13.43	12.36		10.32		11.49
		53,978	16.15	11.16		10.76		10.77
		53,934	3.16	3.60		11.20		5.71
Apr/Spring	97	53,543	7.22	10.58	0.55	15.03	0.12	14.54	0.11
		53,651	12.82	14.65		14.66		15.15
		53,681	18.08	17.45		14.29		15.79
		53,754	12.18	12.18		14.75		15.22
		53,959	10.68	12.00		15.10		15.20
May/Spring	137	53,433	5.90	6.77	0.81	10.83	0.09	10.62	0.08
		53,469	7.50	7.00		10.93		10.63
		53,732	5.62	4.01		11.22		10.33
		53,863	8.21	8.11		11.18		10.70
		53,678	8.29	6.85		11.17		10.59
Jun/Summer	161	53,915	11.72	13.11	0.49/0.53/0.44	14.17	0.43	14.10	0.33
		53,585	8.74	11.36		13.17		11.75
		53,659	10.35	11.72		13.67		12.93
		52,765	29.38	15.42		15.02		7.30
		57,080	18.23	17.78		15.88		14.91
Jul/Summer	209	53,433	3.67	6.19	0.63	14.51	0.09	12.40	0.25
		52,765	33.95	23.62		13.85		17.46
		53,754	16.21	16.23		14.92		14.84
		57,131	13.50	13.75		14.37		14.51
		57,071	13.07	13.69		14.22		14.58
Oct/Autumn	297	53,748	11.77	12.98	0.69	17.20	0.28	16.64	0.27
		57,052	14.20	14.11		18.03		17.21
		53,872	21.58	20.95		13.83		18.73
		53,676	11.78	12.25		16.72		16.28
		53,776	13.82	15.21		16.53		17.28
Nov/Autumn	329	53,512	4.00	7.22	0.69	12.59	0.32	8.10	0.59
		53,553	3.53	5.70		11.57		6.00
		53,585	2.99	6.52		12.44		7.37
		53,780	9.14	10.66		12.18		10.83
		53,868	16.50	15.26		11.95		14.71
Dec/Winter	353	53,478	2.72	0.99	0.66/0.64	8.82	0.30	0.59	0.49
		53,585	2.09	6.35		9.52		7.99
		53,651	3.40	7.47		9.09		8.67
		57,053	14.30	13.29		9.69		12.15
		57,071	16.94	14.41		9.83		13.13

demonstrate that the combined model we applied for RSM retrieval shows the improvements solely over ATI- and TVDI-based models with the minimum difference between the estimated RSM and the observed RSM and maximum R-value (correlation coefficient) in calibration. The comparative analysis indicated that the RSM estimated by the TVDI-based model was closest to the actual RSM measured on-site compared to the ATI-based model, which shows a good agreement with Hsu’s research results [18]. Jiao et al. investigated the relationships between SM and vegetation coverage and suggested that excessive reliance on afforestation was risky for the improvement of vegetation coverage in arid and semi-arid regions and we, thus, need to focus on the soil moisture conditions during large-scale re-vegetation. Thus, the methodology presented in this study has clearly improved the estimation of SM throughout the year of 2017 over the CLP compared to other methods. It highlights the importance of vegetation coverage in soil moisture estimation and the possibility of producing a potentially complete SM estimation map for such a large area by merging subregional maps. Some researchers just performed SM estimation on a small scale at a specific period using the TVDI-based model [38,47,71]. To further this work, more research is required to examine the correlation between other indices and climate regimes as well as using a longer time series MODIS and in situ RSM measurements when they become available [36,82]. It should be noted that there were no in situ RSM measurements in the Ningxia province (Figure 1) and this blank area should be filled when station measurements become available in the future.

In the process of calibration, the 10-fold cross-calibration method was rarely applied to linear regression analysis between satellite-based indices and in situ RSM measurements, and this method indeed improved the accuracy of RSM estimation in terms of validation outcomes with limited in situ measurements. In short, the procedure from identifying NDVI thresholds to modeling could be widely applied because it enables determining an optimal model for RSM estimation.

5.2. Estimation of RSM in the CLP

To show the effectiveness of the selected optimal NDVI thresholds, RSM maps at an 8-day interval were produced for the entire CLP. The limited number of in situ RSM observation stations in the frozen period led to difficulties in ensuring better model calibrations. Regions with the frozen ground or snow-covered are generally masked out [83]. Hence, the average correlation coefficient (R) for verification was generally low at that time. In general, the soil thawing–freezing process occurs from November to March in this area. RSM rises dramatically during the thawing period and decreases rapidly in the freezing period and this phenomenon was not evident in our results. Liu et al. claimed that the TVDI-based method is suitable for obtaining satisfactory outcomes if the effects of both altitude and frozen soil are considered [34]. The main reason for the sudden rise in RSM in the autumn and summer images was very likely to be regional increases in precipitation in the rainy season [84,85,86,87]—e.g., in the west area (Qinghai province) from DOY 241 to 249 and in the northwest area from DOY 265 to 273.

In terms of incomplete RSM maps, we have here some explanations that the overall map was generated by a combination of subregional RSM for one period and it failed to select the highest R-value in certain subregions resulting in no subregional RSM map produced. The reason for that might be the limited availability in situ RSM measurements in those subregions.

5.3. Innovations and Limitations

Studies on RSM estimations and spatiotemporal dynamics on regional or continental scales using MODIS data are important for understanding the overall pattern and developing retrieval model [88]. In our study, an iteration procedure that allows us to identify optimal NDVI thresholds and models for estimating RSM with 10 rounds of 10-fold cross-calibration and validation. We separately used three models (the ATI-based model, the TVDI-based model, and the ATI/TVDI joint model) based on optimal NDVI thresholds in the corresponding subregions for retrieving the 8-day RSM from MODIS imagery with in situ RSM measurements. The combination of different subregional RSM maps for the overall map at large scales, as demonstrated in our case study, provides a new perspective for accuracy RSM estimation.

However, we here would like to highlight some of the weaknesses in this study that can be improved in future work. The MODIS images and in situ RSM measurements from observation stations did not spatially match very well and comparing the RSM estimated from 8-day composite MODIS data with the in situ RSM measurements might decrease the modeling accuracy. In addition, some environmental factors, such as soil type, latitude, topography, and hydrology, can affect the accuracy of RSM estimation and should be considered to improve the retrieval model.

6. Conclusions

In this study, we proposed an iteration procedure that allows us to identify optimal NDVI thresholds by performing 10 rounds of 10-fold cross-calibration and thus to estimate the region-specific relative soil moisture of the Chinese Loess Plateau in 2017 using different models—namely the ATI-based model, the TVDI-based model, and the ATI/TVDI joint model. The key findings and main conclusions are as follows:

• The whole study area was divided into three subregions by selected optimal NDVI thresholds and the overall RSM map was produced by combining subregional RSM maps over the CLP. This combination-based approach provides a new perspective for the mapping of RSM over large geographical areas.
• Spatiotemporal and comparative analysis indicated that the ATI/TVDI joint model had higher applicability (accounting for 36/38 periods) and accuracy than the ATI-based and TVDI-based models. The ATI-based model is only suitable for bare soil and sparsely vegetated areas while the TVDI-based model is more applicable to moderately or highly vegetated areas.
• The 10-fold cross-calibration method in the iteration procedure improves overall reliability and effectiveness in terms of identifying optimal NDVI thresholds.

Our study has innovated the mapping of RSM by applying different estimation models to different land covers; however, it remains necessary to validate the ATI/TVDI joint method with more measurements from various environments. The findings of this study will contribute to decision-makers in agricultural management and formulating policies of the Chinese Loess Plateau.