Estimation of subsurface soil moisture from surface soil moisture in cold mountainous areas

Profile soil moisture (SM) in mountainous areas is important for water resource management and ecohydrological studies of downstream arid watersheds. Satellite products are useful for providing spatially distributed SM information but only have limited penetration depth (e.g., top 5 cm). In contrast, in situ observations can provide measurements at several depths, but only with limited spatial coverage. Spatially continuous estimates of subsurface SM can be obtained from surface observations using multiple methods. This study evaluates methods to calculate subsurface SM from surface SM and its application to satellite SM products, based on a SM observation network in the Qilian Mountains (China) that has operated since 2013. Three different methods were tested to estimate subsurface SM at 10 to 20, 20 to 30, 30 to 50, and 50 to 70 cm, and, in a profile of 0 to 70 cm, from in situ surface SM (0 to 10 cm): the exponential filter (ExpF), the artificial neural network (ANN), and the cumulative distribution function (CDF) matching methods. The ANN method had the lowest estimation errors (RSR), while the ExpF method best captured the temporal variation of subsurface soil moisture; the CDF method is not recommended for the estimation. Meanwhile the ExpF method was able to provide accurate estimates of subsurface soil moisture at 10 to 20 cm and for the profile of 0 to 70 cm using surface (0 to 10 cm) soil moisture only. Furthermore, it was shown that the estimation of profile SM was not significantly worse when an area-generalized optimum characteristic time (T_opt) was used instead of station-specific T_opt for the Qilian Mountains. The ExpF method was applied to obtain profile SM from the SMAP_L3 surface soil moisture product, and the resulting profile SM was compared with in situ observations. The ExpF method was able to estimate profile SM from SMAP_L3 surface data with reasonable accuracy (median R of 0.65). Also, the combination of the ExpF method and SMAP_L3 surface product can significantly improve the estimation of profile SM in mountainous areas compared to the SMAP_L4 root zone product. The ExpF method is useful and has potential for estimating profile SM from SMAP surface products in the Qilian Mountains.

1 Introduction

Soil moisture (SM) is considered to be an essential climate variable (Bojinski et al., 2014) because of its critical role in the water, energy (Jung et al., 2010), and carbon cycles (Green et al., 2019). In particular, knowledge of profile SM is important for runoff modeling (Brocca et al., 2010), water resource management (Gao et al., 2018), drought assessment (Jakobi et al., 2018), and climate analysis (Seneviratne et al., 2010). Methods for SM measurements include ground-based measurements and satellite-based measurements (Dobriyal et al., 2012). Most ground-based methods enable the determination of SM changes with high temporal resolution at different depths but with limited spatial coverage (Jonard et al., 2018). Especially in mountainous regions, measuring SM in situ for a large area is difficult, and thus these measurements are scarce (Ochsner et al., 2013). In addition, strong SM heterogeneity in complex mountainous areas makes SM estimation over large areas more difficult (Williams et al., 2009). By comparison, satellite estimates of SM, such as those from the Soil Moisture Active & Passive (SMAP) mission, provide spatial SM coverage for large areas (Entekhabi et al., 2014; Brocca et al., 2017). Unfortunately, SMAP and other microwave-based SM products from spaceborne sensors only provide SM estimates for a limited depth up to ∼5 cm (Escorihuela et al., 2010). Thus, a gap exists with respect to the availability of subsurface SM information with adequate spatial coverage.

Previous studies have shown that subsurface SM is often related to surface and near-surface SM (Mahmood and Hubbard, 2007; Wang et al., 2017). A variety of methods for estimating subsurface SM from surface SM information have been developed, including data assimilation of remote sensing data into land surface models (Han et al., 2013), physically based methods (Manfreda et al., 2014), (semi-)empirical methods (Albergel et al., 2008), data-driven methods (Kornelsen and Coulibaly, 2014; Zhang et al., 2017a), and statistical methods (Gao et al., 2019). Among them, the application of both data assimilation and physically based methods are limited to data-rich areas due to the large amount of required input data, e.g., soil properties, which are often not available for data-scarce mountainous areas (Jin et al., 2015; Li et al., 2017; Dai et al., 2019). The cumulative distribution function (CDF) matching method is a statistical method developed to adjust systematic differences in different SM datasets (e.g., in situ observations and satellite products) based on observation operators (Drusch et al., 2005; Peng et al., 2017). CDF matching can also be used for upscaling of SM (Han et al., 2012) and estimating subsurface SM from surface SM (Gao et al., 2019). The artificial neural network (ANN) method is an effective and powerful data-driven tool for nonlinear estimation problems and has been widely used to estimate subsurface SM from surface SM measurements (Kornelsen and Coulibaly, 2014; Pan et al., 2017). The exponential filter (ExpF) method is a semi-empirical modeling approach and relies on a two-layer SM balance equation (Wagner et al., 1999). This method has been widely applied with both in situ observations and satellite products, and the performance of the ExpF method for estimating subsurface SM varied considerably over regions with different environmental conditions (Ford et al., 2014; González-Zamora et al., 2016; Tobin et al., 2017; Wang et al., 2017; Zhang et al., 2017a). Ford et al. (2014) found that root zone SM estimated from SMOS satellite products had a mean R² of 0.57 (ranging from 0.00 to 0.86) and 0.24 (ranging from 0.00 to 0.51) for SM networks in Oklahoma and Nebraska, respectively. In addition to surface SM data, the ExpF method requires only one additional parameter (T, the characteristic time) that reflects the combined influence of local conditions on the temporal characteristics of SM (Albergel et al., 2008; Ceballos et al., 2005). Previous studies have shown that T varied among different stations, and several methods have been developed to estimate T (Wagner et al., 1999; Albergel et al., 2008; Brocca et al., 2010; Qiu et al., 2014).

Methods for estimating subsurface SM from surface SM have not previously been evaluated for high and cold mountainous areas using in situ SM observations across a wide area. In the absence of in situ SM observation networks over a wide area, satellite SM products can be an alternative for providing surface SM information for a wide area (Ochsner et al., 2013). Although SM estimation from spaceborne sensors is especially challenging for mountainous regions, some validation studies have shown adequate accuracy (Pasolli et al., 2011; Rasmy et al., 2011; Zhao et al., 2014; Zeng et al., 2015; Zhao and Li, 2015; Colliander et al., 2017; Ullah et al., 2018; Qu et al., 2019; Liu et al., 2019). Nevertheless, the accuracy of profile SM estimation from remotely sensed SM products is currently unknown for mountainous regions.

In this study, we focus on the Qilian Mountains, which is a water source for several key inland rivers with terminal lakes in Northwest China, including the Heihe, Shiyang, and Shule rivers (He et al., 2018). Water scarcity threatens both food and ecosystem security in these endorheic basins (Feng et al., 2019). At the northeastern border of the Tibet–Qinghai Plateau, with its significant role in the Asian monsoon, profile water content in the Qilian Mountains is a key variable in ecohydrological studies on water resources and exchange processes in these basins (Zhao et al., 2013). Therefore, the aim of this study is to use in situ SM observations from 35 stations and remotely sensed SM data from the Qilian Mountains, a prime example of a high and cold mountainous area, to characterize the relationship between surface SM and subsurface SM in order to obtain the spatial distribution of profile SM. We first evaluated the performance of the different methods for estimating subsurface SM. We then employed the best method with SMAP surface SM products to evaluate the utility of this method for estimating profile SM in mountainous regions.

2 Study area

This study was carried out in the upland area of the Heihe River basin, which is a typical terminal lake basin of an arid region (Liu et al., 2018) (Fig. 1). It is located in the Qilian Mountains at the northeastern border of the Qinghai–Tibet Plateau. It covers approximately 2.7×10⁴ km², and the elevation ranges from about 2000 to 5000 m (Yao et al., 2017). The region has an annual precipitation ranging from 200 to 500 mm (Luo et al., 2016), annual potential evapotranspiration ranging from 700 to 2000 mm, and an annual mean temperature ranging from −3.1 to 3.6 ^∘C from 1960 to 2012 (He et al., 2018). The main land covers are grassland, forestland, and sparsely vegetated land (Zhou et al., 2016). The main soil types are Calcic Chernozems, Kastanozems, and Gelic Regosols. The main soil texture classes are silt loam, silt, and sandy loam (Tian et al., 2017, 2019).

Figure 1 (a) Study area and (b) distribution of the SM stations with spatial distribution of annual average precipitation from 2014 to 2016.

3 Data and methods

3.1 Datasets

We established a SM monitoring network in September 2013 in the Qilian Mountains. The network is composed of 35 SM stations distributed over the entire study area (Fig. 1). At each station, SM profiles from 0 to 70 cm were measured by soil moisture probes (ECH2O 5TE, METER Group Inc., USA) at 30 min intervals. These probes were installed at depths of 5 (representing a depth of 0 to 10 cm, SM_5 cm), 15 (10 to 20 cm, SM_15 cm), 25 (20 to 30 cm, SM_25 cm), 40 (30 to 50 cm, SM_40 cm), and 60 cm (50 to 70 cm, SM_60 cm) below the soil surface. Soil-specific sensor calibrations were performed with the direct calibration method using soil samples taken from each station (Cobos and Chambers, 2010; Zhang et al., 2017b). The profile-integrated SM (SM_0–70 cm) was calculated by the method of González-Zamora et al. (2016):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathrm{SM}}_{\text{0–70\hspace{0.17em}cm}}\\ & ={\displaystyle \frac{{\mathrm{SM}}_{\mathrm{5}\mathrm{cm}}\times \mathrm{10}+{\mathrm{SM}}_{\mathrm{15}\mathrm{cm}}\times \mathrm{10}+{\mathrm{SM}}_{\mathrm{25}\mathrm{cm}}\times \mathrm{10}+{\mathrm{SM}}_{\mathrm{40}\mathrm{cm}}\times \mathrm{20}+{\mathrm{SM}}_{\mathrm{60}\mathrm{cm}}\times \mathrm{20}}{\mathrm{70}}}.\end{array}\end{array}$$

The entire dataset used in this study thus consists of six in situ SM time series at depths of 5, 15, 25, 40, 60, and 0 to 70 cm for each of the 35 stations. Due to the influence of soil freezing in winter, the soil moisture time series was limited to the growing seasons (May to October, Tian et al., 2019) of 2014, 2015, and 2016. The measurements were averaged to obtain daily SM, following the approach of Wagner et al. (1999). Data quality management was performed for each station, and data gaps existed in the harsh mountainous environment, as described in detail in Tian et al. (2019). Time series where more than 50 % of observations were missing were excluded from further analysis. The final dataset after processing is presented in Fig. 2. The surface SM measured at 5 cm was used to predict the subsurface SM at depths of 15, 25, 40, and 60 cm and the profile average (0 to 70 cm).

Soil cores were taken to measure soil properties, including soil organic carbon (SOC), saturated hydraulic conductivity (K_S), soil particle composition, and bulk density for each layer during the sensor installation. Detailed descriptions of the soil properties can be found in Tian et al. (2017, 2019). The statistics of the soil physical characteristics are provided in Table 1. Daily reanalysis precipitation product (Chen et al., 2011) and Landsat-based continuous monthly 30 m ×30 m resolution NDVI data for the period 1986 to 2017 (Cihlar et al., 1994; Huete et al., 2002; Wu et al., 2019) were acquired from the National Tibetan Plateau Data Centre (https://data.tpdc.ac.cn/en/, last access: 18 September 2020).

The widely used higher-level SMAP_L3 Global Daily 9 km product for the growing seasons of 2015 to 2017 was used in this study. This product is distributed by NASA (http://nsidc.org/, last access: 18 September 2020) and described by O'Neill et al. (2018). SMAP descending node observations acquired near 06:00 local solar time have been combined with global daily composites in order to reduce the impact of Faraday rotation and to consider the assumption of uniform temperature profiles in the vegetation cover during morning overpasses. It has to be noted that the data are provided on a 9 km grid but that this is a result of a Backus–Gilbert optimal interpolation at brightness temperature level. The actual spatial resolution is coarser (O'Neill et al., 2018). The SMAP_L3 surface soil moisture product was also used to estimate the subsurface soil moisture (layer 2: 10 to 20 cm, layer 3: 20 to 30 cm, layer 4: 30 to 50 cm, layer 5: 50 to 70 cm) and profile soil moisture (0 to 70 cm) during the growing seasons of 2015 and 2016 in the mountainous area.

SMAP_L4 provides estimates of both surface and root zone SM products based on the assimilation of brightness temperature into the NASA land-surface model and has a spatial and temporal resolution of 9 km and 3 h, respectively (Reichle et al., 2017). SMAP_L4 is a widely used root zone SM product (Pablos et al., 2018). Here, the SMAP_L4 data were averaged to a daily resolution in order to compare it with the profile SM estimates from the SMAP_L3 surface product obtained in this study. In particular, the SMAP_L4 SM products with both surface (0 to 5 cm, sm_0–5) and root zone (0 to 100 cm, sm_0–100) information were used to calculate SM of the 0 to 70 cm profile (sm_0–70) using

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathrm{sm}}_{\text{0–100}}=(\mathrm{5}\cdot {\mathrm{sm}}_{\text{0–5}}+\mathrm{95}\cdot {\mathrm{sm}}_{\text{5–100}})/\mathrm{100},\text{(3)}& {\displaystyle }{\mathrm{sm}}_{\text{0–70}}=(\mathrm{5}\cdot {\mathrm{sm}}_{\text{0–5}}+\mathrm{65}\cdot {\mathrm{sm}}_{\text{5–100}})/\mathrm{70}.\end{array}$$

Table 1 Statistics of the physical characteristics of the soil at the 35 soil moisture stations: mean (standard deviation).

Note: K_S is the saturated hydraulic conductivity; SOC is the soil organic carbon.

Figure 2 Daily soil moisture (vol. %) time series during the growing season of 2014 to 2016 for the five layers (layer 1, 0 to 10 cm; layer 2, 10 to 20 cm; layer 3, 20 to 30 cm; layer 4, 30 to 50 cm; layer 5, 50 to 70 cm) in the 35 soil moisture stations. Gaps exist for some stations due to missing data.

3.2 Exponential filter (ExpF) method

The ExpF method predicts the dynamics of subsurface SM using an exponential filter function of the surface SM dynamics (Wagner et al., 1999; Albergel et al., 2008). First, SM (cm³ cm⁻³) is transformed into a soil water index (SWI) with

$$\begin{array}{}\text{(4)}& {\mathrm{SWI}}_i={\displaystyle \frac{{\mathit{\theta }}_i-{\mathit{\theta }}_{i,\min }}{{\mathit{\theta }}_{i,\max }-{\mathit{\theta }}_{i,\min }}},\end{array}$$

where θ_i,min and θ_i,max are the minimum and maximum SM in the time series collected since installation for each layer at each station (Ford et al., 2014). The ExpF method then estimates subsurface SM from surface SM using

$$\begin{array}{}\text{(5)}& {\mathrm{SWI}}_{m,t_n}={\mathrm{SWI}}_{m,t_{n-\mathrm{1}}}+K_{t_n}(m{s}_{t_n}-{\mathrm{SWI}}_{m,t_{n-\mathrm{1}}}),\end{array}$$

where SWI${}_{m,t_{n-\mathrm{1}}}$ and SWI${}_{m,t_n}$ are the predicted subsurface SWI at time t_n−1 and t_n, respectively. ms${}_{t_n}$ is the observed surface SWI at time t_n, and $K_{t_n}$ represents the gain at time t_n calculated by

$$\begin{array}{}\text{(6)}& K_{t_n}={\displaystyle \frac{K_{t_{n-\mathrm{1}}}}{K_{t_{n-\mathrm{1}}}+e^{-\frac{t_n-t_{n-\mathrm{1}}}{T}}}},\end{array}$$

where $K_{t_{n-\mathrm{1}}}$ is the gain at time t_n−1 and T is the characteristic time in days. The equation was initialized with SWI${}_{m,t_{\mathrm{1}}}={\mathrm{ms}}_{t_{\mathrm{1}}}$ and $K_{t_{\mathrm{1}}}=\mathrm{1}$ (Albergel et al., 2008). This method is particularly useful as T is the only unknown parameter. The optimum T (T_opt) was determined by optimization using the highest Nash–Sutcliffe score for each specific depth at each station.

3.3 Artificial neural network (ANN) method

The ANN method is a data-driven method to predict subsurface SM from surface SM (Zhang et al., 2017a). If properly trained, an ANN can describe nonlinear relationships between dynamics of SM at different depths (Kornelsen and Coulibaly, 2014). The commonly used feed-forward ANN (with one hidden layer and 10 neurons, Levenberg–Marquardt algorithm, Ford et al., 2014) was used in this study. The ANN modeling was carried out using MATLAB (neural network time series tool, R2017b, The MathWorks). The output of the ANN was calculated using

$$\begin{array}{}\text{(7)}& y=f\left[W_{\mathrm{2}}g\left(W_{\mathrm{1}}X+b_{\mathrm{1}}\right)+b_{\mathrm{2}}\right],\end{array}$$

where y is the output (the estimated subsurface soil moisture), f and g are the activation functions of the hidden layer and the input layer (the surface soil moisture), respectively, W₁ and W₂ are the weights of the input layer and the hidden layer, respectively, and b₁ and b₂ are the biases of the input layer and the hidden layer, respectively. The tangent sigmoid function was used as the activation function as it has shown good performance in hydrological studies (Yonaba et al., 2010). As suggested by Zhang et al. (2017a), 70 % of the data were selected for training the ANN and the remaining 30 % were used for validation. A separate ANN model was developed for every depth combination and every site.

3.4 Cumulative distribution function (CDF) matching method

In this study, the following procedure for CDF matching was used.

(1) Rank the surface (θ₁) and the subsurface SM (θ₂) time series.

(2) Calculate the difference between the two observation time series:

$$\begin{array}{}\text{(8)}& {\mathrm{\Delta }}_i={\mathit{\theta }}_{\mathrm{1},i}-{\mathit{\theta }}_{\mathrm{2},i}.\end{array}$$

(3) Use a cubic polynomial fit to relate the difference (Δ) to surface SM (θ₁) as recommended by Gao et al. (2019):

$$\begin{array}{}\text{(9)}& \widehat{\mathrm{\Delta }}=K_{\mathrm{0}}+K_{\mathrm{1}}\cdot {\mathit{\theta }}_{\mathrm{1}}+K_{\mathrm{2}}\cdot {\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}+K_{\mathrm{3}}\cdot {\mathit{\theta }}_{\mathrm{1}}^{\mathrm{3}},\end{array}$$

where $\widehat{\mathrm{\Delta }}$ is the predicted difference between surface and subsurface SM, and K_i ($i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}$) are parameters.

(4) Calculate CDF-matched subsurface SM (θ_CDF) with

$$\begin{array}{}\text{(10)}& {\mathit{\theta }}_{\mathrm{CDF}}={\mathit{\theta }}_{\mathrm{1}}-\widehat{\mathrm{\Delta }}.\end{array}$$

Similarly to the ANN method, 70 % of the data were used to calibrate the approach, and the remaining 30 % of the data were used for validation of the CDF matching method.

3.5 Statistical analysis

Boxplots were used to show the scatter of the data. The difference between data in different groups was examined using a one-way analysis of variance (ANOVA) with the post hoc Bonferroni test when the normality and homogeneity of variance of the datasets were satisfied. The Kruskal–Wallis ANOVA with a post hoc Dunn's test was used in cases where these conditions were not satisfied (Lange et al., 2008). The statistical analysis was performed in SPSS (SPSS 18.0, SPSS Inc.) and Matlab (R2017b, The MathWorks). The significance level was 0.05 for all statistical tests.

4 Results and discussion

4.1 Comparison of different methods

The ExpF method estimates subsurface SM based on SWI, while the ANN and CDF methods are based on volumetric soil moisture. Following Moriasi et al. (2007), the Nash–Sutcliffe efficiency (NSE), the ratio of RMSE to the standard deviation of the observations (RSR, an error statistic that normalizes the RMSE), and the Pearson correlation coefficient (R) were used to evaluate the performance of different methods with different units. To ensure that the comparison between the three methods is made under the same conditions, we divide the datasets into training data (the first 70 % of the data) and validation data (the remaining 30 % of the data) for all three methods. Figure 3 and Table 2 summarize the metrics (NSE, RSR, and R) for the subsurface SM estimates at different depths derived by the three different methods for the growing seasons of 2014, 2015, and 2016. Results show that ANN performed better than ExpF for the individual layers (layers 1 to 5) in terms of both NSE and RSR (Table 2 and Fig. 3), while ExpF performed better than ANN in estimating soil moisture for the entire soil profile. Additionally, the comparison of the performances between the ExpF and ANN methods was non-significant (p > 0.05) for all the layers, but ExpF showed a significantly (p < 0.05) higher R value compared to ANN for all the layers (with a median value of 0.97, 0.93, 0.84, 0.74, and 0.96 for layers 2, 3, 4, and 5, and profile SM, respectively). The good performance for R suggests that the ExpF method had the best ability to describe the temporal variability in SM. Furthermore, Table 2 and Fig. 3 indicate that CDF provided the worst performance among the three methods and thus cannot be recommended.

As expected, all metrics showed that the performance of the three methods decreased with depth. The results indicate that for two out of the three statistical measures (i.e., RSR and NSE), the ANN method was statistically superior to the other two methods. Specifically, the ANN method resulted in the lowest estimation error, while the ExpF method was better able to capture the SM dynamics. A similar finding was reported by Zhang et al. (2017a), who found that the ExpF method had a significantly higher correlation coefficient along with a higher mean bias compared to the ANN method. Furthermore, the ExpF method is a simpler approach as it only needs one parameter (T_opt) and can thus be easily applied in data-scarce mountainous areas, while the establishment of the ANN method is much more complicated. In addition, the ExpF method is a process-based method, while ANN is a machine learning method. Therefore, the ExpF method was used to estimate the subsurface SM in the remainder of this study.

Figure 3 Boxplot of the metrics (NSE, RSR, R) to compare the subsurface SM estimation using the surface SM by the three methods (ExpF, ANN, CDF) with the observations of the 35 stations during the growing seasons of 2014 to 2016. Different letters above the box indicate the significant difference (p < 0.05) among the different methods.

Table 2 The median of the performance (RSR, R, and NSE) of the three different methods (ExpF, ANN, and CDF) for estimating the subsurface SM using the surface SM for each layer of 35 stations during the growing seasons of 2014, 2015, and 2016.

4.2 Evaluation of T_opt for the ExpF method

4.2.1 Variation of T_opt with depth

In the method comparison, the first 70 % and the remaining 30 % of data were selected as training and validation data, respectively, to ensure the comparison was under the same condition. However, for the standard procedure of the ExpF method in earlier studies, the entire dataset is always used to derive the Topt and validate the ExpF method (e.g., Wagner et al., 1999; Albergel et al., 2008; De Lange et al., 2008; Ford et al., 2014; Wang et al., 2017). Thus, the ExpF method is evaluated and analyzed using the entire dataset as well (performance of the ExpF method using the entire dataset was shown in Table 3 and Fig. S1 in the Supplement). Results indicate that the performances of ExpF in both layer 2 and profile are significantly higher than that of other layers. Moreover, results also indicate that the ExpF method showed good performance for layer 2 and profile SM (with median NSE > 0.65, median RSR < 0.60, Moriasi et al., 2007).

Table 3 The statistics (mean±standard deviation and median) of the performance (RSR, R, and NSE) of the ExpF method for estimating the subsurface SM using the surface SM for each layer of 35 observation stations during the growing seasons of 2014, 2015, and 2016.

The accuracy of the ExpF method varied with the selected T value, and higher T values resulted in more stable estimations of SM time series (Wagner et al., 1999; Albergel et al., 2008). Furthermore, it was found that each station had an optimum T (T_opt) as determined by the best match with observations in terms of NSE. The variation of NSE with T (ranging from 0 to 68 d) for different layers for each station is shown in Fig. 4 and Table 4. The sensitivity of high values of NSE to changes in T decreased with increasing depth, indicating that the range of T values with high NSE was larger deeper in the soil. This was also observed in previous studies (e.g., Wang et al., 2017).

Figure 4 Variation of NSE with T of the exponential filter method for different layers at each station during the growing seasons of 2014, 2015, and 2016. The vertical axis is the NSE value. The frequency distribution curve and histogram show the distribution of T_opt with depth for all stations.

Results of a two-way ANOVA showed that the difference of T_opt is not significant between different years (p=0.06), while differences were significant between layers (p < 0.001). Furthermore, T_opt increased with depth from layer 2 to layer 5. The median of T_opt ranged from 1.5 d for layer 2 to 12.5 d for layer 5. The median T_opt for profile SM was 3.5 d. Significant differences in T_opt were obtained for layer 2, layer 3, and layer 4, but the difference between layers 4 and 5 was not significant. The increase in T_opt with depth has already been observed in many studies and is related to the greater temporal stability of SM in deeper soil layers (Wang et al., 2017; Tian et al., 2019).

Table 4 The statistics of T_opt (day) for each layer and different year for all stations.

Note: SD represents the standard deviation. This summary represents the statistical result of the 3 years. Letters in the summary row indicate significant difference between respective layers: the same letter in each column indicates that the difference is nonsignificant, while different letters indicate a significant difference between the two layers (p < 0.05).

4.2.2 Evaluation of alternative methods for T_opt estimation

Previous studies have used various methods to estimate T_opt. For example, Albergel et al. (2008) and Ford et al. (2014) found that using a single representative value for T_opt (e.g., average or median) for all stations did not significantly reduce the accuracy of the SM estimates. Wagner et al. (1999) recommended a common value of T_opt=20 (d) to estimate root zone SM, and this value has been widely adopted (e.g., Lange et al., 2008; Muhammad et al., 2017). Qiu et al. (2014) proposed to estimate T_opt using the station-specific long-term mean NDVI using $T_{\mathrm{opt}}=-\mathrm{75.263}\times \mathrm{NDVI}+\mathrm{68.171}$ (R=0.5, p < 0.01). This approach has also been applied in another study (Tobin et al., 2017).

Here, we evaluated four different methods to estimate T_opt in our study region for estimating profile SM (0 to 70 cm, SWI) from surface SM (5 cm, SWI). In the first method, T_opt was estimated from the NDVI-based regression of Qiu et al. (2014) to provide T_Qiu. In the second method, T_opt was set to 20 d as recommended by Wagner et al. (1999) to provide T_Wagner. In the third method, an area-generalized T_opt was obtained from the median value for the profile SM in our study region (3.5 d) to provide T_general. In the fourth method, the original station-specific T_opt parameter for profile SM was used (T_specific). The accuracy of the SM estimates obtained using the different methods to estimate T_opt was again evaluated using NSE, R, and RMSE (Fig. 5). The performance metrics show that T_specific performed best (mean RSR of 0.58, R of 0.88, and NSE of 0.61), followed by T_general (mean RSR of 0.61, R of 0.85, and NSE of 0.58), T_Wagner (mean RSR of 0.79, R of 0.69, and NSE of 0.32), and T_Qiu (mean RSR of 0.89, R of 0.59, and NSE of 0.17). However, the difference in performance between T_specific and T_general is not significantly different. The T_Wagner and T_Qiu approaches performed worse, and the metrics (NSE, R, RSR) are significantly (p < 0.001) lower than those of the T_general and T_specific methods. Our results suggest that a site-specific T_opt significantly improves the performance of the ExpF method compared to the use of the universal T_opt recommended by Wagner et al. (1999) or the regression of Qiu et al. (2014). Similarly, Lange et al. (2008) also found a significant improvement when using a station-specific T_opt instead of T_opt=20 d. It should be mentioned that the estimation depth in the method of Wagner et al. (1999) was 0 to 100 cm, while that of our study was 0 to 70 cm. This may partly explain the poor performance of the T_Wagner approach in this study. The use of an area-generalized T_opt (3.5 d) is a suitable alternative to T_opt estimation in our study area and provides similar estimation performance. Other studies have also found a good performance when using an area-generalized T_opt (e.g., Albergel et al., 2008; Brocca et al., 2010; Ford et al., 2014).

Figure 5 The boxplot of NSE, Pearson's R, and RSR for the T_opt generated from different schemes. The different letters above each box indicate the significant difference for different schemes.

4.3 Estimating profile soil moisture using SMAP

The ExpF method is suitable for estimating the profile SM from the surface SM and the median of T_opt is suitable for estimating subsurface soil moisture. Thus, in this section, we evaluate the utility of the ExpF method (with the median of T_opt from SMAP) in combination with SMAP surface products for estimating subsurface SM in mountainous areas.

4.3.1 Assessment of the SMAP surface SM product

The observed surface SM of each station was compared with the SMAP_L3 soil moisture product that overlapped with the corresponding station for the growing seasons of 2015 and 2016 for all stations to evaluate the accuracy of the SMAP measurements (Pablos et al., 2018). The root mean square error (RMSE), mean bias error (MBE), unbiased RMSE (ubRMSE), and R were adopted as metrics to evaluate accuracy. The relationship between the SMAP_L3 SM data product and the in situ observations at 5 cm depth are presented in Fig. 6. Clearly, the larger deviation from linearity in the relationship is due to the scale discrepancy between the relatively large satellite footprints and the point location of in situ SM measurements. Nevertheless, the statistical metrics still indicate a significant relationship between the SMAP_L3 SM data product and the in situ observations at 5 cm depth. The time series of the two datasets for each station are provided in Fig. S2. Figures 6 and S2 show that the performance was low at two stations (D13 with R of 0.18, D15 with R of 0.08) with scrubland and relatively high soil moisture. The poor performance at scrubland sites is consistent with results presented by Zhang et al. (2017b) for this study region. Results showed that the MBE varied from −0.23 to 0.07 cm³ cm⁻³ with a median of −0.021 cm³ cm⁻³. This indicates that SMAP underestimated surface SM over the study region, which is consistent with previous studies in the area (Chen et al., 2017; Zhang et al., 2017b). The RMSE varied between 0.026 and 0.250 cm³ cm⁻³ between sites with a median value of 0.052 cm³ cm⁻³. After removing the bias, the SMAP product had a median ubRMSE of 0.036 cm³ cm⁻³ (range from 0.024 to 0.083 cm³/cm³). Therefore, the SMAP product achieved the accuracy requirement of 0.04 cm³ cm⁻³ (Chan et al., 2016) in this study area. The R value ranged from 0.075 to 0.81 with a median value of 0.59. The relationship between SMAP-derived and in situ observed surface SM was significant (p < 0.05) at all but one station. This suggests that the SMAP surface product can represent the temporal dynamics of the observed surface SM time series.

Figure 6 The SMAP_L3 surface SM (cm³ cm⁻³) versus in situ observations at the surface (5 cm) for the 35 soil moisture stations. Color indicates station. The averaged metrics (RMSE, MBE, R, ubRMSE) are for all 35 stations during the growing seasons of 2015 and 2016.

Figure 7 Comparisons of SMAP_L3 estimated–observed subsurface SWI for all stations during the growing seasons of 2015 to 2016. The smoothed color density in the scatter plot shows the density of points more clearly. The dash and solid lines are the best-fitted curve and “y=x” line, respectively.

Table 5 Performance metrics (RSR, R, NSE) for the comparison of SMAP estimated and observed SWI at different layers for the 35 stations during the growing seasons of 2015 to 2016.

Note: the different letters after the layers indicate that the difference is significant at p < 0.05 (Kruskal–Wallis ANOVA).

Table 6 Statistics of the metrics (RSR, R, NSE) of the comparisons of estimated–observed profile SWI datasets, SMAP_L3-observed surface SWI datasets, SMAP_L3-observed profile SWI datasets, and SMAP_L4-observed profile SWI datasets for the 35 stations during the growing seasons of 2015 and 2016.

Note: e.g., estimated–observed PSWI means the comparison of the estimated profile SWI and observed profile SWI.

Figure 8 The spatial distribution of the monthly averaged profile SWI product estimated from the SMAP_L3 surface product during the growing seasons from 2015 to 2017. The title of each subplot provides the month and year.

4.3.2 SMAP-based estimation of subsurface soil moisture

For the estimation of subsurface soil moisture from the SMAP_L3 surface product, the site-specific T_opt was calculated based on the best match between SMAP estimations and in situ observations in terms of NSE. The median values of T_opt for the layers 2, 3, 4, and 5 and the profile are 7, 12, 22, 35, and 10 d, respectively. The subsurface SWI estimated from the combination of SMAP surface SM with the ExpF method (with the median values of T_opt) were compared with the in situ observations. A comparison of the subsurface SWI time series for different layers at each station are provided in Figs. S3 to S7. Figure 7 shows the measured SWI plotted against the predicted SWI. The performance metrics of these comparisons for each layer are summarized in Table 5.

As expected, the estimation accuracy of subsurface SM decreased with depth. The ANOVA results showed that the subsurface SM estimation accuracy for layer 2 (median value of RSR =0.92, R=0.69, NSE =0.18) and profile SM (RSR =0.92, R=0.65, NSE =0.14) was significantly higher than for layer 4 (RSR =1.12, R=0.31, NSE $=-\mathrm{0.13}$) and layer 5 (RSR =1.17, R=0.34, NSE $=-\mathrm{0.15}$) (p < 0.05). The NSE values were positive for layer 2 and profile SM, while the NSE values for the other layers were negative. The negative MBE shows that subsurface SM was underestimated. The relationship between SMAP-derived and in situ observed subsurface SM for layer 2 and profile SM was significant (p < 0.01) at all but one station (D15). Thus, the SMAP surface product and ExpF method can be used to estimate the subsurface SM in the study area, especially for layer 2 (10 to 20 cm) and profile (0 to 70 cm) SM.

As suggested by Ford et al. (2014), we partitioned the error in the SMAP-based estimation of profile SWI (“SMAP-observed profile SWI”, Fig. S8c) into errors associated with the ExpF method and errors due to SMAP observation differences to gain some insight into the error sources of SMAP-based estimates of profile SWI. For this, profile SWI estimated using the ExpF method from observed surface SWI was compared with in situ observed profile SWI (“estimated–observed profile SWI”) to assess errors of the ExpF method (Fig. S8a). In addition, SMAP-based and in situ observed surface SWIs (“SMAP-observed surface SWI”) were compared to assess inherent errors of the SMAP product (Fig. S8b). RMSE, R, and MAE were used as the metrics to assess accuracy. The results of this analysis are summarized in Table 6.

Figure S8 and Table 5 show that the SMAP-observed SWI had lower performance metrics for surface SWI (median values of RSR, R, and NSE are 1.01, 0.59, and −0.07, respectively) than for profile SWI (median values of RSR, R, and NSE are 0.88, 0.72, and 0.19, respectively), which was similar to the results obtained from the Nebraska SM network (Ford et al., 2014). This may be because the profile SWI was estimated based on the SMAP surface SWI and T_opt, which was determined by optimization using the maximum NSE. This may have improved the performance of profile SWI estimation. In addition, the performance metrics for SMAP–observed SWI comparisons for both surface and profile SWI were significantly (p < 0.001) lower than those of estimated–observed profile SWI (median values of RSR, R, and NSE are 0.68, 0.90, and 0.64, respectively). Thus, the major error in SMAP-based profile SWI estimates stems from the SMAP satellite product and is not derived from the ExpF method, which is also supported by previous studies (e.g., Ford et al., 2014; Pablos et al., 2018). As mentioned before, the scale mismatch between point measurements and satellite footprints will introduce additional errors in the validation of the satellite-derived subsurface products (Jin et al., 2017).

Subsequently, the SMAP_L4 and SMAP_L3 estimated profile SWIs were compared to the in situ observed profile SWI (see Fig. S9 and Table 5). Table 5 shows that the performance of profile soil moisture estimation using the SMAP_L3 surface product and the ExpF method (median RSR, R, and NSE of 0.92, 0.65, and 0.14, respectively) was significantly (p < 0.01) better than that of the SMAP_L4 product (median RSR, R, and NSE of 1.25, 0.55, and −0.3, respectively). The low performance of the SMAP_L4 profile product may be associated with uncertainty in the meteorological driving forces and the soil parameters in the NASA catchment model for cold mountainous areas (Reichle et al., 2017; Zhao et al., 2018; Dai et al., 2019). Thus, our results suggest that combining the exponential filter method with the SMAP_L3 product significantly improves the estimation of profile SM for the data-scarce cold arid mountainous areas.

Finally, the spatial distribution of profile soil moisture during the growing seasons of 2015, 2016, and 2017 was obtained using the median value of T_opt and the SMAP_L3 product to get the spatial distribution of profile SM in the study area (Fig. 8). Profile SM is higher in the southeast and lower in the northwestern part of the study area. This distribution coincides with the spatial distribution of precipitation and surface SM. The temporal variations of profile SWI, surface SWI, and precipitation are shown in Fig. S10. Figure S10 shows that the temporal variation of the SM profile corresponded well to the occurrence of precipitation: profile SM increased from May (mean SM of 0.27) to September (0.533) and then decreased until October (0.304). Profile SWI_SMAP (define SWI_SMAP before using it) was lower than surface SWI_SMAP from May to August, while profile SWI_SMAP was higher than surface SWI_SMAP from September to October. This can be attributed to the higher sensitivity of surface SM dynamics to precipitation and evapotranspiration (ET). During the months of September and October, less precipitation and higher ET caused a faster decrease in surface SM compared to profile SM.

Previous studies have shown the difficulty of applying the ExpF method to satellite products in mountainous areas, where complex topography (Paulik et al., 2014), snow, and soil freezing (Ford et al., 2014; Pablos et al., 2018) cause large errors and poor performance of the filtering method (Albergel et al., 2008). Ford et al. (2014) found an improvement of performance after removing the effects of snow from the data in the SCAN network, USA. In contrast, the present study showed that the ExpF method is useful in estimating profile SM from SMAP surface products in the growing season in high and cold mountainous areas, based on in situ SM observations.

5 Conclusions

We used three methods (the exponential filter (ExpF), the artificial neural network (ANN), and the cumulative distribution function matching (CDF) methods) to calculate subsurface SM from in situ surface SM observations at 5 cm depth in the Qilian Mountains (China). We also evaluated the utility of the ExpF method to estimate profile SM from SMAP surface products in the study area. Our main findings are the following.

1. With increasing depth of the predicted soil layer, the accuracies of all three methods decreased. The ExpF methods showed good performance for the estimation of SM down to 20 cm and profile.

2. The ANN method exhibited the lowest estimation error, while the ExpF approach captures the temporal variation of subsurface SM better than other methods.

3. The area-generalized T_opt value of the ExpF method can be used in the study area to estimate the subsurface SM without significantly reducing the performance compared to a station-specific T_opt.

4. Subsurface SM derived from the SMAP_L3 surface SM product using the ExpF method showed less deviation from the in situ observations compared to the SMAP_L4 root zone product for the study area.

We anticipate that our findings can improve the estimation of subsurface SM for large regions in mountainous areas, which in turn will support ecohydrological research and water resource management in inland river basins.