Data worth analysis within a model-free data assimilation framework for soil moisture flow

Conventional data worth (DW) analysis for soil water problems depends on physical dynamic models. The widespread occurrence of model structural errors and the strong nonlinearity of soil water flow may lead to biased or wrong worth assessment. By introducing the nonparametric data worth analysis (NP-DWA) framework coupled with the ensemble Kalman filter (EnKF), this real-world case study attempts to assess the worth of potential soil moisture observations regarding the reconstruction of fully data-driven soil water flow models prior to data gathering. The DW of real-time soil moisture observations after Gaussian process training and Kalman update was quantified with three representative information metrics, including the trace, Shannon entropy difference and relative entropy. The sequential NP-DWA framework was examined by a number of cases in terms of the variable of interest, spatial location, observation error, and prior data content. Our results indicated that, similarly to the traditional DW analysis based on physical models, the overall increasing trend of the DW from the sequential augmentation of additional observations within the NP-DWA framework was also susceptible to interruptions by localized surges due to never-experienced atmospheric conditions (i.e., rainfall events). The difference is that this biased DW in the former is caused by model structural errors triggered by contrasting scenarios, which is difficult to be compensated for by assimilating more prior data, while this performance degradation in the NP-DWA can be effectively alleviated by enriching training scenarios or the appropriate amplification of observational noise under extreme meteorological conditions. Nevertheless, a substantial expansion of the prior data content may cause an unexpected increase in the DW of future potential observations due to the possible introduction of ensuing observation noises. Hence, high-quality and representative small data may be a better choice than unfiltered big data. Compared with the observations in the surface layer with the strongest time variability, the soil water content in the middle layer robustly exhibited remarkable superiority in the construction of model-free soil moisture models. We also demonstrated that the DW assessment performance was jointly determined by 3C, i.e., the capacity of potential observation realizations to capture actual observations, the correlation of potential observations with the variables of interest and the choice of DW indicators. Direct mapping from regular meteorological data to soil water content within the NP-DWA mitigated the adverse effects of nonlinearity-related interference, which thus facilitated the identification of the soil moisture covariance matrix, especially the cross-covariance.

1 Introduction

As one of the few directly observable hydrological variables, soil water content (SWC) exhibits critical importance in optimal water resource management, irrigation and drainage schemes, fertilizer application, and crop production in agriculture (Liu et al., 2011; Akhtar et al., 2019; Dobriyal et al., 2012; Gu et al., 2021). Various data assimilation (DA) approaches (Dunne and Entekhabi, 2005; Li and Ren, 2011; Reichle et al., 2008; Song et al., 2014) have been established to reconstruct the spatiotemporal dynamics of SWC from noisy or partial observations. The core of these traditional parametric filters is their reliance on repeated forward integrations of an explicitly known physical model of unsaturated flow, such as the HYDRUS (Šimùnek et al., 2006), Soil and Water Assessment Tool (SWAT) (van Dam and Feddes, 2000), and Ross models (Ross, 2003; Zha et al., 2013).

Currently, the ever-increasing availability of multi-source data from remote sensing (Montzka et al., 2011; Shi et al., 2011), ground-based measurements (Li et al., 2018; Shuwen et al., 2005; Yang et al., 2000) and numerical modeling has paved the way for the development of fully data-driven techniques within the DA framework. In particular, recent advances in machine-learning-based DA schemes (Brajard et al., 2020, 2021; Yamanaka et al., 2019) offer exciting new opportunities for extracting patterns and insights of soil moisture dynamics from data (Ju et al., 2018; Li et al., 2020; Liu et al., 2020; Wang et al., 2021a). For instance, Kashif Gill et al. (2007) proposed a hybrid DA methodology that combined support vector machines and the ensemble Kalman filter (EnKF) for soil moisture dynamics. Li et al. (2020) compared the performance of a physical-based model with DA and machine learning methods in terms of the simulation of soil water dynamics under synthetic and real-world conditions. Wang et al. (2021a, b) further attempted to learn unknown relationships between SWC, as well as its spatio-temporal gradients and highly accessible data, via the Gaussian process (GP) regression.

Notwithstanding the success of these model-free DA schemes built on machine learning for unsaturated flow, essential caveats and limitations have hampered their further adoption and impact. First, the amount of data required to infer nonlinear relationships in unsaturated flow problems may be overwhelming (Hughes, 1968), thus greatly increasing the data collection budget. Subsequently, addressing the abovementioned explosive data growth is also a challenging and time-demanding task requiring an extensive computational infrastructure. Second, the performance and quality of the knowledge extracted by machine learning algorithms are highly dependent on the quality and suitability of data (García-Gil et al., 2019). Unfortunately, data gathering is rarely perfect, and data corruption often occurs (Wang et al., 2018). The identification of the multi-source SWC data quality or measurement error is not an easy task. This limitation can instead create extra uncertainties in DA systems (Kisekka et al., 2015). Third, it is the diversity of scenarios contained in prior data rather than its volume that is more decisive for the generalization ability of machine learning methods (Wang et al., 2020). Direct data fusion without screening may instead induce accidental correlations in learning algorithms, thereby diminishing their generalization ability (García et al., 2016). To avoid the overloaded monitoring cost due to redundant data, it is essential to develop a framework to assess the worth of alternative sampling strategies prior to data collection.

The data worth, sometimes called the data information content or data impact, of a design is often defined as its individual capacity to reduce uncertainty associated with a prediction goal or to maximize some related measure of data utility. Over the past decades, two main types of sophisticated DW analysis frameworks have been proposed to identify the most informative monitoring strategy in hydrology, namely, one type based on sensitivity analysis (Dausman et al., 2010; Fienen et al., 2010; Hill and Tiedeman, 2006) and the other within a fully Bayesian framework (Dai et al., 2016; Neuman et al., 2012; Nowak et al., 2012). The former approaches are computationally fast, but these methods require model calibration and assume linear models (Finsterle, 2015). The latter methods are derived based on the law of the total possibility without assumptions of the model and of the distributions of observations and model parameters. Nevertheless, both well-established frameworks are predicated on the availability of the underlying physical models. For example, Man et al. (2016) evaluated the expected value of alternative SWC sampling strategies with respect to the estimation of soil hydraulic parameters in the HYDRUS-1D model, while Finsterle (2015) examined the worth of datasets potentially applicable to the calibration of geothermal reservoir models. Within such parametric data worth analysis frameworks, however, the strong nonlinearity of soil water problems (Lannoy et al., 2006; Leube et al., 2012; Yeh et al., 1985) and the prevalence of model structural errors (Zhang et al., 2019) are highly likely to lead to biased data worth (Wang et al., 2018, 2020). Ultimately the reliability of the optimal design of monitoring networks based on such evaluations is greatly compromised.

Fortunately, the superiority of data-driven algorithms in handling nonlinearities and structural errors in unsaturated flow has been well demonstrated in our previous studies (Wang et al., 2021a, b). With the explosive growth of big data, how to evaluate the worth of multi-source data in this new data-mining approach is becoming a critical issue. Several recent works in the field of statistical learning have bloomed in identifying and removing irrelevant and redundant information from big data, such as feature selection (Chandrashekar and Sahin, 2014; Hall, 1999) and instance reduction (Amro et al., 2021; Olvera-López et al., 2010). To the best of our knowledge, few studies have systematically evaluated the worth of future observations regarding the construction of fully data-driven models prior to data gathering. As a follow-up study to Wang et al. (2018, 2020), one major contribution of this study is the first embedding of a purely data-driven model into the Bayesian data worth analysis framework, referred to as the nonparametric data worth analysis (NP-DWA). Similarly to traditional DW analysis, the proposed NP-DWA consists of prior, posterior and pre-posterior stages (Dai et al., 2016). The pre-posterior analysis evaluates the anticipated worth of future observations regarding the construction of purely data-driven models, for which possible distributions are predicted in advance by conditioning on prior data.

There is a consensus in the field of statistics that “the highest accuracy results that an inductive learning system can achieve depend on the quality of data and the appropriate selection of a learning algorithm for the data” (Pechenizkiy et al., 2006). In other words, once the algorithm is specified, the significance of data noise to learning accuracy as almost the only factor should not be overlooked. Considering the powerful ability of dealing with observational noises of the ensemble Kalman filter (EnKF) (Hamilton et al., 2017; Li et al., 2018), another innovation of this study is the introduction of EnKF into our NP-DWA framework. In conventional DW analysis, the worth of data is primarily embodied in its ability to be utilized or calibrated to adjust physical parameters (Dai et al., 2016; Finsterle, 2015; Man et al., 2016). In the proposed NP-DWA, nevertheless, future observations are first used to construct data-driven models in the forecast step and then are sequentially assimilated with the Kalman update in the analysis step. Ultimately, its combined capacity to reduce system uncertainties in these two ways is defined as its worth. Furthermore, as a typical sequential DA scheme, the EnKF facilitates the dynamic models and also its hyperparameters to be updated in real time so the data utility in relation to the modeling system can be detected instantaneously. Eventually, the sampling scheme can be dynamically adjusted to save the monitoring and analysis costs.

Most previous studies are based on synthetic cases, and data worth analysis in the context of dynamically evolving soil moisture profiles was still poorly studied in a real-world case. For nonlinear problems, nevertheless, the estimation variance and more sophisticated measures of data utility depend on the actual values of measurements, which are still unknown prior to collection (Leube et al., 2012). It will be more convincing to investigate the data worth regarding the reconstruction of fully data-driven models under real-world cases for unsaturated flow. With the aid of observed data retrieved from three typical stations with different climate regimes, we aim to shed light on the following questions: (1) as opposed to the traditional way of utilizing data (to calibrate physical parameters), is the worth of observations capable of being accurately quantified by NP-DWA in this new, purely data-driven approach? (2) Given multiple prediction objectives, how does the DW (in the form of various indices) evolve under different hydrometeorological conditions in the determination of fully data-driven soil moisture dynamics? (3) How does the proposed NP-DWA respond to the presence of multiple levels of data noise? We intend for this study to provide guidance in the design of future monitoring strategies within the fully data-driven soil water flow models for real-world problems.

The remainder of this paper is organized as follows: Sect. 2 first summarizes the experimental data and methods. Thereinto, the principles of Bayesian DW analysis, nonparametric DA and the hybrid framework are detailed. Section 3 presents the results, and a discussion is contained in Sect. 4. Finally, conclusions are outlined in Sect. 5.

2 Methodology

In Wang et al. (2021a), a nonparametric sequential data assimilation scheme (Kalman-GP) has been proposed based on the filtering equations of EnKF and data-driven modeling with GP. On top of that, this paper further develops a nonparametric data worth analysis framework to assess the potential worth of future observations in the reconstruction of dynamical soil water flow models prior to data collection. In our proposed NP-DWA framework, the mapping from the input vector x and state variable of interest y is approximated by the GP regression model. In this study, the target output y is the soil moisture, while the entries of x include the corresponding observation time, depth, daily precipitation and air temperature. Thus, here, the dimension of input x is d=4. Similarly to the Kalman-GP in Wang et al. (2021a), EnKF is implemented to update the forecast yield by GP models by assimilating real-time observations of soil moisture. On the one hand, this fusion can effectively reduce the risk of unreasonable spatio-temporal interpolation in GP models, ultimately enhancing the robustness of such purely data-driven models. On the other hand, by combining with the Kalman update, the forecast cross-covariance between the state and the predictions corresponding to available observations constrained the otherwise high error covariances of state variables at unobserved locations, which resulted in a significantly reduced uncertainty for this hybrid method relative to GP alone. More details can be found in Wang et al. (2021a). In addition, to demonstrate the generalizability of the NP-DWA, other machine learning algorithms (including support vector machine and random forest) and DA schemes (i.e., particle filter) are also employed to replace GP and EnKF, respectively. The specific details of these algorithms are not presented here, but the corresponding results are provided in the Supplement.

2.1 Construction of GP dynamic models

As stated in Wang et al. (2021a), N GP models are constructed independently at each time step in a sequential manner as new data are recorded. N is the ensemble size. Thereinto, the observed time series is corrupted by the prescribed observation noises satisfying Gaussian distribution to obtain N sets of training data. Here, we only take the procedure of building the mth GP model at any time t=k as an example, i.e., ${\mathrm{GP}}_k^m$ (m=1, 2, …, N). The superscript m will be omitted below for ease of expression.

At t=k, the input ${\mathbf{X}}_{\mathrm{1}:(k-\mathrm{1})}$ and output ${\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}$ of training data in GP_k can be expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\mathbf{X}}_{\mathrm{1}:(k-\mathrm{1})}=& \left[{\mathit{x}}_{\mathrm{1}}^{\mathrm{1}},{\mathbf{x}}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },{\mathit{x}}_{\mathrm{1}}^{q_{\mathrm{1}}},\mathrm{\dots },{\mathit{x}}_p^{\mathrm{1}},{\mathit{x}}_p^{\mathrm{2}},\mathrm{\dots },{\mathit{x}}_p^{q_p},\right.\\ & {\left.\mathrm{\dots },{\mathit{x}}_{k-\mathrm{1}}^{\mathrm{1}},{\mathit{x}}_{k-\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },{\mathit{x}}_{k-\mathrm{1}}^{q_{k-\mathrm{1}}}\right]}^T\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}=& \left[y_{\mathrm{1}}^{\mathrm{1}},y_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },y_{\mathrm{1}}^{q_{\mathrm{1}}},\mathrm{\dots },y_p^{\mathrm{1}},y_p^{\mathrm{2}},\mathrm{\dots },y_p^{q_p},\right.\\ & {\left.\mathrm{\dots },y_{k-\mathrm{1}}^{\mathrm{1}},y_{k-\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },y_{k-\mathrm{1}}^{q_{k-\mathrm{1}}}\right]}^T\end{array}\end{array}$$

where ${\mathbf{X}}_{\mathrm{1}:(k-\mathrm{1})}$ and ${\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}$ denote a collection of all available x and y from t=1 to (k−1), respectively; q_p denotes the number of available observations at t=p (p=1, 2, …, k−1). In this paper we assume that the number of available observations at each time step is identical, i.e., $q_{\mathrm{1}}=q_{\mathrm{2}}=\mathrm{\dots }=q_{k-\mathrm{1}}=q$. Hence the dimensions of matrix X_k−1 and vector y_k−1 are $q(k-\mathrm{1})\times d$ and q(k−1), respectively.

As defined in Rasmussen (2004) and Rasmussen and Williams (2006), a GP model can be fully specified by a mean function μ(x) and covariance function $k(\mathit{x},{\mathit{x}}^{\prime })$, i.e., $G\left(\mathit{x}\right)\sim N\left(\mathit{\mu }\right(\mathit{x}),k(\mathit{x}{\mathit{x}}^{\prime }\left)\right)$. In this study, a linear mean function and an anisotropic squared exponential covariance function are specified (Zhang et al., 2019) as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\mu }\left(\mathit{x}\right)={\mathit{\beta }}^T\mathit{x},\text{(4)}& {\displaystyle }& {\displaystyle }k(\mathit{x},{\mathit{x}}^{\prime })={\mathit{\sigma }}^{\mathrm{2}}\exp \left[-\sum _{l=\mathrm{1}}^d{\displaystyle \frac{{\left(x_l-x_l^{\prime }\right)}^{\mathrm{2}}}{{\mathit{\tau }}_l^{\mathrm{2}}}}\right],\end{array}$$

where β is a vector containing d linear coefficients, i.e., $\mathit{\beta }=\mathit{\{}{\mathit{\beta }}_{\mathrm{1}},{\mathit{\beta }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\beta }}_d\mathit{\}}$. σ² controls the marginal variance in the output, and $\mathit{\tau }=\mathit{\{}{\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\tau }}_d\mathit{\}}$ determines the dependence strength in each of the component directions of x.

Then, the hyperparameters of the GP_k, $\mathrm{\varnothing }=\mathit{\{}\mathit{\beta },{\mathit{\sigma }}^{\mathrm{2}},\mathit{\tau }\mathit{\}}$, can be inferred from the training datasets $\mathit{\{}{\mathbf{X}}_{\mathrm{1}:(k-\mathrm{1})},{\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}\mathit{\}}$ via log marginal likelihood maximization:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}L=& \log p\left({\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}\mathrm{|}{\mathbf{X}}_{\mathrm{1}:(k-\mathrm{1})},\mathrm{\varnothing }\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}-\mathit{\mu }\right)}^T\\ & {\mathbf{\Sigma }}^{-\mathrm{1}}\left({\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}-\mathit{\mu }\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\log \left|\mathbf{\Sigma }\right|-{\displaystyle \frac{n}{\mathrm{2}}}\log \mathrm{2}\mathit{\pi },\end{array}\end{array}$$

where μ denotes the prior mean vector with the dimension of q(k−1), and Σ denotes the covariance matrix whose elements in the ith row and jth column constitute ${\mathrm{\Sigma }}_{ij}=k({\mathit{x}}_i,{\mathit{x}}_j)$ (i=1, 2, …, q(k−1); j=1, 2, …, q(k−1)). The GPML MATLAB toolbox (version 4.2), as documented in Rasmussen and Williams (2006), was adopted for GP inference in this study (http://www.gaussianprocess.org/gpml/code/matlab/doc/, last access: 15 July 2023).

In this study, the entire soil moisture profile at t=k is expected to be forecasted. Assuming that the total number of nodes of the vertical one-dimensional soil profile is N_n, then the input at the current time step is ${\mathbf{X}}_k^{*}=[{\mathit{x}}_k^{\mathrm{1}},{\mathit{x}}_k^{\mathrm{2}},\mathrm{\dots },{\mathit{x}}_k^{N_n}{]}^T$ with the dimension of N_n×d. The corresponding output vector ${\mathit{y}}_k^{\mathrm{f}}$ with the dimension of N_n can be calculated as follows:

$$\begin{array}{}\text{(6)}& {\mathit{y}}_k^{\mathrm{f}}={\mathit{\mu }}^{*}+{\mathbf{\Sigma }}^{{*}^T}{\mathbf{\Sigma }}^{-\mathrm{1}}\left({\mathit{y}}_{\mathrm{1}:(k-\mathrm{1})}-\mathit{\mu }\right),\end{array}$$

where μ^* denotes the prior mean of $\mathit{\mu }\left({\mathbf{X}}_k^{*}\right)$ with the dimension of N_n, and Σ^∗ is calculated as ${\mathbf{\Sigma }}_i^{*}=k({\mathit{x}}_k^i,{\mathit{x}}_j)$ (i=1, 2, …, N_n; j=1, 2, …, q(k−1)).

As a collection of ${\mathit{y}}_k^{\mathrm{f}}$ from N GP models, the resultant forecasted state vector ${\mathit{Y}}_k^{\mathrm{f}}$ at t=k can be represented as follows:

$$\begin{array}{}\text{(7)}& {\mathbf{Y}}_k^{\mathrm{f}}={\left[{\mathit{y}}_{k,\mathrm{1}}^{\mathrm{f}},{\mathit{y}}_{k,\mathrm{2}}^{\mathrm{f}},\mathrm{\dots },{\mathit{y}}_{k,m}^{\mathrm{f}},\mathrm{\dots },{\mathit{y}}_{k,N}^{\mathrm{f}}\right]}^T,\end{array}$$

where ${\mathit{y}}_{k,m}^{\mathrm{f}}$ denotes the forecasted state vector of interest for ${\mathrm{GP}}_k^m$ (m=1, 2, …, N); the dimension of ${\mathbf{Y}}_k^{\mathrm{f}}$ is N×N_n; the superscript “f” denotes forecast.

2.2 The Kalman update in a nonparametric data assimilation scheme

In the analysis step of the EnKF, for any ensemble member m at t=k, the state vector can be updated by combing GP model predictions and observations ${\mathit{d}}_k^{\mathrm{obs}}$:

$$\begin{array}{}\text{(8)}& {\mathit{y}}_{k,m}^{\mathrm{a}}={\mathit{y}}_{k,m}^{\mathrm{f}}+{\mathbf{K}}_k\left({\mathit{d}}_{k,m}^{\mathrm{obs}}-\mathbf{H}{\mathit{y}}_{k,m}^{\mathrm{f}}\right),\end{array}$$

where ${\mathit{y}}_{k,m}^{\mathrm{a}}$ denotes the improved estimates for realization m by conditioning on the observed information at t=k; H is the observation operator with the dimension of q×N_n, which represents the relationship between the state and observation vectors; the superscript “a” indicates analysis, and ${\mathit{d}}_{k,m}^{\mathrm{obs}}$ with the dimension of q denotes the observation vector at t=k for the mth ensemble member of ${\mathit{d}}_k^{\mathrm{obs}}$. It should be emphasized that the relationship between observations at t=k, ${\mathit{d}}_k^{\mathrm{obs}}$ and their true values ${\mathit{y}}_k=[y_k^{\mathrm{1}},y_k^{\mathrm{2}},\mathrm{\dots },y_k^q{]}^T$ can be expressed as follows:

$$\begin{array}{}\text{(9)}& {\mathit{d}}_k^{\mathrm{obs}}={\mathit{y}}_k+{\mathit{\epsilon }}_k,\end{array}$$

where ε_k with the dimension of N×q represents the measurement error vector which is assumed to be zero-mean Gaussian with R_k; R_k denotes the error covariance matrix of the observations with the dimension of q×q.

The Kalman gain at t=k, K_k, with the dimension of N_n×q can be defined as follows:

$$\begin{array}{}\text{(10)}& {\mathbf{K}}_k={\mathbf{C}}_k^{\mathrm{f}}{\mathbf{H}}^T{\left({\mathbf{HC}}_k^{\mathrm{f}}{\mathbf{H}}^T+{\mathbf{R}}_k\right)}^{-\mathrm{1}},\end{array}$$

where ${\mathbf{C}}_k^{\mathrm{f}}$ with the dimension of N_n×N_n is the covariance matrix of the state vector at t=k, which can be approximated as follows:

$$\begin{array}{}\text{(11)}& {\mathbf{C}}_k^{\mathrm{f}}\approx {\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _{m=\mathrm{1}}^N\left\{\left[{\mathit{y}}_{k,m}^{\mathrm{f}}-〈{\mathbf{Y}}_k^{\mathrm{f}}〉\right]{\left[{\mathit{y}}_{k,m}^{\mathrm{f}}-〈{\mathbf{Y}}_k^{\mathrm{f}}〉\right]}^T\right\},\end{array}$$

where $\langle {\mathbf{Y}}_k^{\mathrm{f}}\rangle$ denotes the ensemble mean of ${\mathbf{Y}}_k^{\mathrm{f}}$.

2.3 Nonparametric data worth analysis framework

Following the methodologies of Neuman et al. (2012) and Dai et al. (2016), data worth analysis of future monitoring networks within the aforementioned NP-DWA framework also consists of three stages. The whole workflow of the NP-DWA framework is depicted in Fig. 1.

Figure 1 The workflow of nonparametric data worth analysis framework coupled with ensemble Kalman filter (EnKF).

2.3.1 Prior stage

At the prior stage ($\mathrm{0}<t\le T_p$), the integration of GP dynamic models and EnKF with an ensemble size of N=N₁ is implemented to sequentially train and assimilate the prior data via Eqs. (1)–(11). Here, all available prior datasets from t=0 to t=T_p are denoted as vector $\mathit{A}={\mathit{y}}_{\mathrm{1}:T_p}={\mathit{d}}_{\mathrm{1}:T_p}^{\mathrm{obs}}$ with the dimension of qT_p, while the corresponding GP input is denoted as matrix ${\mathbf{X}}_{\mathrm{1}:T_p}$ with the dimension of qT_p×d. Then, a set of N₁ hypothetical observations can be generated, denoted as ${\mathbf{B}}_{k,i}={\mathbf{H}}_k{\mathit{y}}_{k,i}^{\mathrm{f}}$ ($k=T_p+\mathrm{1},T_p+\mathrm{2},\mathrm{\dots },T_{\mathrm{t}}$; i=1, 2, …, N₁), via Eq. (6). T_t is the total simulation time. Moreover, prior prediction statistics (mean and covariance) of the posterior vector Y_k, i.e., E(Y|A) and Cov(Y|A), can be yielded conditional on {A}, which can be denoted as E₁ and C₁, respectively, for the sake of simplicity.

2.3.2 Pre-posterior stage

At the pre-posterior stage ($T_p+\mathrm{1}<t\le T_{\mathrm{t}}$), for each possible data B_k,i at t=k, N₂ realizations satisfying a Gaussian distribution are further generated. The ensemble mean is the value of B_k,i, while the variance is the measurement error. Since this method is recursive, the time index k is omitted in the following equations. Then, the integration of GP models and EnKF is again implemented through a set of N₂ Monte Carlo realizations for each of the N₁ hypothetical observations. This allows us to calculate prediction statistics of the posterior state vector Y_ij (i=1, 2, …, N₁; j=1, 2, …, N₂), i.e., E(Y_i|A, B_i) and Cov(Y_i|A, B_i), conditional on {A, B_i}. Finally, quantities E_B|AE(Y|A, B), E_B|ACov(Y|A, B) and Cov_B|AE(Y|A, B) can be yielded by averaging over the collection of N₁×N₂ realizations. It should be emphasized that E_B|AE(Y|A, B) and E_B|ACov(Y|A, B) represent the pre-posterior prediction mean and uncertainty after the addition of future possible data B, which can be denoted as E₂ and C₂, respectively.

To quantify the expected data worth of potential measurements, three commonly considered information metrics, namely the trace (Tr), Shannon entropy difference (SED) and relative entropy (RE), are introduced in this study. Tr and SED offer scalar indices to measure the decrease in variance and covariance, respectively, while the RE comprehensively quantifies both mean and covariance effects.

(1) Trace

As a scalar indicator (Dai et al., 2016), Tr quantifies the DW in terms of variance reduction as follows:

$$\begin{array}{}\text{(12)}& \mathrm{Tr}=\mathrm{Tr}\left({\mathbf{C}}_{\mathrm{1}}\right)-\mathrm{Tr}\left({\mathbf{C}}_{\mathrm{2}}\right),\end{array}$$

where Tr(*) denotes the trace (sum of the diagonal entries) of a matrix.

(2) Shannon entropy difference

According to Shannon (1949), the Shannon entropy (SE) of PDF p(x) can be defined as follows:

$$\begin{array}{}\text{(13)}& \mathrm{SE}\left(p\right)=-\int p\left(x\right)\ln p\left(x\right)\mathrm{d}x,x\in \mathit{ϵ}R.\end{array}$$

The SED between the prior and pre-posterior PDFs can also be considered to quantify the information content extracted from additional observations. Assuming that these two PDFs are both Gaussian in the EnKF model, the SED can be expressed in terms of covariance reduction (Xu, 2007) as follows:

$$\begin{array}{}\text{(14)}& \mathrm{SED}={\displaystyle \frac{\ln \det \left({\mathbf{C}}_{\mathrm{1}}\right)}{\mathrm{2}}}-{\displaystyle \frac{\ln \det \left({\mathbf{C}}_{\mathrm{2}}\right)}{\mathrm{2}}}={\displaystyle \frac{\ln \det \left({\mathbf{C}}_{\mathrm{1}}{\mathbf{C}}_{\mathrm{2}}^{-\mathrm{1}}\right)}{\mathrm{2}}},\end{array}$$

where det(*) is the determinant of a matrix.

(3) Relative entropy

Similar to the SED, the RE also provides a measure of the information content of the pre-posterior PDF with respect to the prior PDF. In addition to uncertainty reduction, the influence of future data on the mean behavior of PDFs is considered (Singh et al., 2013; Zhang et al., 2015). Considering that the prior and pre-posterior PDFs are n-dimensional Gaussian functions, the RE can be defined as follows:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathrm{RE}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathbf{E}}_{\mathrm{2}}-{\mathbf{E}}_{\mathrm{1}}\right)}^T{\mathbf{C}}_{\mathrm{1}}^{-\mathrm{1}}\left({\mathbf{E}}_{\mathrm{2}}-{\mathbf{E}}_{\mathrm{1}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\ln \det \left({\mathbf{C}}_{\mathrm{1}}{\mathbf{C}}_{\mathrm{2}}^{-\mathrm{1}}\right)+\mathrm{Tr}\left({\mathbf{C}}_{\mathrm{2}}{\mathbf{C}}_{\mathrm{1}}^{-\mathrm{1}}\right)-n\right].\end{array}\end{array}$$

Finally, the expected DW of B_k can be estimated in the form of the above three indices prior to data gathering. Similar procedures are repeated until the final time t=T_t is reached.

2.3.3 Posterior stage

At the posterior stage ($T_p+\mathrm{1}<t\le T_{\mathrm{t}}$), the available actual dataset B^′ is incorporated into the GP training datasets and assimilated in a sequential manner. The actual mean and covariance of posterior state vector Y, i.e., $E(\mathbf{Y}\mathrm{|}\mathbf{A},{\mathbf{B}}^{\prime })$ and $\mathrm{Cov}(\mathbf{Y}\mathrm{|}\mathbf{A},{\mathbf{B}}^{\prime })$, respectively, are obtained conditional on $\mathit{\{}\mathbf{A},{\mathbf{B}}^{\prime }\mathit{\}}$. The reference data worth in the form of the various indices can be calculated via Eqs. (12)–(15), where E₂ and C₂ are replaced with $E(\mathbf{Y}\mathrm{|}\mathbf{A},{\mathbf{B}}^{\prime })$ and $\mathrm{Cov}(\mathbf{Y}\mathrm{|}\mathbf{A},{\mathbf{B}}^{\prime })$, respectively.

3 Description of experimental data and model setup

3.1 Data sources and site description

Three typical sites, namely Falkenberg (52.1669^∘ N, 14.1241^∘ E), Cape_Charles_5_ENE (37.2907^∘ N, 75.9270^∘ W – hereafter referred to as Cape), and DAHRA (15.4035^∘ N, 15.4320^∘ W), were selected from the International Soil Moisture Network (ISMN, 2023) to evaluate the performance of the proposed NP-DWA framework under different soil types and climatic regimes. According to the dominant fraction of clay, silt and sand for two layers (topsoil: 0.0–0.3 m, subsoil: 0.3–1.0 m) provided by the ISMN, we use the United States Department of Agriculture (USDA) soil texture classification and classified the soil at three sites. The soil at Falkenberg is sandy loam, and the DAHRA soil is loamy sand. The topsoil and subsoil at Cape are clay loam and loamy clay, respectively. At these three sites, the in situ volumetric SWC was operationally measured with the TRIME-EZ (IMKO), the Stevens Hydraprobe II SDI-12 (Stevens Water Inc.) and the ThetaProbe ML2X (Delta-T Devices) instruments, respectively. The measurement depths were (1) 0.08, 0.15, 0.30, 0.45, 0.60 and 0.90 m at the Falkenberg site; (2) 0.05, 0.10, 0.20, 0.50 and 1.00 m at the Cape site; and (3) 0.05, 0.10, 0.50 and 1.00 m at the DAHRA site. The measurement error was artificially specified to be 0.02 cm³ cm⁻³ unless otherwise specified.

Apart from soil water measurements at different depths, the daily precipitation and air temperature at the height of 2 m were obtained from the ISMN. At each site, 200 d time series (from 15 January and 2 August 2005 at the Falkenberg site, from 24 April to 9 November 2004 at the Cape site and from 9 April to 25 October 2011 at the DAHRA site) were collected in this study, as shown in Fig. 2. Having a continental climate, the Falkenberg receives frequent but less intense precipitation during the simulation period. The Cape has a humid subtropical climate with the highest rates of rainfall among the three sites, and there were a few rainstorm events during the study period (e.g., up to 150 mm d⁻¹ on 8 September 2014). The region of DAHRA has a tropical climate with well-defined dry and humid seasons. The early stage of the simulation is in its dry season, with little to no rainfall. The late stage is in its humid season when frequent but less intense rainfall events occur and the daily average air temperature is about 30 ^∘C.

Figure 2 The temporal evolutions of soil moisture at various depths, daily rainfall and mean daily air temperature (AT) at 2 m height at Falkenberg, Cape and DAHRA, respectively. Note that the red area indicates the pre-posterior or posterior stage.

3.2 Model simulation setup and case design

The key parameters of this study are summarized in Table 1. Each site is represented by a one-dimensional soil column with a height of 1 m, which is discretized into 2 cm grids with the local refinement of 1 cm monitoring-depth intervals, i.e., z=0.15 m and 0.45 m at the Falkenberg site and z=0.05 m at the Cape and DAHRA sites. At each time step, N₁=50 GP-based dynamic models of unsaturated flow are constructed. The GP model input x includes the observation time, depth, daily precipitation and air temperature, while the output is the corresponding soil moisture. The state vector y comprises the soil moisture for all nodes at each site, and the trained and assimilated observations d^obs refer to the available soil moisture at all observed depths (as described in Sect. 3.1).

Table 1 The summary of key parameters.

Table 2 The summary of designed test cases and main characteristics.

We illustrate our approach based on a set of real-world test cases, as listed in Table 2. The performances of the three indices, namely, Tr, SED and RE, in data worth quantification are compared at all three sites. In all test cases of this study, prior data for training GP includes the soil water content at all observed depths during the prior stage (from t=1 to t=T_p), i.e., z=0.08, 0.15, 0.30, 0.45, 0.60 and 0.90 m at Falkenberg; z=0.05, 0.10, 0.20, 0.50 and 1.00 m at Cape; and z=0.05, 0.10, 0.50 and 1.00 m at the DAHRA. At the pre-posterior or posterior stage in the NP-DWA, the worth of potential observations regarding the retrieval of three quantities of interest, including ${\mathit{\theta }}_{\mathrm{1.00}}^{\mathrm{ave}}$, ${\mathit{\theta }}_{\mathrm{0.60}}^{\mathrm{ave}}$ and ${\mathit{\theta }}_{\mathrm{0.30}}^{\mathrm{ave}}$, is evaluated. Here, ${\mathit{\theta }}_{\mathrm{1.00}}^{\mathrm{ave}}$, ${\mathit{\theta }}_{\mathrm{0.60}}^{\mathrm{ave}}$ and ${\mathit{\theta }}_{\mathrm{0.30}}^{\mathrm{ave}}$ represent the average soil moisture in the top 1.00, 0.60 and 0.30 m, respectively. A comparison among cases TC1, TC2 and TC3 is designed to investigate the data worth response of surface (θ_S), middle (θ_M) and deep (θ_D) SWC regarding the above different prediction objectives. Specifically, θ_S refers to soil moisture at z=0.08 m at Falkenberg and at z=0.05 m at Cape and DAHRA, respectively. θ_M refers to soil moisture at z=0.45 m at Falkenberg and at z=0.50 m at Cape and DAHRA, respectively. θ_D refers to soil moisture at z=0.90 m at Falkenberg and at z=1.00 m at Cape and DAHRA, respectively. The prior datasets entering these cases comprise SWC, daily precipitation and air temperature over the first 80 d, as shown in the gray areas of Fig. 2. The subsequent 20 d data (red areas in Fig. 2) are augmented as additional data for reference DW assessment in the posterior stage.

As stated in Pechenizkiy et al. (2006) and Zhu and Wu (2004), the maximum accuracy of statistical learning algorithms mainly depends on the quality of training data in addition to the inherent bias in the algorithm itself. In other words, the magnitude and accuracy of the expected worth of driving data in machine-learning-based DA may be closely related to the noise level. Thus, two additional test cases (TC4 and TC5) are considered to evaluate the performance of the proposed NP-DWA framework under different measurement errors. The soil moisture measurement error variance values of 0.01² and 0.04² are artificially specified in TC4 and TC5, respectively, to be compared to a value of 0.02² in TC1-1.

Figure 3 The probability distributions (dotted curved line) of potential observation realizations, as well as their mean (dotted vertical line) and the corresponding actual soil water content (SWC) observation (solid line) in the surface, middle and deep layers on the 81st, 90th and 99th days at Falkenberg, Cape and DAHRA, respectively.

Moreover, test cases TC6 and TC7 differ from test case TC1-1. These test cases are designed to investigate the influence of the prior data content on data worth analysis, which facilitates the determination of the required prior information content to ensure the accuracy of data worth assessment. The 80 d prior data in test case TC1-1 are reduced backward in time to 40 d in test case TC6 and augmented forward to 180 d in test case TC7. In addition, test cases TC1-1, TC8 and TC9 consider the composite DW of different combinations of monitoring schemes. The comprehensive contributions of the surface SWC jointly with the middle and/or deep ones are compared with its individual contribution.

3.3 Evaluation setup

To compare the relative differences in data worth estimation accuracy under the various test scenarios, the mean absolute percentage error (MAPE) between the expected and reference data worth in the form of Tr, SED and RE is defined as follows:

$$\begin{array}{}\text{(16)}& \mathrm{MAPE}={\displaystyle \frac{\mathrm{1}}{\mathrm{Tt}-T_p}}\sum _{k=T_p+\mathrm{1}}^{T_{\mathrm{t}}}\left|{\displaystyle \frac{{\mathrm{DW}}_k^{\mathrm{Expect}}-{\mathrm{DW}}_k^{\mathrm{Refer}}}{{\mathrm{DW}}_k^{\mathrm{Refer}}}}\right|,\end{array}$$

where ${\mathrm{DW}}_k^{\mathrm{Expect}}$ and ${\mathrm{DW}}_k^{\mathrm{Refer}}$ denote the expected and reference DW values, respectively, at time step t=k.

4 Results and discussions

4.1 Optimal monitoring location for the multiple predictive objectives (TC1,TC2 and TC3)

Figure 3 shows the probability distributions of the generated potential observation realizations, as well as their ensemble mean and the corresponding actual observations of the surface (θ_S), middle (θ_M) and deep (θ_D) soil moisture at three sites. Only the results on the 81st, 90th and 99th days are presented here. Overall, the N₁=50 potential realizations could capture the actual SWC observations with acceptable accuracy. Specifically, the forecasted middle SWC exhibited a considerably more robust capturing performance with sustained better proximity of potential and actual θ_M throughout the simulation period. This occurred especially pronouncedly at the Cape site. For example, both surface and deep layers at Cape may be at risk of a poor fit of potential observations to measurements (Fig. 3d and w), while the generated middle SWC is always fairly well approximated to the corresponding actual values in Fig. 3m–o.

Figure 4 The expected data worth of potential soil moisture observations in the surface, middle and deep layers in the form of trace (Tr), Shannon entropy difference (SED) and relative entropy (RE), respectively, regarding the retrieval of average soil moisture in the top 0.30, 0.60 and 1.00 m at three sites.

Based on the above potential observations, their expected data worth regarding the retrieval of ${\mathit{\theta }}_{\mathrm{1.00}}^{\mathrm{ave}}$, ${\mathit{\theta }}_{\mathrm{0.60}}^{\mathrm{ave}}$ and ${\mathit{\theta }}_{\mathrm{0.30}}^{\mathrm{ave}}$ can be quantified in the form of Tr, SED and RE, as depicted in Fig. 4. Meanwhile, for ease of analysis, Fig. 5 compares the covariance matrices of entire soil moisture profiles in the prior stage, pre-posterior stage, and posterior stage. Only the results from the 81st to the 90th day at Falkenberg are revealed here. It can be observed that, despite an overall increasing trend over time, the values of expected DW were prone to local spikes due to changes in the atmospheric boundary conditions such as rainfall. First of all, this general trend of increase should be attributed to the sequential augmentation of potential observations based on existing prior data, resulting in the cumulative values of DW over time. However, abrupt changes in external forcing, such as unexperienced rainfall events on the 88^th day at the Falkenberg, could trigger temporal extrapolation of statistical learning (Li et al., 2020; Minns and Hall, 1996; Xu and Valocchi, 2015), which in turn led to a surge in prior predictive uncertainty, i.e., C₁=Cov(Y|A) (the first column of Fig. 5). Fortunately, joint GP training and sequential assimilation of real-time potential observations can effectively lower the risk of such irrational extrapolation (Wang et al., 2021a, b), allowing these temporal mutations to be substantially attenuated at the pre-posterior stage (i.e., $C_{\mathrm{2}}=E_{B\mathrm{|}A}\mathrm{Cov}(Y\mathrm{|}A,B)$) (the second to the fourth column of Fig. 5). This uncertainty reduction brought about the fusion of additional data became significantly larger when external forcing encountered mutations, which ultimately led to the localized surge in DW during rainfall events. We recall that such DW surges induced by contrasting scenarios also occurred in traditional DW analysis based on physically motivated models. As stated in Wang et al. (2018), “If one model undergoes different and contrasting scenarios, the model structural error is likely to appear since model parameters updated or calibrated under one scenario have not been examined under another scenario”. Although in different ways, the expected DW of future monitoring strategies within both the traditional parametric and proposed nonparametric DA frameworks heavily depends on the coverage of additional data by available prior scenarios (Li et al., 2020).

Figure 5 The covariance matrices of soil moisture profiles from the 81st to the 90th day at Falkenberg before (column 1) and after potential (columns 2–4) and corresponding actual (columns 5–7) soil moisture observations in the surface, middle and deep layers were fused, respectively.

Figure 6 The MAPEs between expected and reference data worth in the form of Tr, SED and RE of potential soil moisture observations in the surface (S), middle (M) and deep (D) layers, respectively, regarding the retrieval of average soil moisture in the top 0.30, 0.60 and 1.00 m at three sites.

Moreover, Fig. 4 also suggests that the optimal observation depth shifted as the prediction target changed. As expected, the surface SWC θ_S produced higher Tr, SED and RE values regarding the estimation of ${\mathit{\theta }}_{\mathrm{0.30}}^{\mathrm{ave}}$. As the depth range of the average SWC to be estimated was extended downward, the data worth advantages of θ_M and θ_D began to emerge. Surprisingly, the potential middle SWC still exhibited a considerably higher superiority even in ${\mathit{\theta }}_{\mathrm{1.00}}^{\mathrm{ave}}$ estimation. In other words, the soil moisture in the middle layer has the most robust advantage in data worth. This may be due to the fact that the integration of surface or deep SWC only reduced the uncertainty within the corresponding depth ranges (the second and fourth columns of Fig. 5), whereas the augmentation of θ_M significantly decreased the covariance matrices of the entirety of the SWC profiles (the third column of Fig. 5). This selection result of the optimal monitoring location seemingly contradicts previous findings within the traditional parametric DW analysis where the surface observations with the largest temporal variation always produced the greatest data worth, as reported in Dai et al. (2016) and Wang et al. (2018). This discrepancy is likely to depend on the different mechanisms that characterize soil moisture dynamics in the vertical direction between the two approaches. The traditional parametric unsaturated-flow model follows the law-of-mass-conservation-based physical governing equations (i.e., the Richardson–Richards equation, Richards, 1931; Richardson, 1922), as well as their physical properties, to simulate the soil water infiltration process. The strongest time-varying nature of surface SWC was conducive to the effective updating of the physical parameters in EnKF, eventually generating the maximum data worth (Wang et al., 2018). However, the spatial prediction performance of data-driven methods substantially hinged on the similarity of data between different depths. Theoretically, there occurs an inherent delayed response of soil moisture profiles to rainfall events, which has been well-documented experimentally (Wierenga et al., 1986; Bresler et al., 1971; Vauclin et al., 1979). This causes the temporal changes in surface and deep SWC to be naturally asynchronous, thus rendering their representativeness in characterizing the whole soil moisture profile somewhat limited. Ultimately, the complete reliance on statistical and information-theoretic measures allowed the most representative middle SWC to establish the most robust superiority in DW.

It can also be seen from Fig. 4 that, when using different information indices (i.e., Tr, SED and RE) to quantify the data worth, the optimal observation location selected is identical regardless of soil textures and climatic regimes. This conclusion is generally in line with Wang et al. (2018) and Man et al. (2016). Furthermore, to quantify the data worth assessment accuracy, Fig. 6 depicts the MAPE between the expected and reference data worth in the form of Tr, SED and RE of alternative monitoring schemes at different depths. It can be observed that the surface SWC yielded the smallest MAPE when retrieving ${\mathit{\theta }}_{\mathrm{0.30}}^{\mathrm{ave}}$ regardless of the metric type. Nevertheless, during the estimation of ${\mathit{\theta }}_{\mathrm{0.60}}^{\mathrm{ave}}$ and ${\mathit{\theta }}_{\mathrm{1.00}}^{\mathrm{ave}}$, the expected data worth of θ_M more accurately and robustly approached the reference counterparts with overall smaller MAPEs. We recall that this ranking of DW estimation accuracy was exactly in line with the ranking of the magnitude of their expected DW in Fig. 4. To be specific, a comparison of Figs. 4 and 6 reveals that potential observations with a larger expected DW are prone to a higher DW estimation accuracy due to its more robust ability to imitate the actual observations (Fig. 3).

4.2 Effects of observation noise (TC1-1, TC4 and TC5)

Figure 7 shows the probability distributions of the potential observation ensemble, as well as their mean and the corresponding actual observations of the surface SWC under different SWC noise levels. Similarly, only the results of the 81st, 90th and 99th days are displayed. It can be observed that a higher noise level was not always detrimental but rather expanded the distribution width along the SWC axis and produced a flatter curve. The risk of failure in the generated realizations capturing the real observations was thus reduced. Even on the 81st day at Falkenberg, for example, the increase in SWC error variance from 0.01² to 0.04² facilitated a better agreement between the potential and actual surface soil moisture, as revealed in Fig. 7a, j and s. Similar phenomena can also be found via a comparison of Fig. 7e and n.

Figure 7 The probability distributions (dotted curved line) of potential observation realizations, as well as their mean (dotted vertical line) and the corresponding actual soil water content (SWC) observations (solid line) in the surface layer on the 81st, 90th and 99th days at three sites under different measurement error variances, these being 0.01² (TC4), 0.02² (TC1-1) and 0.04² (TC5), respectively.

Figure 8 shows the temporal evolution and time-averaged MAPE of the expected and reference data worth in the form of three information indices under various noise levels. Some interesting findings can be obtained. (1) Overall, the potential SWC data corrupted by a lower noise level yielded larger data worth with higher accuracy. (2) Nevertheless, the occurrence of rainfall events triggered a futile DW increase while also rendering the potential observations with appropriately magnified observation errors more valuable. For instance, a properly inflated observation error of 0.02² on the 88th day at the Falkenberg site resulted in a notably higher data worth than that of 0.01², as highlighted by the dashed ellipse boxes in Fig. 8a and d. Furthermore, this increase in data worth resulting from noise amplification was particularly evident in the form of Tr over the other two metrics, as depicted in Fig. 8a–c and d–i. At DAHRA, potential observations with an observation error of 0.02² even produced a significantly higher Tr value than that of 0.01² throughout almost the entire simulation period (Fig. 8c). (3) As opposed to Tr and SED indices focusing only on the system uncertainty (variance or covariance), the expected RE, as a comprehensive mean-covariance-type metric, was often more challenging to approach its reference counterparts with the largest MAPE at all sites, as shown in Fig. 8j–l.

Figure 8 The temporal evolution (a–i) and time-averaged MAPEs (j–l) of the expected and reference data worth in the form of Tr, SED and RE at three sites, respectively, under different measurement error variances, these being 0.01² (TC4), 0.02² (TC1-1) and 0.04² (TC5), respectively.

4.3 Effects of prior data content (TC1-1, TC6 and TC7)

Figure 9a–i depicts the temporal evolution of the expected and reference data worth of the surface SWC in cases TC6, TC1-1 and TC7 with the 40, 80 and 180 d prior data content, respectively. Under normal circumstances, an increase in the available prior data content inevitably entails a shrinkage in the DW of subsequent data due to the possibility of information redundancy. However, this seems to be valid only for a modest increase in prior data (from 40 to 80 d) within our NP-DWA framework. The substantial augmentation in available data content from 80 to 180 d instead resulted in a notably higher DW of the additional data (Fig. 9). Even more unexpectedly, this DW growth was prevalent across sites, regardless of the soil types and climatic regimes. To clarify this anomaly, Fig. 10 further shows the predicted covariance matrices of soil moisture profiles conditional on {A} in the prior stage and {A, B} in the pre-posterior stage in cases TC6, TC1-1 and TC7, respectively. Only the results from the 81st to the 90th day at Falkenberg are presented here. Our previous studies have demonstrated that, although the mean values of potential samples can approach actual observations well in fully (Wang et al., 2021a) or partially (Zhang et al., 2019) data-driven dynamical systems, their ensemble was apt to suffer from considerable uncertainty (Wang et al., 2021b). Unfortunately, augmented prior data, despite its potential to enrich available GP training scenarios, failed to prevent the non-convergence of N_e=50 GP samples. In contrast, the additional noise associated with prior data supplementation could exacerbate the increase in the prior prediction uncertainty (i.e., C₁), as illustrated by a comparison between the first three columns of Fig. 10. It should be highlighted that the fusion of B enabled a notable reduction in the pre-posterior uncertainties (i.e., C₂) in the data assimilation system to a comparable level (the last three columns of Fig. 10), even with different prior data content. The gradual widening of the gap between C₁ and C₂ eventually yielded the highest data worth with the maximum amount of prior data in test case TC7 (Fig. 9). This seems to warn us that uncontrolled expansion of big data within fully data-driven systems may not be beneficial. The adverse effects of extra noise may overshadow its original superiority in generalization capability. Access to high-quality and representative small data may constitute the key to the successful application of fully data-driven algorithms for reshaping soil moisture dynamics.

Figure 9 The temporal evolution (a–i) and time-averaged MAPEs (j–l) of the expected and reference data worth in the form of Tr, SED and RE at three sites for cases TC6, TC1-1 and TC7 with 40, 80 and 180 d prior data content, respectively.

Figure 10 The covariance matrices of soil moisture profiles from the 81st to the 90th day at Falkenberg in the prior and pre-posterior stages for cases TC6, TC1-1 and TC7 with 40, 80 and 180 d prior data content, respectively.

Furthermore, Fig. 9j–l depicts the time-averaged MAPE in the expected and reference data worth in cases TC6, TC1-1 and TC7, respectively. A comparison of Fig. 9a–i and j–l reveals some interesting findings: (1) similarly to the results in Sect. 4.1, the potential measurements with the largest expected (or reference) data worth in TC7 are apt to possess the highest estimation accuracy of data worth. (2) Local variations in data worth at different sites respond slightly differently to the augmentation of prior data content. For instance, even with 180 d of available historical data, the DW spike induced by the unexperienced rainfall event on the 88th day at the Falkenberg has not been eliminated or diminished (Fig. 9a, d and g). However, similar DW surges on the 82nd day at Cape were successfully mitigated as the amount of prior data content increased from 40 d (TC6) to 80 d (TC1-1) (Fig. 9e and h). This is because the prior data at Falkenberg, even if augmented to 180 day, did not cover the rainfall event on the 88th day (Fig. 2b), whereas the 80 d training data at Cape already included the scenario on the 82nd (Fig. 2d). These results agree with the conclusions reported in Wang et al. (2020) that the diversity of scenarios in the training data is more decisive than the data volume regarding the performance of data-driven methods. However, it is worth noting that the response of our nonparametric framework to the prior data augmentation is not in line with that of the physical-model-based DW analysis framework (Wang et al., 2018). In the latter, insufficient prior scenarios within the traditional framework may trigger unresolved model structural errors. When these lead to a deteriorating DW assessment performance, this cannot be compensated for by assimilating more prior data. In contrast, the elimination of dependence on physical governing equations shields our model-free DA schemes from model structural errors. Continuous enrichment of prior scenarios in the NP-DWA directly assures effective coverage of potential future scenarios, yielding the mitigation of DW local surges and an improved DW estimation accuracy. (3) Although inferior to Tr and SED, the estimation accuracy of RE is generally acceptable, especially when prior data is expanded to 180 d. This is certainly a remarkable improvement over the rather poor performance of RE in traditional parametric data worth analysis (Wang et al., 2018, 2020). This progress should be attributed to the radical abandonment of physical models in the NP-DWA, which prevented adverse effects of the high nonlinearity of soil water flow in the propagation of uncertainties from input to output (i.e., soil moisture in this study). Direct mapping from regular meteorological data to SWC facilitated the identification of the soil moisture covariance matrix from potential observations.

4.4 Effects of potential observational combinations (TC1-1, TC8 and TC9)

Figure 11a–i compares the expected (and reference) data worth of three combinations of potential observations at different depths at the three sites. It can be seen that the composite data worth of the alternative monitoring schemes exhibited an increasing pattern as the depth range of the observed SWC continued to expand downward. Nevertheless, the response of the different data worth indicators and study sites to this vertical expansion of potential observations varied slightly. Further integration of θ_D in TC9 did not cause a marked increase in Tr but yielded notably greater SED and RE values, especially at the DAHRA site (Fig. 11c, f and d). This is undoubtedly due to the extra consideration of the latter two indicators for the non-diagonal elements of the covariance matrix or/and the behavior of the mean. Moreover, the joint fusion of potential θ_S and θ_M failed to result in a sustained increase in Tr and RE at DAHRA while creating a significant increase in composite DW at the other two sites. This could be attributed to the sandy soil texture at DAHRA (with the fraction of sand up to 90 % and K_s=3.22 m d⁻¹), resulting in the almost synchronous responses of the SWC at z=0.05 m and 0.50 m to the atmospheric boundary conditions (Fig. 2e) and thus triggering possible data redundancy.

Figure 11 The temporal evolution (a–i) and time-averaged MAPEs (j–l) of the expected and reference data worth in the form of Tr, SED and RE at three sites for cases TC1-1 (surface soil moisture), TC8 (surface and middle soil moisture) and TC9 (surface, middle and deep soil moisture), respectively.

Figure 11j–l further show the estimation accuracy of expected data worth for the above three potential observation combinations. Surprisingly, the increase in the number of potential observations, while making it more difficult to capture actual SWC data, ends up significantly improving the accuracy of the data worth assessment. As shown in Fig. 11j–l, as more potential observations along the vertical direction were evaluated, the MAPEs between expected data worth and its reference counterparts decreased continuously. This phenomenon actually breaks the misconceptions about the data worth assessment accuracy in previous studies, i.e., that an excellent fit of potential observations is equivalent to high-precision estimates of the corresponding data worth. For the sake of explanation, Fig. 12 shows the predicted covariance matrices of soil moisture profiles in cases TC1-1, TC8 and TC9 from the 81st to the 90th day at Falkenberg conditional on {At}, {A, B} and $\mathit{\{}\mathbf{A},{\mathbf{B}}^{\prime }\mathit{\}}$, respectively. It can be found that, compared to TC1-1, which only reduces the uncertainties in the surface SWC, the integration of observations at multiple depths clearly reduces the uncertainties in the entirety of the SWC profiles to a considerably lower level. This ultimately facilitates better proximity between expected and reference covariance matrices, as revealed in the fourth and seventh columns in Fig. 12. The above results suggest that the accuracy of data worth assessment of potential observations does not only depend on their capacity to capture actual measurements but is also closely related to their correlation with the variable of interest. We recall that similar phenomena also exist in the preceding test cases. For example, the weaker correlation between surface SWC observations and ${\mathit{\theta }}_{\mathrm{1.00}}^{\mathrm{ave}}$ led to deterioration in the DW estimation performance with the largest MAPE values (Fig. 6a, d and g) even if the actual surface observations could be suitably reproduced (Fig. 3a–c). Therefore, to enhance the reliability of data worth assessment, a strategy wherein potential observations at multiple depths were simultaneously incorporated into existing DA systems was recommended in this study.

Figure 12 The covariance matrices of soil moisture profiles from the 81st to the 90th day at Falkenberg before (column 1) and after potential (columns 2–4) and actual observations (columns 5–7) for cases TC1-1 (surface soil moisture), TC8 (surface and middle soil moisture) and TC9 (surface, middle and deep soil moisture) were fused, respectively.

5 Conclusions

Conventional data worth analysis for soil water problems depends on physical dynamic models. Due to the widespread occurrence of model structural errors, it may lead to biased or wrong worth assessment. The strong nonlinearity of unsaturated flow further deteriorates the DW assessment performance in the retrieval of soil moisture profiles. This study proposed a nonparametric data worth analysis method within a fully data-driven modeling framework. The information extracted from real-time soil moisture data after GP training and Kalman update was quantified with three representative types of indicators, i.e., variance- (Tr), covariance- (SED) and mean-covariance-type (RE) indicators. With the aid of a series of real-world cases, the ability and challenge of the NP-DWA in terms of the variables of interest, spatial location, observation error and prior data content were assessed. The following conclusions were drawn:

1. The proposed NP-DWA framework enabled an accurate assessment of the data worth of potential observations regarding the reconstruction of purely data-driven soil water flow models prior to data collection. Similarly to traditional DW analysis based on physical models, the overall increasing trend of the DW from the sequential augmentation of additional observations within the proposed NP-DWA framework was also susceptible to interruptions by localized surges due to never-experienced atmospheric conditions. The difference is that this biased DW in the traditional parametric method is caused by model structural errors triggered by contrasting scenarios, which is difficult to be compensated for by assimilating more prior data, while the adverse effects of anomalous GP extrapolation in our NP-DWA could be suitably avoided by the enrichment of training scenarios in prior data. Moreover, the appropriate amplification of observational noise under extreme meteorological conditions also facilitated the alleviation of these biased estimates by enhancing the generalization capacity of dynamic models.

2. The optimal observation depth shifted as the prediction target varied. In contrast to the notably higher DW of surface SWC observations within the conventional DW analysis framework based on physical models, middle SWC observations tended to exhibit considerably higher robustness in the construction of model-free soil moisture dynamic models. This should be attributed to the ability of the SWC in the middle layer to effectively reduce the predictive uncertainty of the entirety of the soil moisture profiles due to its optimal representativeness. The inherent delayed response of soil moisture profiles to rainfall events allowed this advantage of middle SWC prevalent across sites, even becoming increasingly pronounced with increasing delay effect.

3. Although the addition of prior data content could greatly improve the estimation accuracy of the expected DW, the ensuing observation noise could substantially increase the uncertainty in a purely data-driven DA system, leading to potentially higher data worth of subsequent observations. Hence, high-quality and representative small data may be regarded as a better choice than unfiltered big data.

4. The performance of data worth assessment was jointly determined by the 3Cs, i.e., the capacity of potential observation realizations to capture actual observations, the correlation of potential observations with the predicted variables of interest, and the choice of DW quantitative indicators. Furthermore, the direct mapping from regular meteorological data to SWC in our nonparametric method facilitated the identification of the soil moisture covariance matrix (especially the cross-covariance) due to its alleviation of the high nonlinearity of soil water flow problems. Hence, satisfactory estimation accuracy could also be achieved even with covariance-related data worth metrics (i.e., the SED and RE).