An Improved Soil Moisture Downscaling Method Based on Soil Properties and Geographical Divisions over the Loess Plateau

Abstract

As the contradiction between vegetation growth and soil moisture (SM) demand in arid zones gradually expands, accurately obtaining SM data is crucial for ecological construction. Remote sensing products limit small-scale studies due to the low resolution, and the emergence of downscaling solves this problem. This study proposes an improved semi-physical SM downscaling method. The effects of environmental factors on SM in different geographical zones (Windy Sand Hills, Flood Plains, Loess Yuan, Hilly Loess, Earth-rock Hills and Rocky Mountain) were analyzed using Random Forests. Vegetation and topographic factors were incorporated into the traditional downscaling algorithm based on the Mualem–van Genuchten model by setting weights, yielding 250 m resolution SM data for the Loess Plateau. This study found the following: (1) The Normalized Difference Vegetation Index (NDVI) was the most important environmental factor in all divisions except the Flood Plain, and the Digital Elevation Model (DEM) was second only to the NDVI in the overall importance evaluation, both of which positively influenced SM. (2) SM variability increased and then decreased when SM was below 0.4 cm³/cm³, but showed a quadratic growth trend when exceeding this threshold. The Rocky Mountain division exhibited the highest variability under the same SM. (3) Validation showed that the improved algorithm, based on geographic divisions to analyze factors importance and interpolation of coarse-scale SM and variability, had the highest accuracy, with an average R of 0.753 and an average ubRMSE of 0.042 cm³/cm³. The improved algorithm produced higher resolution, more accurate SM data, and offered insights for downscaling studies in arid regions, meeting the region’s high-resolution SM needs.

1. Introduction

Arid areas are characterized by scarce precipitation, infertile soils, fragile ecology, and high sensitivity to climate change and anthropogenic activities [1,2,3]. In the context of global warming, the area of arid zones and the degrees of aridity continue to increase [4], and ecological security is facing more serious challenges [5]. Soil moisture (SM), as a prominent part of the water cycle in terrestrial ecosystems, is a major constraint on vegetation growth in arid zones [6]; large-scale vegetation construction has led to frequent negative SM balances [7], and the phenomenon of soil desiccation threatens regional ecological stability and the results of vegetation restoration [8]. Therefore, efficient and accurate acquisition of SM information is crucial for agricultural production guidance [9], drought monitoring [10], and sustainable ecosystem development [11].

The traditional methods used for SM acquisition include the drying method, Time-Domain Reflectometry (TDR) method, Frequency-Domain Reflectometry (FDR) method, Resistive Sensor method, and Digital Tensiometer method [12,13]. Despite the high accuracy of the results [14], they require considerable manpower and material resources, the observation points are sparse and poorly representative, and they are only applicable to SM measurements at the point scale [15]. Large-scale SM data are obtained mainly by inversions of remote sensing (RS) data, which have the advantages of wide coverage, long time series, and low acquisition cost. Based on sensor wavelength, the main inversion technologies can be categorized into optical RS and microwave RS [16]. Optical inversion estimates SM using the reflective characteristics of visible and infrared bands. However, limited by the weak penetration capability of visible light through clouds and surface vegetation, and the indirect correlation between optical data and SM, the inversion accuracy is constrained. In contrast, microwave RS not only possesses certain cloud penetration capabilities, but also demonstrates superior accuracy in inversion algorithms based on backscattering coefficients and brightness temperature, making it the mainstream technology for large-scale SM observation, although it remains limited to surface SM [17]. With the development of land surface process models, large-scale multilayer SM data with high temporal resolution can be obtained with the help of surface data assimilation systems and reanalysis information [18]. Nevertheless, the spatial resolution of existing microwave and reanalysis SM products remains insufficient (GLDAS 0.25°, ERA5-Land 0.1°) [19] to meet the refined requirements of medium- and small-scale hydrological studies. Therefore, utilizing surface factors to construct quantitative downscaling relationship models and applying them at high-spatial-resolution scales, while maintaining the accuracy of SM products and enhancing their spatial resolution, has become a popular subject of research.

In recent years, numerous scholars have extensively explored the downscaling techniques for SM products [20,21,22,23,24,25], and the methods can be summarized as spatial interpolation [26,27], statistical regression [28,29], physical modeling [30,31], data assimilation [32,33], and machine learning [34,35,36]. The spatial interpolation and statistical regression methods are commonly used because of their simplicity and easy access to auxiliary data [27,37,38], but the regression equations are not applicable to large-scale regions due to the lack of physical basis [39]; data assimilation and physical modeling techniques are complex in form, and some of the auxiliary data are difficult to obtain through RS [40,41], which prevents them from being widely applied [42]; and the machine learning method has a strong nonlinear fitting ability and can obtain the downscaled results with higher accuracy when the training samples are sufficient [43], but lack interpretability [44]. Although different downscaling methods have their own advantages and disadvantages, these studies have focused on downscaling models and the selection of ancillary data, while neglecting the SM transport mechanism in the downscaling process.

SM presents intense spatial and temporal variability, owing to the interactions among climate, topography, vegetation, and soil properties [45,46,47], and characterizing its spatial variability is crucial for downscaling studies of SM products. Qu et al. [48], based on the Mualem–van Genuchten (MvG) model [49], used soil texture data to stochastically analyze one-dimensional unsaturated gravity flow and derived a closed expression describing the variability at a certain SM. Unlike establishing statistical relationships between SM and auxiliary factors, this expression analyzed the dynamics of SM within a unit grid based on soil texture and hydraulic parameters and provided a new idea for the decomposition of SM products at the coarse scale. Montzka et al. [50] estimated the standard deviation of the Soil Moisture Active Passive (SMAP) coarse-grid using the physical method of Qu to obtain 1 km SM data for the Upper Rhine Valley in Germany, but the auxiliary data chosen were a constant Field Capacity (FC), which could not describe the change in SM in time. Cao et al. [51] replaced auxiliary data with dynamic Apparent Thermal Inertia (ATI) data on this basis, and also obtained 1 km downscaled SMAP SM data in the alpine region of northwestern China.

However, the above studies have the common shortcoming of downscaling SM products by using only a single type of auxiliary data and the SM variability, ignoring the joint effect of multiple environmental factors on the distribution of SM in space. Therefore, this study aims to improve the traditional methods and uses the Loess Plateau, a typical ecologically fragile region in the world with aridity and low rainfall, complex topography, and significant differences in vegetation distribution, as the study area, and selects typical vegetation and topographic factors as auxiliary data. Random Forests (RFs) were used to analyze the impact of environmental factors on SM across geographical divisions and set the weights of each variable. Combined with the standard deviation data for SM on the ERA5-Land grid calculated by the physical model, a downscaling method that integrates a variety of environmental factors was established to downscale the ERA5-Land SM product to 250 m and was validated using the measured data.

2. Materials and Methods

2.1. Study Area and Data

2.1.1. Study Area

The Loess Plateau is situated in the north–central part of China, with a total area of 62.85 × 10⁴ km², and it is in the transition zone from the plain to the plateau [52]. The climate is transitional, with a warm–temperate zone with a semimoist climate in the southeast to the arid–semiarid zone with a medium–temperate zone in the northwest [53], and the main types of landforms are hills, loess, terraces, plains, and tufo-stony mountains, etc. [54]. The average annual temperature ranges from 3.6 to 14.3 °C [55] and decreases with increasing latitude. Annual precipitation ranges from 200 to 800 mm [56], decreasing gradually from southeast to northwest, with uneven spatial distribution and large seasonal changes, mainly concentrated in summer, accounting for 50–65% of the annual precipitation [57]. Based on topography and geology, the region can be divided into six geographical divisions, including Earth-rock Hills, Rocky Mountain, Loess Yuan, Hilly Loess, Windy Sand Hills, and Flood Plains (Figure 1).

2.1.2. ERA5-Land Data

ERA5-Land is the land surface product dataset of ERA5 (the Fifth Generation of European Reanalysis) and contains global data at high spatial resolution (0.1°) and temporal resolution (1 h) for more than 50 variables since 1950 (https://cds.climate.copernicus.eu/, accessed on 12 February 2025). The dataset provides SM data for four layers: 0–7, 7–28, 28–100, and 100–289 cm [58]. A variety of SM products have been evaluated by scholars in typical areas of northern China, and the overall performance of the ERA5 product was the best in terms of site assessment [59]. To match the measured data, this study used 0–7 cm SM data as the research object, used the 2011–2013 vegetation growing season (April–September) as the research period, and hourly data were averaged for each day to obtain daily resolved SM data.

2.1.3. Soil Property Data

In this study, we used the OpenLandMap dataset (https://openlandmap.org/, accessed on 12 February 2025), which provides global high-resolution soil properties, including the soil bulk weight, clay content, sand content, organic matter content, and field capacity at 250 m spatial resolution, in the study area. The soil texture information was converted into soil hydraulic parameters via the soil transfer function that was developed by Wösten [60], which, in turn, was solved for the MvG model parameters (Ks, θs, θr, α, and n) that are used to compute the SM variability in the ERA5-Land grid.

2.1.4. Auxiliary Data

Considering that the Loess Plateau is rich in geomorphic types, gullies and ravines, and diverse vegetation types, the vegetation and topographic factors strongly influence the dynamic changes in and spatial distribution of SM [61,62,63,64]. Therefore, the typical Normalized Difference Vegetation Index (NDVI), Digital Elevation Model (DEM), slope, and aspect indices were selected as auxiliary data, which have been widely used in related studies and have shown important roles [65,66,67,68]. The NDVI data were selected from MODIS for a 16-day synthetic spatial resolution of 250 m for the vegetation index product, MOD13Q1, and the SRTM 90 m resolution elevation data were selected for the DEM data, both of which were provided free of charge by NASA Earth Science Data (https://earthdata.nasa.gov/, accessed on 12 February 2025). The DEM data were aggregated to 250 m resolution in ArcGIS, and slope and aspect data were generated at this resolution. The aspect data were converted from compass values of 0–360° to aspect indexes of 0–1 according to the method of Robert and Cooper [69], and the conversion equation was as follows:

$$\mathrm{T}\mathrm{R}\mathrm{A}\mathrm{S}\mathrm{P}=\left\{1-cos\left[\left(\pi /180\right)\left(aspect-30\right)\right]\right\}/2$$

(1)

where TRASP is the transformation of the aspect, with larger values indicating drier and hotter environments, where 0 represents the north–northeast direction and 1 represents the south–southwest direction.

2.1.5. Site Observation Data

The measured site data are from 1992 to 2013 and are based on data from 732 sites in China in which soil relative humidity and soil characteristic parameter observations were observed together to establish a volumetric SM dataset (https://www.scidb.cn/, accessed on 12 February 2025), a temporal resolution of one month (from the average of the observations on the 8th, 18th, and 28th of each month). SM at each site was measured gravimetrically with strict quality controls and was converted to volumetric SM content by the law of mass–volume relationships [70]. In this study, the 0–10 cm SM observations, which are closest to the ERA5-Land dataset, were selected to validate the downscaled results. Observations with missing data were not used for comparison and evaluation, and 148 stations that had continuous data during the study period were ultimately selected, covering all geographic divisions and climate types of the Loess Plateau.

2.2. Method

2.2.1. ERA5-Land Grid Soil Moisture Variability Calculations

The SM variability corresponding to the coarse grid-averaged SM can be estimated via the MvG model, which describes the soil and water conservation curve equations and hydraulic conductivity equations as follows:

$$S_e={\displaystyle \frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}}={\displaystyle \frac{1}{{(1+{\left(\alpha \left|h\right|\right)}^n)}^m}},\hspace{1em}h<0$$

(2)

$$m=1-{\displaystyle \frac{1}{n}}$$

(3)

$$K\left(S_e\right)=K_s{S}_e^L{\left[1-{\left(1-S_e^{{\displaystyle \frac{1}{m}}}\right)}^m\right]}^2,h<0$$

(4)

where $S_e$ is the effective saturation (-), $\theta$ is the SM (cm³/cm³), ${\theta }_r$ is the residual SM (cm³/cm³), ${\theta }_s$ is the saturated SM (cm³/cm³), $\alpha$ is the inlet value parameter (cm⁻¹),$h$ is the pressure head (cm), $n$(-) is the pore-size distribution parameter, $K_s$ is the saturated permeability coefficient (cm/d), K is the unsaturated permeability coefficient (cm/d), and L is the void channel parameter (L = 0.5).

Qu et al. [48] found the relationship between MvG hydraulic parameters and SM variability with the equation shown below:

$$\begin{array}{c}{\sigma }_{\theta }^2=b_0^2\{b_1^2{\sigma }_{\alpha }^2+b_2^2\left[{\displaystyle \frac{{\sigma }_f^2{\rho }_f}{\left(1+a_2{\rho }_f\right)a_2}}+{\displaystyle \frac{a_1{\sigma }_{\alpha }^2{\rho }_{\alpha }}{\left(1+a_2{\rho }_a\right)a_2}}+{\displaystyle \frac{a_3{\sigma }_n^2{\rho }_n}{\left(1+a_2{\rho }_n\right)a_2}}\right]+b_3^2{\sigma }_n^2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+b_4^2{\sigma }_{{\theta }_s}^2+2b_1{b}_2\left(-{\displaystyle \frac{a_1{\sigma }_{\alpha }^2{\rho }_{\alpha }}{1+a_2{\rho }_a}}\right)+2b_2{b}_3\left(-{\displaystyle \frac{a_3{\sigma }_n^2{\rho }_n}{1+a_2{\rho }_n}}\right)\}\end{array}$$

(5)

where $\rho$ is the vertical correlation length of each parameter and where $a_i$ and $b_i$ are the parameter values obtained through the calculation of each hydraulic parameter, which are calculated as follows:

$$a_1={\displaystyle \frac{\left(\frac{5}{2}-\frac{1}{2\overline{n}}\right){\left(\overline{a}\overline{h}\right)}^{\overline{n}}}{1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}}{\displaystyle \frac{\overline{n}}{\overline{a}}}$$

(6)

$$a_2={\displaystyle \frac{\left(\frac{5}{2}-\frac{1}{2\overline{n}}\right){\left(\overline{a}\overline{h}\right)}^{\overline{n}}}{1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}}{\displaystyle \frac{\overline{n}}{\overline{h}}}$$

(7)

$$a_3={\displaystyle \frac{\left(\frac{5}{2}-\frac{1}{2\overline{n}}\right){\left(\overline{a}\overline{h}\right)}^{\overline{n}}}{1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}}\ln \left(\overline{a}\overline{h}\right)+{\displaystyle \frac{\ln \left(1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}\right)}{2{\overline{n}}^2}}-{\displaystyle \frac{2}{{\overline{n}}^2-\overline{n}}}$$

(8)

$$b_0=\left(\overline{{\theta }_s}-\overline{{\theta }_r}\right)\left\{{\displaystyle \frac{\overline{a}\overline{h}}{\left[1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}\right]{\left(\overline{a}\overline{h}\right)}^{\overline{n}}\overline{n}}}\right\}$$

(9)

$$b_1={\displaystyle \frac{{\overline{n}\left(\overline{a}\overline{h}\right)}^{\overline{n}}+1-\overline{n}}{\overline{a}}}-{\displaystyle \frac{\left[{\overline{n}\left(\overline{a}\overline{h}\right)}^{\overline{n}}+1\right]{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}{1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}}{\displaystyle \frac{\overline{n}}{\overline{a}}}$$

(10)

$$b_2={\displaystyle \frac{{\overline{n}\left(\overline{a}\overline{h}\right)}^{\overline{n}}+1-\overline{n}}{\overline{h}}}-{\displaystyle \frac{\left[{\overline{n}\left(\overline{a}\overline{h}\right)}^{\overline{n}}+1\right]{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}{1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}}{\displaystyle \frac{\overline{n}}{\overline{h}}}$$

(11)

$$b_3=-{\displaystyle \frac{1}{\overline{n}}}-\ln \left(\overline{a}\overline{h}\right)-\ln \left(\overline{a}\overline{h}\right){\displaystyle \frac{\left[{\overline{n}\left(\overline{a}\overline{h}\right)}^{\overline{n}}+1\right]{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}{1+{\left(\overline{a}\overline{h}\right)}^{\overline{n}}}}$$

(12)

$$b_4={\overline{n}\left(\overline{a}\overline{h}\right)}^{\overline{n}}+1$$

(13)

$$f=\ln K_s$$

(14)

The above closed expression establishes the relationships among the mean SM variability, SM content, vertical correlation length of the MvG parameter, mean value of the hydraulic parameter, and standard deviation of the hydraulic parameter.

2.2.2. Downscaling Methods

Based on the original SM data to calculate the SM variability at low resolution, combined with the spatial distribution characteristics of the high-resolution auxiliary data at the scale of 250 m, the SM product is decomposed to obtain high-resolution SM data. The auxiliary data need to be normalized before the calculations, and the downscaling formula is as follows:

$$\widetilde{{\theta }_{i,j}}=\overline{\theta }+{\sigma }_{\theta }\left(\overline{\theta }\right){\displaystyle \frac{P_{i,j}-\overline{P}}{{\sigma }_p}}$$

(15)

where $\widetilde{{\theta }_{i,j}}$ is the high-resolution SM; $\overline{\theta }$ is the ERA5-Land SM at a low resolution; ${\sigma }_{\theta }\left(\overline{\theta }\right)$ is the standard deviation of the SM at the particular SM value; $P_{i,j}$ is the high-resolution auxiliary data; $\overline{P}$ is the average value within a unit image element at a low resolution for the auxiliary data; and ${\sigma }_p$ is the standard deviation of the auxiliary data within a unit image element at a low resolution. In Montzka’s study [50], Field Capacity (FC) was used for auxiliary data. However, FC is relatively stable and does not change significantly within a short period of time, and its values in spatially neighboring areas are usually close to each other, so the spatial detail information of the downscaled SM is relatively limited. In view of the good linear correlation between Apparent Thermal Inertia (ATI) and SM [71,72], Cao et al. [51] used dynamic ATI instead of FC for downscaling. In this study, we attempt to improve this method by adding several types of auxiliary data (NDVI, DEM, SLOPE, TRASP, and FC) to the downscaling process to obtain more accurate SM data.

By identifying the complex linear or nonlinear features that exist between downscaling factors and SM [73], Equation (15) can be decomposed according to the contribution of each variable. Using the relative importance score of auxiliary data, the RF can assess the influence of environmental variables on SM. In this study, %IncMSE was chosen as the evaluation criterion, and a high %IncMSE score for a variable indicates that the value of the variable will significantly increase the prediction error of the model when it is randomly replaced; thus, a higher score indicates that the variable is more important. The weights of each auxiliary data were set according to the %IncMSE scores, and the effects of multiple variables on SM were integrated into the downscaled results. The RF can be used in RStudio, and the improved equation is as follows:

$$\widetilde{{\theta }_{i,j}}=\overline{\theta }+{\sigma }_{\theta }\left(\overline{\theta }\right)\left[a\right.{\displaystyle \frac{A_{i,j}-\overline{A}}{{\sigma }_A}}\left.+b{\displaystyle \frac{B_{i,j}-\overline{B}}{{\sigma }_B}}+c{\displaystyle \frac{C_{i,j}-\overline{C}}{{\sigma }_C}}+d{\displaystyle \frac{D_{i,j}-\overline{D}}{{\sigma }_D}}+e{\displaystyle \frac{E_{i,j}-\overline{E}}{{\sigma }_E}}\right]$$

(16)

where $A_{i,j},B_{i,j},C_{i,j}{,D}_{i,j}{,\mathrm{a}\mathrm{n}\mathrm{d}E}_{i,j}$ are the high-resolution auxiliary data, NDVI, DEM, SLOPE, TRASP, and FC, respectively; $\overline{A},\overline{B},\overline{C},\overline{D},and\overline{E}$ are the mean values of the corresponding auxiliary data per unit image element at the low-resolution; ${\sigma }_A,{\sigma }_B,{\sigma }_C,{\sigma }_D,and{\sigma }_E$ are the standard deviations of the corresponding auxiliary data per unit image element at the low-resolution; and a, b, c, d, and e are the weighting factors, which are positively or negatively determined by the correlations among the environmental factors and the SM.

In order to compare the differences in accuracy and spatial details between the different downscaling methods, this study utilized the traditional method (using only FC as auxiliary data), the traditional interpolation-based method, the improved method incorporating multiple auxiliary data, the improved interpolation-based method, and the further improved method based on interpolation and geographic divisions for downscaling.

2.2.3. Assessment and Validation

In order to match the dates of the SM measurement and calculation methods at the validation station, the downscaled results on the 8th, 18th, and 28th of each month in April–September of 2011, 2012, and 2013 were selected for the mean value calculation to obtain monthly scale SM downscaled data for validation. Typical evaluation metrics, including the correlation coefficient (R), mean bias (Bias), root mean square error (RMSE), and unbiased root mean square error (ubRMSE), were selected to evaluate the accuracy of the traditional and improved SM downscaling algorithms [74]. In addition, statistical gain metrics, including G_PREC and G_RMSE, were used to evaluate the changes in the errors before and after downscaling [75]. The above metrics can be calculated as follows in RStudio:

$$R={\displaystyle \frac{E\left[\left({SM}_s-E\left[{SM}_s\right]\right)·\left({SM}_g-E\left[{SM}_g\right]\right)\right]}{{\sigma }_s·{\sigma }_g}}$$

(17)

$$Bias=E\left[{SM}_s\right]-E\left[{SM}_g\right]$$

(18)

$$RMSE=\sqrt{E\left[{\left({SM}_s-{SM}_g\right)}^2\right]}$$

(19)

$$ubRMSE=\sqrt{E\left[{\left(\left({SM}_s-E\left[{SM}_s\right]\right)-\left({SM}_g-E\left[{SM}_g\right]\right)\right)}^2\right]}$$

(20)

$$G_{PREC}={\displaystyle \frac{\left|1-R_L\right|-\left|1-R_H\right|}{\left|1-R_L\right|+\left|1-R_H\right|}}$$

(21)

$$G_{RMSE}={\displaystyle \frac{{RMSE}_L-{RMSE}_H}{{RMSE}_L+{RMSE}_H}}$$

(22)

where $E\left[·\right]$ represents the mean operator, s represents ERA5-Land product and downscaled results, g represents site-observed SM, and L and H represent low- and high-resolution SM, respectively. The gain index ranges from −1 to 1, with larger values implying more significant accuracy gains after downscaling.

3. Results

3.1. Importance Analysis of the Downscaling Factors

Considering that the effect of environmental variables on SM is closely related to the study period and area [76], and that the Loess Plateau has a large range, including a variety of terrain, soil, and vegetation types, the importance analysis of environmental factors based on geographical divisions could be set to a more regional characterized weight value by using RF’s variable importance measures (

Table 1. Weights of environmental factors for each geographical division.

Environmental Factors		NDVI	DEM	TRASP	SLOPE	FC
Windy Sand Hills	R	0.512 *	0.077 *	−0.107 *	−0.088 *	0.268 *
	weight	0.334	0.157	0.165	0.169	0.175
Flood Plain	R	0.505 *	−0.353 *	−0.103 *	0.349 *	0.491 *
	weight	0.241	0.135	0.184	0.248	0.192
Loess Yuan	R	0.364 *	0.325 *	−0.133 *	0.072 *	−0.148
	weight	0.386	0.167	0.140	0.151	0.156
Hilly Loess	R	0.539 *	0.276 *	−0.165 *	0.383 *	0.098 *
	weight	0.405	0.164	0.125	0.131	0.174
Earth-rock Hills	R	0.322 *	0.058 *	−0.116 *	0.027 *	0.054
	weight	0.336	0.161	0.170	0.162	0.170
Rocky Mountain	R	0.389 *	0.367 *	−0.157 *	0.301 *	0.229 *
	weight	0.289	0.230	0.158	0.142	0.181
Overall	R	0.509 *	0.342 *	−0.135	0.366 *	0.284 *
	weight	0.301	0.235	0.146	0.140	0.178

* p-value < 0.01.

).

During the study period, the ranking results of the downscaling factors importance analysis differed slightly among geographical divisions (Figure 2), but NDVI was the most important impact factor in all factors, suggesting that it dominated the changes in SM. A high NDVI indicates sufficient SM and good vegetation growth. Topographic relief affects the spatial distributions of air temperature, precipitation and vegetation and then indirectly affects the SM. The results of the correlation analysis of the three types of terrain factors in different geographical divisions varied greatly, especially in the Windy Sandy Hills and Earth-rock Hills areas. DEM correlation with SM was very weak, but in the overall analysis of the correlation with SM, it was second only to NDVI. Notably, DEM showed a negative correlation with SM in the Flood Plain, which may be related to the presence of large areas of agricultural land in the region, where the terrain is flat for easy accumulation of precipitation and some of the agricultural land also has irrigation facilities. TRASP index reflects the degrees of heat and dryness at the surface, and different aspects correspond to different hydrothermal conditions and vegetation statuses. Generally, sunny slopes receive more solar radiation, and SM evaporation is stronger than that on shady slopes. Additionally, the FC selected in the traditional downscaling algorithm did not show an outstanding importance or correlation ranking, which suggests that the environmental factors introduced in this study were important for improving the downscaling algorithm.

3.2. Estimation of Soil Moisture Variability

Using the MvG parameters and the mean and variance of each parameter, the standard deviation of SM that express the SM variability in each ERA5-Land grid were estimated with Equation (5). Typically, grids with complex soil textures correspond to larger standard deviations. We randomly selected ten grids in each geographic division to analyze their variability characteristics (Figure 3). Most of the grids show similar variation trends: at SM of 0–0.4 cm³/cm³, the standard deviations for SM increased and then decreased with increasing SM and peaked at about 0.2 cm³/cm³; when the SM was greater than 0.4 cm³/cm³, the variability in SM tended to increase significantly. This convex distribution and renewed growth trend were influenced by hydraulic parameters, and the variability was most sensitive to the pore size distribution parameter n [48], followed by f, θ_s, and α. In environments with homogeneous soil textures, sandy soils have a more uniform SM distribution due to larger pores and better water permeability, whereas clay soils are more prone to localized ponding or drying due to small water transport rates [77].

Figure 4 shows the SM variability distributions under specific SM (i.e., 0.1 cm³/cm³, 0.2 cm³/cm³, 0.3 cm³/cm³, and 0.4 cm³/cm³). Corresponding to the SM variability curves, the standard deviations for SM exhibited very similar spatial distribution characteristics at SM of 0.1 cm³/cm³ and 0.3 cm³/cm³; the standard deviations for SM were generally greater when the SM was 0.2 cm³/cm³. In terms of regional distribution, the Rocky Mountain subdivision has a diverse soil texture composition (Figure 1), and its variability under different SM conditions is higher (also confirmed in Figure 3). In contrast, the central Loess Hills region exhibited lower variability due to its single soil type composition, i.e., dominated by large, continuously distributed Loessial soil [78,79,80]. The above findings are similar to those of previous researchers [48,50,81], indicating that this SM variability estimation method had good applicability in the Loess Plateau and could be used as a basis for downscaling.

3.3. Spatial Distribution of Downscaled Soil Moisture

In most regions, the downscaled results obtained using both the traditional and improved algorithms were able to enrich the details of the original product while preserving its spatial distribution patterns. Since the downscaled SM was calculated from the ERA5-Land product, the grid edges of the original product were still more obvious in some areas, with large differences in SM in the coarse grids. Similar characteristics exist when other downscaling methods are used when the resolution difference between the original and downscaled images is large [31,82]. Therefore, the low-resolution SM and variability data were interpolated before downscaling to make the downscaled image smoother.

In order to show the effect of the improved algorithm more clearly, Figure 5 compares the details of the spatial distribution of the downscaled results and auxiliary data in different regions. It is evident that the downscaled results based on interpolation effectively alleviated the “edge effect” of the coarse grid and show richer detail features. In addition, the downscaled results have strong spatial heterogeneity in the Hilly Loess, the Rocky Mountain, and the Loess Yuan (Region A and B), while their spatial texture show high similarity with NDVI and DEM, reflecting the dominant roles of these two factors gained due to the importance analysis and the setting of the weighting factors. However, the downscaled results show a slight lack of detail in the northwest Windy Sand Hills (Region C), because the lower SM and more homogeneous soil texture in this area caused low SM variability and a small degree of spatial differentiation of the high-resolution auxiliary data, which was nevertheless consistent with the actual situation.

Figure 6 shows the 250 m resolution downscaled results obtained via the improved downscaling algorithm based on subdivision and interpolation. The spatial distribution of the downscaled results was relatively consistent with that of the original ERA5-Land SM, generally showing an increasing trend from northwest to southeast, similar to the spatial distribution pattern for precipitation, which was the consequence of the integrated effects of climate, topography, vegetation type, and soil type. The northwestern part of the Loess Plateau is a desert steppe zone with sparse vegetation and loose, porous soil, and therefore poor water retention; the southeastern part is affected by monsoons with frequent precipitation and a high vegetation cover, which can effectively attenuate the evaporation of soil water. On the timescale, because precipitation on the Loess Plateau is mostly concentrated in summer and SM is mainly replenished by precipitation, SM from July to September was significantly higher than that from April to June, except for the Northwest Arid Zone, where SM was generally low in different seasons.

3.4. Accuracy Evaluation

In this study, the raw product and downscaled data from April to September 2011–2013 were validated and evaluated by using station observation data.

Table 2. Statistics of the mean values of the site validation metrics.

Method	R	Bias (cm3/cm3)	RMSE (cm3/cm3)	ubRMSE (cm3/cm3)
Orig	0.718	0.045	0.082	0.044
TD	0.716	0.044	0.081	0.044
TDI	0.718	0.044	0.081	0.044
ID	0.718	0.040	0.078	0.044
IDI	0.721	0.040	0.077	0.044
IDI(GD)	0.753	0.040	0.071	0.042

Orig is the original ERA5-Land product, TD is the traditional downscaling algorithm, TDI is the traditional downscaling algorithm based on interpolation, ID is the improved downscaling algorithm, IDI is the improved downscaling algorithm based on interpolation, and IDI(GD) is the improved downscaling algorithm based on geographical division and interpolation.

shows the downscaled results with the statistics of the mean values of the evaluation indices at all of the sites, which indicates that different downscaling methods could improve the accuracy of the raw products. For different downscaling algorithms, the downscaled results based on interpolation were better than the results of direct downscaling, and the improved downscaling algorithm based on geographical division and interpolation had the highest accuracy. The improved downscaling methods improved the R values of the raw products to some extent, and, compared with the traditional downscaling algorithms, the improved downscaling algorithms based on interpolation further reduced the bias and RMSE values, but there was almost no improvement in the ubRMSE, and only the geographical division-based downscaling algorithm reduced the ubRMSE, which proved the effectiveness of the geographic partitioning strategy.

Figure 7 compares the distribution of R, Bias, RMSE, and ubRMSE for sites in each geographical division (since there is only one site in the Earth-rock Hills, it was combined into Rocky Mountain for error statistics). In terms of R and ubRMSE, Loess Yuan possessed the highest R and the lowest ubRMSE with less dispersion, followed by Hilly Loess. This may be related to the fact that vegetation and topography in these two regions scored high in importance and possessed stronger correlations with SM, making the downscaled results closer to the measured SM.

Figure 8 illustrates the spatial distribution of the validation results of the improved algorithm based on geographical division and interpolation. The correlation between downscaled and measured SM was overall high in the south and low in the north, and was broadly similar to the SM content spatial distribution. The correlation coefficients were greater than 0.7 at 92 sites and greater than 0.9 at 19 sites. Combined with the distribution of the mean bias, the downscaled results reveal underestimations of low-SM areas, most of which were distributed in the northwestern arid zone, while a few were located in the southeastern part of the Loess Plateau; additionally, the results reveal overestimations of high-SM areas, which were distributed mainly in the central and southwestern parts of the Loess Plateau. In addition, 89 sites had an RMSE less than 0.08 cm³/cm³, 106 sites had an ubRMSE less than 0.05 cm³/cm³, and, in general, the downscaled SM was able to satisfy the accuracy needs of the study area.

The gain distribution of the improved downscaling algorithm based on geographical division and interpolation (Figure 9) illustrated that most sites showed synergistic improvements in terms of the gains of RMSE and R. A small number of sites presented improved RMSE values, but decreased R values. Overall, the improvement in SM accuracy was greater in the southeastern and southwestern regions. Combined with the error distribution, the G_RMSE gain was distributed mainly between RMSE of 0.07–0.17 cm³/cm³, which was able to reduce the higher RMSE, but there was no significant correlation between the G_PREC and R distributions. The above validation results show that the improved downscaling algorithm based on geographical division and interpolation was superior to the traditional algorithm and could not only significantly increase the spatial resolution of ERA5-Land product, but also effectively improve its accuracy.

4. Discussion

4.1. Selection and Impact of Environmental Factors

Even though previous researchers estimated the spatial variability in SM and realized the downscaling in SM by using soil texture data and the MvG model, the effects of vegetation and topographic factors on SM were neglected. To make the downscaling process more scientific, we expanded the original single auxiliary factor, considering that the weight of each factor decreases with an increasing number of factors and the dilution effect of this weight distribution weakens the individual contribution of each factor. Therefore, we added only four typical environmental variables to the previous downscaling algorithm, which were the vegetation factor, NDVI, and the three terrain factors DEM, SLOPE, and TRASP [83,84,85].

In order to clarify how the selected factors affect SM, we plotted the partial dependence of each factor in geographic subdivisions (Figure 10). The higher the importance of an environmental factor, the greater the slope of its corresponding curve, which was consistent with the results of the correlation analysis. When NDVI was less than 0.3, an increase in NDVI significantly enhanced SM; however, its effect slowed down when NDVI exceeded 0.4. This phenomenon reflected the limiting effect of SM on vegetation growth and the nutrient-supporting effect of increased vegetation on SM [86]. Adequate SM supports healthy vegetation growth, while increased vegetation rooting and surface cover enhance precipitation retention and infiltration, so higher NDVI generally corresponds to higher SM.

The positive correlation between SM and DEM is mainly because the increase in altitude is accompanied by a decrease in temperature, which reduces the evapotranspiration potential and thus reduces SM loss [87]. Meanwhile, topographic uplift effect increases precipitation on windward slopes [88], providing more replenishment of SM at higher elevations. However, in the Flood Plain area, agricultural activities and anthropogenic disturbances are more frequent, resulting in a negative correlation between DEM and SM. The increase in SLOPE positively affected SM accumulation through pooling, whereas the increase in TRASP represented an intensification of the dry thermal environment, making SM more susceptible to evaporation. Notably, because the aggregation process greatly smoothed the auxiliary variables such as TRASP and SLOPE, which have high spatial heterogeneity, fine-scale topographic features were significantly masked, resulting in the TRASP failing the significance test with a 99% confidence interval in the overall analysis, as well as the extremely weak correlation between SLOPE and SM in some areas. However, related studies have shown that fine-scale topographic data are essential for generating high-resolution SM data [89]. Therefore, it is still necessary to retain the TRASP factor as an auxiliary variable.

4.2. Improvement and Applicability of the Model

This study employed RF to investigate the significance of environmental factors on SM, incorporating vegetation and topographic factors into the downscaling process. This approach enhanced the traditional downscaling method that relied on a single auxiliary variable, thereby improving the accuracy of the downscaled SM. Furthermore, this study integrated the interpretability advantages of physical models with the capability of machine learning models in exploring complex multi-factor relationships, providing a reference for the integration of downscaling methods. The improved algorithm based on geographical division and interpolation increased the positive G_RMSE of the traditional algorithm from 83 to 100 sites, with mean values from 0.008 to 0.022, and the positive G_PREC from 73 to 84, with a mean value of 0.002 to 0.004, indicating that the improved algorithm improved the ERA5-Land SM accuracy more significantly in both RMSE and R. Moreover, compared to the overall downscaled results, the geographical division-based downscaled results possessed lower RMSE and higher G_RMSE, which is the only result with reduced ubRMSE, indicating that setting up weights with regionality can further increase the accuracy of downscaled SM.

The ERA5-Land product, as the first term of Equations (15) and (16), is the most critical factor in determining the downscaled results. Therefore, the accuracy of the downscaled results depends heavily on the performance of the ERA5-Land SM, which has the same error distribution as the original product, i.e., underestimations of SM in dry areas and overestimations of SM in wet areas. In addition, SM variability, which serves as the basis for the decomposition of the raw products, determines the differences in the numerical distributions of the downscaled results. If the SM variability in a region is relatively low, the calculated results will be similar to the values of the original products, even if the high-resolution ancillary data have some variability. Especially in the desert and grassland zones, the accuracy of the downscaled results is highly similar to the original SM products due to the generally low SM and homogeneous soil texture [90], which results in low SM variability, and the auxiliary data, such as NDVI and DEM, are also highly homogeneous in this region. However, in the Rocky Mountain area and the southeastern humid region, with complex soil texture [91] and higher SM variability, and the spatial distribution of vegetation and topography varies significantly, downscaled results with more detailed spatial texture and higher accuracy can be obtained. Therefore, the method is more suitable for downscaling in areas with high spatial heterogeneity. Furthermore, in Equations (15) and (16), downscaling factors were used for coarse-scale SM decomposition, with underlying assumptions based on positive and negative correlations between the factors and SM. Therefore, the units of downscaling factors must have a clear physical meaning. For example, land cover data cannot be used as a downscaling factor, otherwise the downscaled SM will lose their scientificity and reliability.

4.3. Limitations and Prospects

Although this study improved the downscaling method and achieved better results, there are still some limitations. First, the high-resolution soil property datasets obtained from OpenLandMap were employed, which have some uncertainties in the number of samples used in the production process and in the degree of accuracy [92]. Having access to high-precision high-resolution soil property data would improve the computational accuracy of the MvG parameters and thus more accurately characterize the SM variability of the ERA5-Land product grid. During the validation process, although there is a scale mismatch between the validation data and the downscaled results, the observation stations are uniformly distributed throughout the Loess Plateau, and the correlation is improved while the error is reduced after downscaling, so the validation of the downscaled results can be considered to have good validity. In addition, the monthly scale station observations are synthesized by averaging, and the reliability of the validation results can be further improved if daily scale station observations are available, thus avoiding the secondary calculation of the downscaled results. In the future, surface-scale SM inversion using airborne sensors (e.g., hyperspectral, thermal infrared sensors) can be performed as ideal validation data, but it requires substantial human and material resources. In the practical application of this method, it is recommended that a systematic accuracy assessment and comparative analysis of multiple SM products be conducted first. Adopting SM products with higher precision can significantly improve the accuracy of downscaling SM. In future research, regional environmental characteristics can be fully considered to select regionally representative environmental variables as downscaling factors, which will help to obtain more targeted downscaling results.

5. Conclusions

In this study, the spatial variability of coarse-scale SM was estimated based on the MvG model, and the importance of multiple environmental factors within each subdivision was calculated using RF, and the downscaling algorithm was improved to obtain ERA5-Land SM downscaling data with 250 m spatial resolution in the Loess Plateau, which was validated using measured SM data. The improved downscaling algorithm adopted in this study can be a reference for subsequent studies on SM downscaling in arid and semi-arid areas, while the obtained downscaled SM can provide data support for agricultural development, vegetation construction, and ecological monitoring in this region. The following conclusions can be drawn:

(1) In the analysis of the importance of the selected environmental factors, NDVI and DEM ranked in the top two places overall in terms of the importance of influencing SM, although the results varied slightly in different geographical divisions. These two factors have the greatest influence on the spatial and temporal distribution of SM, which is superior to the FC selected using the traditional algorithm and is more applicable to the downscaling of SM.

(2) Using the PTF empirical formula to transform the soil texture into the hydraulic parameters of the MvG model, we successfully estimated the coarse-scale SM variability in the study area, which exhibited a convex distribution pattern when the SM was less than 0.4 cm³/cm³.

(3) By interpolating coarse-scale SM and variability, the edge effect due to the large difference in spatial resolution can be alleviated. A comparison of the accuracy of multiple downscaled results using site-observed SM revealed that the improved downscaling algorithm based on geographical division and interpolation had the highest accuracy and was significantly better than the traditional algorithm, confirming the effectiveness of the algorithm proposed in this study.