Uncertainty estimation for a new exponential-filter-based long-term root-zone soil moisture dataset from Copernicus Climate Change Service (C3S) surface observations

Soil moisture is a key variable in monitoring climate and an important component of the hydrological, carbon, and energy cycles. Satellite products ameliorate the sparsity of field measurements but are inherently limited to observing the near-surface layer, while water available in the unobserved root-zone controls critical processes like plant water uptake and evapotranspiration. A variety of approaches exist for modelling root-zone soil moisture (RZSM), including approximating it from surface layer observations. While the number of available RZSM datasets is growing, they usually do not contain estimates of their uncertainty. In this paper we derive a long-term RZSM dataset (2002–2020) from the Copernicus Climate Change Service (C3S) surface soil moisture (SSM) COMBINED product via the exponential filter (EF) method. We identify the optimal value of the method's model parameter T, which controls the level of smoothing and delaying applied to the surface observations, by maximizing the correlation of RZSM estimates with field measurements from the International Soil Moisture Network (ISMN). Optimized T-parameter values were calculated for four soil depth layers (0–10, 10–40, 40–100, and 100–200 cm) and used to calculate a global RZSM dataset. The quality of this dataset is then globally evaluated against RZSM estimates of the ERA5-Land reanalysis. Results of the product comparison show satisfactory skill in all four layers, with the median Pearson correlation ranging from 0.54 in the topmost to 0.28 in the deepest soil layer. Temporally dynamic product uncertainties for each of the RZSM product layers are estimated by applying standard uncertainty propagation to SSM input data and by estimating structural uncertainties in the EF method from ISMN ground reference measurements taken at the surface and at varying depths. Uncertainty estimates were found to exhibit both realistic absolute magnitudes and temporal variations. The product described here is, to the best of our knowledge, the first global, long-term, uncertainty-characterized, and purely observation-based product for RZSM estimates up to 2 m depth.

1 Introduction

Soil moisture (SM) is an essential climate variable (ECV) that is crucial for understanding and modelling the Earth's climate and an important control of hydrological, energy, and carbon fluxes (GCOS, 2022; Dorigo et al., 2021a). Global monitoring of SM is necessary for a variety of applications such as meteorological modelling (Albergel et al., 2008); monitoring drought (Tobin et al., 2017); and modelling groundwater recharge (Bouaziz et al., 2020), runoff, and catchment response to storms (Brocca et al., 2010).

In situ SM measurements are considered to provide the most accurate SM data but can differ greatly in measuring equipment and usually lack estimates of their uncertainties (Dorigo et al., 2011). Widely distributed SM field measurements are available from centralized platforms such as the International Soil Moisture Network (ISMN) (Dorigo et al., 2021b). While being essential for satellite and model product calibration and validation, in situ measurements lack the spatial coverage necessary for large-scale applications, especially in the Global South (see Fig. A1; Dorigo et al., 2021a; Mishra et al., 2020). Quasi-global SM information is available from modelled and satellite products, but their spatial resolution is very coarse (usually tens to hundreds of square kilometres) and usually insufficient to resolve the significant spatio-temporal heterogeneity of SM, which poses challenges to large-scale monitoring (Brocca et al., 2010). Global land surface model products provide gap-free and long-term SM estimates at various depths and chosen time intervals but are computationally expensive and may depend on many auxiliary inputs that are not always available globally or with sufficient quality or resolution (Mishra et al., 2020; Albergel et al., 2008). In contrast, remote sensing retrievals are available only at satellite overpass times and are unreliable under various conditions, including frozen ground, dense vegetation, and radio frequency interference (RFI) (Gruber et al., 2019; Dorigo et al., 2017). Moreover, microwaves used for SM retrieval mainly contain information on water content in the surface layer, hampering their usability for studying or modelling processes in the soil root zone. Root-zone soil moisture (RZSM), often defined as the water present in the top metre of the soil column (Mishra et al., 2020; Baldwin et al., 2017; de Lange et al., 2008), is a component of the Global Climate Observing System (GCOS) ECV portfolio and a necessary variable for closing the water cycle (GCOS, 2016, 2022). RZSM also represents the water available for plant water uptake and thus affects evapotranspiration rates (Martens et al., 2017; Ford et al., 2014; Albergel et al., 2008) and plays a critical role in agricultural productivity forecasting (Wang et al., 2017) and drought monitoring (Vreugdenhil et al., 2022; Tobin et al., 2017).

The existing link between SM dynamics in the surface layer and the root zone (Albergel et al., 2008; Wang et al., 2017; Ford et al., 2014; Sure and Dikshit, 2019) allows for the estimation of RZSM from surface SM (SSM) observations via a variety of hydrological models. These include relatively simple two-layer approaches approximating RZSM as a function of SSM (Manfreda et al., 2014), compound process-based models requiring sophisticated parameter calibration (Bouaziz et al., 2020), and immensely complex and computationally expensive land surface models requiring many auxiliary inputs (Muñoz Sabater et al., 2021; Rodell et al., 2004). Satellite-based SSM observations can also be assimilated into a land surface model to produce estimates of RZSM with global coverage, as in the case of the SMAP L4 RZSM product (Reichle et al., 2017). An alternative, less complex approach that approximates RZSM solely from SSM estimates – and can thus be readily applied to satellite retrievals – is the so-called exponential filter (EF) method (Wagner et al., 1999; Albergel et al., 2008). In essence, the EF method approximates conditions in the root zone by smoothing and delaying SSM, which is generally characterized by greater fluctuations (Beck et al., 2009; Mahmood and Hubbard, 2007). Even though the coupling strength between the surface and root-zone layers decreases with depth (Mahmood and Hubbard, 2007; Ford et al., 2014; Mishra et al., 2020), and the skill of the method in predicting RZSM has been demonstrated to deteriorate accordingly (Paulik et al., 2014; Brocca et al., 2010; Sure and Dikshit, 2019), it is still widely used due to its relatively good performance and independence of ancillary inputs as well as its low computational cost and overall simplicity. However, the EF method is susceptible to prolonged data gaps in SSM data and thus requires an adequate number of input observations within a time interval consistent with the temporal scale of RZSM dynamics.

Regardless of the method used to derive RZSM estimates, most products do not provide information about the magnitude of random errors such as the standard deviation of their distribution, hereinafter referred to as uncertainties (Gruber et al., 2020). Two approaches have been proposed to characterize the time-variant quality of RZSM estimates derived with the EF method. The first approach, reported in Bauer Marschallinger (2018) and also utilized in this study, is a quality flag that is derived from the number of valid SSM estimates available within a specific time window preceding a specific EF-based RZSM estimate. The second approach, proposed by De Santis and Biondi (2018), uses the standard law of uncertainty propagation (Taylor, 1997) in order to characterize the random error variances of EF-based RZSM estimates. This approach takes into account the uncertainties in both the SSM input data and the EF model parameter but does not consider the model structural error in the EF method (Beven, 2005). The latter, due to the simplistic nature of the EF method and the limited surface–root-zone coupling, can also contribute significantly to the uncertainty budget and thus must not be neglected when characterizing product errors.

In this paper, we propose to estimate the model structural uncertainty in the EF using in situ measurements of surface and root-zone SM from the ISMN. We then use these estimates together with the law for the propagation of uncertainties (similar to De Santis and Biondi, 2018) to produce a global, fully error-characterized RZSM dataset for four soil layers (0–10, 10–40, 40–100, and 100–200 cm) between 2002 and 2020, taking C3S soil moisture as input to the model. While other EF-based datasets exist (e.g. the SMOS L4 product), they offer limited spatio-temporal coverage and lack quantitative uncertainty information (Al Bitar and Mahmoodi, 2020; Bauer-Marschallinger et al., 2018). The focus and novelty of this paper are in quantifying, rather than reducing, the EF model's known limitations by providing methodology for comprehensive uncertainty estimation for the EF method. Additionally, to the best of our knowledge, our dataset is, as yet, the longest available observation-based, error-characterized global RZSM product.

2 Datasets and data pre-processing

2.1 C3S surface soil moisture

Global input satellite surface observations were obtained from the Copernicus Climate Change Service (C3S) surface soil moisture COMBINED product v202012, hereinafter referred to as C3S SSM. C3S SSM is a merged product that combines satellite SSM retrievals from 4 active and 10 passive microwave sensors into a daily global dataset on a regular 0.25^∘ grid, expressed in volumetric units (m³ m⁻³) (C3S, 2020). Invalid retrievals due to frozen ground, dense vegetation, RFI, and other factors are masked out. Although the C3S product provides SSM data from 1978 onward, their quality and spatio-temporal coverage have increased significantly in more recent periods, when sensors measuring in frequency domains better suitable for SSM retrieval became available. Therefore, only C3S SSM data for the period 2001–2020 were used in this study. Note that data from the first year of this period were used only as the model adjustment period and not included in later analyses.

The uncertainty estimates provided for the merged SSM retrievals in the C3S SSM product were computed by means of triple collocation analysis (TCA) (Gruber et al., 2017). More specifically, (stationary) uncertainties were estimated for each satellite sensor separately and used to calculate the merging weights. Uncertainties in the merged SSM estimates were then calculated from the law for the propagation of uncertainties to account for the quality improvement due to the merging. Note that the distinctive life spans of the satellite missions therefore also lead to distinctive changes in the data quality of the merged product. These sudden changes in product uncertainty are hereinafter referred to as structural breaks (Preimesberger et al., 2021). As more and newer sensors provide better-quality retrievals, mean uncertainty values after each structural break typically decrease (Gruber et al., 2017). This is apparent, for example, in the shift in C3S SSM uncertainty values after the introduction of AMSR-E in 2002 (van der Schalie et al., 2017; Gruber et al., 2019).

C3S data are readily available from the Copernicus Climate Data Store (CDS), and detailed information on the C3S dataset and its underlying ESA CCI v5 merging algorithm can be found in the relevant documentation (C3S, 2020; Dorigo et al., 2021c).

2.2 Soil moisture field measurements

Field measurements for optimizing the model parameters of the EF method and for estimating its uncertainties were obtained from the International Soil Moisture Network (ISMN) for the period 2002–2020 (Dorigo et al., 2021b). Only data from sensors with a measuring depth ≤200 cm and internally flagged as reliable (Dorigo et al., 2013) were considered. Measuring depths of SM sensors placed vertically in a depth range, e.g. 10–40 cm, refer to their mean measuring depth. Data from multiple sensors installed at the same location and depth were averaged. ISMN data, typically available as hourly readings, were aggregated to mean daily values to match the temporal sampling of satellite observations. Furthermore, we only used ISMN stations where at least 100 data points concurrent with C3S SSM retrievals were available. Notably, approximately 80 % of the selected ISMN time series originate from North America and Europe (Fig. A1), and the availability of data declines with depth.

2.3 ERA5-Land soil moisture

ERA5-Land (E5L) is a multi-decadal climate reanalysis with an extensive portfolio of land variables computed by the assimilation of ERA5 atmospheric variables into the H-TESSEL land surface model (Muñoz Sabater et al., 2021). Modelled SM data are available for four depth layers (0–7, 7–28, 28–100, and 100–289 cm) on a regular 0.1^∘ grid and are accessible via the Copernicus Climate Data Store (CDS) (Muñoz Sabater, 2019, 2021). We used E5L for a product intercomparison with the RZSM product developed in this study, carried out for the period 2002–2020 within the Quality Assurance for Soil Moisture framework (QA4SM; https://qa4sm.eu, last access: 28 August 2023), which automatically resamples and matches observations of the compared datasets and delivers a wide range of validation metrics.

3 Methods

3.1 Exponential filter

The EF method (Wagner et al., 1999) relies on a simple two-layer water balance model where the only considered exchange between the surface layer and the reservoir below it is infiltration. The method assumes that the fluxes from the surface to the sub-surface layers are proportionate to the difference in SM content between both layers. In this study, we utilize the recursive formulation of the method (Albergel et al., 2008):

$$\begin{array}{}\text{(1)}& \mathrm{RZSM}\left(t_n\right)=\mathrm{RZSM}\left(t_{n-\mathrm{1}}\right)+K_n\cdot \left(\mathrm{SSM}\right(t_n)-\mathrm{RZSM}(t_{n-\mathrm{1}}\left)\right),\end{array}$$

where t_n and t_n−1 denote timestamps (in days) of the current and previous SSM observations, respectively. Conditions in the root zone are approximated by a weighted combination of the new input SSM observation and past model estimates, with more recent estimates receiving higher weights on a timescale defined by the method's only parameter T (temporal length, typically in days). Weights are controlled by the gain term K, which ranges from 0 to 1 and is calculated as follows:

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

At initialization, when no preceding estimates are available, the EF calculation is started with K₀=1 and RZSM(t₀)=SSM(t₀).

Temporal variability in the root zone is generally smaller than at the surface; hence the T value and its associated level of smoothing applied to the SSM data increase with depth (Wagner et al., 1999; Paulik et al., 2014; Wang et al., 2017; Beck et al., 2009; Mahmood and Hubbard, 2007). The optimal T value (T_opt, the value that leads to the best possible representation of RZSM at a certain location using the EF) has been related to differences in utilized SSM sensors (Bouaziz et al., 2020; Sure and Dikshit, 2019), SSM sampling frequency (Brocca et al., 2010; Pellarin et al., 2006), and land surface features (Albergel et al., 2008; de Lange et al., 2008). In particular, T acts as a conglomerate proxy for various environmental factors assumed to govern the infiltration process (e.g. soil texture, evapotranspiration, and climate), but past research on the importance of the exact driving factors is inconclusive and even contradictory (Wang et al., 2017; Bouaziz et al., 2020). To optimize the T parameter, numerous control factors have been tested (Bouaziz et al., 2020; Mishra et al., 2020; Stefan et al., 2021), and ever more sophisticated methods have been employed, including machine learning approaches (Grillakis et al., 2021). Other limitations of the method include generally poorer performance in arid zones and when soil texture is not homogeneous throughout the soil column (Ford et al., 2014; Yang et al., 2022).

Due to the high spatio-temporal heterogeneity of SM (Famiglietti et al., 2008) and its surface–root-zone coupling – and hence the difficulty in properly estimating the T parameter accurately – an uncalibrated value of T=20 has sometimes been used to describe all of the water content in the first 100 cm of the soil column (Wagner et al., 1999; de Lange et al., 2008). Results obtained by using a constant value T=20 were similar to those obtained with T values calibrated for soil texture (de Lange et al., 2008). Limited sensitivity of the EF to T due to different environmental factors was also observed by other studies, which supports choosing a single value for T_opt to represent a particular depth for large areas or even globally (Albergel et al., 2008; Brocca et al., 2010, 2011; Grillakis et al., 2021).

3.1.1 RZSM quality flags

Prolonged temporal data gaps will cause K to increase and may cause the EF to put excessive weight on new SSM input. In the extreme case, a very long data gap (whose duration depends on the chosen T value) can reset the EF to the initial state of K_n=1 and RZSM(t_n)=SSM(t_n) (see above). We run a 1-year adjustment period (2001) for K to reach an equilibrium state and utilize the EF quality flag (qflag) described in Bauer Marschallinger (2018) to avoid such re-initializations due to frequent and/or persistent data gaps. The qflag is recursively calculated for each RZSM estimate and reflects the availability of SSM input data in the preceding time period.

$$\begin{array}{}\text{(3)}& \mathrm{qflag}\left(t_n\right)=\left\{\begin{array}{ll}\mathrm{1}+\mathrm{qflag}\left(t_{n-\mathrm{1}}\right)\cdot e^{-\frac{t_n-t_{n-\mathrm{1}}}{T}},& \text{if SSM at}{t}_n\\ & \text{is available}\\ \mathrm{qflag}\left(t_{n-\mathrm{1}}\right)\cdot e^{-\frac{t_n-t_{n-\mathrm{1}}}{T}},& \text{if SSM at}{t}_n\\ & \text{is unavailable}\end{array}\right.\end{array}$$

The quality flag calculation is initialized with qflag(t_n)=1. A normalization factor of ${\sum }_{j=\mathrm{0}}^{\mathrm{\infty }}{\mathrm{e}}^{-\frac{j}{T}}$ is used to express the calculated flag values in percentages, with higher values indicating a greater density of SSM data available for calculation. If the quality flag falls below a T-specific threshold, RZSM estimates are masked out. The thresholds used here have been interpolated from those empirically determined by Bauer Marschallinger (2022) for a set of discrete T values (35 %, 40 %, 45 %, 50 %, 55 %, 60 %, 65 %, and 70 % for the T values 2, 5, 10, 15, 20, 40, 60, and 100, respectively). If input data are unavailable, but satisfactory data density has been achieved in the preceding days, the latest RZSM estimate is propagated forward until new input data become available, or the quality flag drops below its respective threshold. In the latter case, the output value is masked out. Importantly, even if new SSM input becomes available to the EF after prolonged data gaps, RZSM estimates derived from it remain masked until the qflag exceeds the aforementioned threshold again.

3.1.2 T-parameter optimization

We optimize T for a particular depth of the soil column by maximizing the correlation between the satellite-based RZSM estimates and the in situ measurements (Paulik et al., 2014; Grillakis et al., 2021). Satellite and in situ data are matched in space by means of the nearest-neighbour method. The impact of the spatial mismatch error between the large footprint of the satellite-based product and point-scale field measurement is mitigated by excluding time series that exhibit a correlation coefficient (Pearson's r) lower than 0.5 (Grillakis et al., 2021) or that are not statistically significant (p≥0.05).

EF calculations are repeated for T values of 1–100, and T_opt is selected for each of the available ISMN time series based on the highest correlation coefficient. We then group T_opt values based on the measurement depth of the respective in situ sensor into four bins corresponding to the RZSM target layers. These depth layers, chosen to be 0–10, 10–40, 40–100, and 100–200 cm, were defined to reflect those in common model-based RZSM products (Rodell et al., 2004; Muñoz Sabater et al., 2021). Finally, the median value of T_opt from each bin is chosen to compute a global RZSM product from the C3S SSM dataset.

A cross-validation is carried out to verify that T_opt values were not over-fitted to the local ISMN site conditions. Therefore, the sample set is randomly divided into five subsets of equal size (per bin), and then each of the subsets was used once to validate the method fit to the remaining four bins.

3.2 Uncertainty estimation

3.2.1 Baseline method

In De Santis and Biondi (2018), the standard law for the propagation of uncertainties is applied to the EF method. We use this approach as a baseline for our analyses. The recursive formulation of this baseline method is as follows:

$$\begin{array}{}\text{(4)}& \mathit{\sigma }\left({\mathrm{RZSM}}_n\right)=\sqrt{{\mathrm{\Delta }}_n^{\mathrm{2}}+\left({\displaystyle \frac{\partial {\mathrm{RZSM}}_n}{\partial T}}{)}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\right(T)},\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\mathrm{\Delta }}_n^{\mathrm{2}}=K_n^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\left({\mathrm{SSM}}_n\right)+(\mathrm{1}-K_n{)}^{\mathrm{2}}{\mathrm{\Delta }}_{n-\mathrm{1}}^{\mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\partial {\mathrm{RZSM}}_n}{\partial T}}& ={\displaystyle \frac{K_n}{T}}\left[G_n\right({\mathrm{RZSM}}_{n-\mathrm{1}}-{\mathrm{RZSM}}_n)\\ & +e^{-\frac{t_n-t_{n-\mathrm{1}}}{T}}{\displaystyle \frac{T}{K_{n-\mathrm{1}}}}{\displaystyle \frac{\partial {\mathrm{RZSM}}_{n-\mathrm{1}}}{\partial T}}],\end{array}\end{array}$$

with G_n defined as

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

σ(RZSM) and σ(T) denote the uncertainty in the RZSM estimates and the EF model parameter T, respectively. The equation is initialized as Δ₀=σ(SSM₀), $\partial {\mathrm{RZSM}}_{\mathrm{0}}/\partial T=\mathrm{0}$, and G₀=0. Uncertainties in the SSM input data are considered by the Δ term, which also takes into account the effect of possible prolonged input data gaps dependent on the T value. The Jacobian term $\partial \mathrm{RZSM}/\partial T$ assumes high values proportional to the latest SSM input variability on a timescale related to the T parameter. This is reflected in significant changes in the RZSM value associated with wetting or drying of the soil.

3.2.2 T-parameter uncertainty

De Santis and Biondi (2018) used an arbitrary value of σ(T) equal to 10 % of locally calibrated T_opt. This is in line with other studies on SM uncertainty propagation (Parinussa et al., 2011; Pathe et al., 2009), who used this uncertainty percentage for parameters without well-defined accuracy. In our study, we determine T_opt values based on a limited number of available in situ time series and apply these values to estimate RZSM globally. Consequently, σ(T) is likely to be greater due to a variety of environmental conditions not accounted for or underrepresented in the available in situ sample. We therefore propose the median absolute deviation (MAD) of T_opt (Sect. 3.1.2) as a more appropriate proxy for σ(T). In this case, the MAD is preferred over the variance because the sampling distribution of T_opt is both non-Gaussian and bounded (Leys et al., 2013).

3.2.3 EF model structural uncertainty

Recall that the standard law for the propagation of uncertainty (which is used in the baseline method) does not account for model structural uncertainty in the EF, which, due to the simplistic nature of the method and the limited surface–root-zone coupling, can account for a significant portion of the overall uncertainty budget.

We propose to estimate model structural uncertainty (σ(EF)) from in situ data using stations that operate sensors both at the surface and in the root zone. At these stations, we derive RZSM estimates from the SSM measurements using the EF method and then compare them to actual RZSM station measurements. For this analysis, the T value was optimized for each station and depth individually to minimize its influence on the estimation of σ(EF). This provides direct estimates for σ(EF) as

$$\begin{array}{}\text{(8)}& \mathit{\sigma }\left(\mathrm{EF}\right)=\mathrm{ubRMSD}({\mathrm{RZSM}}_{\mathrm{EF}},{\mathrm{RZSM}}_{\mathrm{ISMN}}),\end{array}$$

where ubRMSD denotes the unbiased root-mean-square difference. Note that “unbiased”, in this case, refers not only to a correction for bias in the mean (as is most commonly done) but also to a correction for bias in variance, which also constitutes an unintended systematic component in the RMSE (Gupta et al., 2009). Only sites with measurements from more than a single depth and at least one sensor within the surface layer (≤10 cm) were selected. Time series with negative correlation between EF-based RZSM estimates and in situ RZSM measurements were disregarded. As a result, a total of 1509 in situ sites were considered. Note that the EF model structural uncertainty computed at the point scale is assumed to be representative of the coarse scale as well.

Finally, the EF structural uncertainties obtained from Eq. (8) add to the propagated RZSM uncertainty budget (Eq. 4) as

$$\begin{array}{}\text{(9)}& \mathit{\sigma }\left({\mathrm{RZSM}}_n\right)=\sqrt{{\mathrm{\Delta }}_n^{\mathrm{2}}+\left({\displaystyle \frac{\partial {\mathrm{RZSM}}_n}{\partial T}}{)}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\right(T)+{\mathit{\sigma }}^{\mathrm{2}}(\mathrm{EF})}.\end{array}$$

4 Results and discussion

In this section, we first show results of the T-parameter point-scale optimization. Next, we compare the gridded RZSM product globally to E5L. We then discuss the estimates for EF model structural uncertainties. Finally, we compare our RZSM uncertainty estimates with those obtained with the baseline method.

4.1 T-parameter optimization

After filtering out unreliable data (see Sect. 2.2), 3901 ISMN time series from 67 different measuring depths between 0 and 200 cm were available for the T-optimization process. Figure 1 shows the distribution of T_opt values binned into our four chosen RZSM layers (0–10, 10–40, 40–100, and 100–200 cm). The median T_opt values for these layers were 6, 15, 48, and 70 d, increasing with soil depth as expected (Paulik et al., 2014; Wang et al., 2017). These median T_opt values were then used to compute RZSM globally.

Figure 1 T_opt values calibrated with 3901 in situ time series and binned according to RZSM layers 1–4. Median values (represented by orange lines) from each bin were used to compute a global RZSM product. Median absolute deviations (MAD(T_opt)) were used to estimate RZSM uncertainties.

A fivefold cross-validation was performed to verify the robustness of this approach. The variability in median T_opt values per soil layer increases with depth but remains negligibly small in all layers, with 6, 15–16, 47–50, and 67–72 for soil layers 1–4, respectively (Fig. 2a). Subsequently, the five median T_opt values derived from the training subsets were used to estimate RZSM for the different layers of the respective validation sets (Fig. 2c) and resulted in Pearson's r of 0.64–0.67, 0.64–0.65, 0.57–0.6, and 0.48–0.6 for soil layers 1–4, respectively. When evaluating each training set directly, correlations were 0.65–0.66, 0.65, 0.58–0.59, and 0.53–0.56 for soil layers 1–4, respectively (Fig. 2b).

Figure 2 Cross validation results showing the spread in T_opt values (a) and agreement of the training (b) and validation (c) sets with in situ data.

The little variability between the validation and training sets suggests that T_opt values are not over-fitted to ISMN site conditions and can be used robustly in other regions as well. Notably, the spread in median T_opt values increases with soil depth, while the correlation scores decrease. This indicates reduced reliability of the method in deeper soil layers, which is in line with the assumption that the coupling between the surface and root-zone SM decreases with depth. Note, however, that results for deeper layers are also affected by the smaller sample sizes at greater depths.

4.2 Global RZSM product quality assessment

A global SM dataset spanning the period 2002–2020 was computed using the EF method and T parameters optimized at point-scale with the approach described in Sect. 4.1. Figure 3a–e show correlation maps of each of the RZSM product layers as well as the input C3S SSM dataset with E5L; median correlation coefficients are summarized in Fig. 3f.

The spatial patterns observed in the C3S SSM data (Fig. 3a) are strikingly similar in RZSM layer 1 (Fig. 3b), with observed slight to moderate deterioration in performance over the high latitudes (>60^∘ N). This is not surprising given that both products differ only by a small degree of smoothing applied to RZSM layer 1 and are compared to the same E5L layer (0–7 cm). RZSM layers 2 and 3 (Fig. 3c–d), respectively, are compared to E5L layers 7–28 and 28–100 cm, respectively, and largely preserve good performance in regions where the input C3S SSM product also performs well, i.e. Europe (bar Scandinavia), the Caspian and Aral sea basins, the eastern United States, India, Southeast Asia, South America, sub-Saharan Africa, and Australia. At the same time deterioration of performance is observed at high latitudes and in arid environments such as the Sahara desert and the Arabian Peninsula, where the reduced strength of coupling between the surface and root-zone dynamics hinders the EF performance (Yang et al., 2022). Areas of good and poor performance visible in RZSM layers 1–3 are not apparent in RZSM layer 4 (Fig. 3e), where agreement with the reference E5L is very heterogeneous. A few regions where the good performance observed in shallower layers was preserved are India, Southeast Asia, and the eastern United States.

The median Pearson correlations between the RZSM product and the E5L reference layers (0–7, 7–28, 28–100, and 100–289 cm) were 0.54 (RZSM layer 1), 0.47 (RZSM layer 2), 0.41 (RZSM layer 3), and 0.28 (RZSM layer 4), respectively (Fig. 3f). The C3S SSM data included in this comparison show only marginally better performance against the E5L reference (0–7 cm) than the RZSM layer 1, with a median Pearson correlations of 0.55 and 0.54, respectively.

Figure 3 Spatial correlation maps of the C3S SSM and RZSM products with E5L SM (a–e).

The results are also consistent with the assumption of the EF model that SM dynamics decrease with depth and that T_opt ought to increase accordingly, as was also found by other studies (Wagner et al., 1999; Paulik et al., 2014; Wang et al., 2017; Beck et al., 2009; Mahmood and Hubbard, 2007). At the same time, the maximum correlation values decrease with depth, confirming the diminishing coupling between the surface and root-zone layers, as also found at the in situ station level (Fig. 2) and demonstrated by others (Paulik et al., 2014; Brocca et al., 2010; Sure and Dikshit, 2019).

4.3 EF model structural uncertainty

Figure 4 shows estimates for the model structural uncertainties (Sect. 3.2.3) obtained at all available in situ sites, binned into the four RZSM product layers. Their median values (represented by orange lines and annotated) were used as estimates for σ(EF). Note that in situ measurement errors were assumed to be negligible and thus did not influence ubRMSD estimates, which likely causes model structural uncertainties to be overestimated. Also, structural uncertainties are assumed to be constant in time.

Figure 4 The ubRMSD between propagated RZSM from in situ SSM using the EF model and measurements of RZSM at the same location and the same depth, calculated at 1509 different sites. The median ubRMSD value for each bin (represented by orange lines and annotated) represents σ(EF) for the respective T_opt.

As anticipated, an increase in σ(EF) corresponds to the growing distance between the surface and the root-zone measurements, demonstrating the decreasing coupling strength between both layers. Note that σ(EF) shows significant variability within RZSM layers, which is likely, at least to some degree, related to variations in local conditions. However, as with the T-parameter optimization, we estimate structural uncertainties based only on a limited number of in situ stations and therefore use the median to extrapolate globally.

4.4 RZSM uncertainty budget calculation

Figure 5 compares RZSM uncertainty estimates obtained from the baseline method (De Santis and Biondi, 2018) with those from the approach proposed here. Figure 5a shows a time series of RZSM uncertainties from the baseline method at an arbitrary location. Figure 5b shows the effect of changing σ(T) from 10 % of T_opt to the median absolute deviation of T_opt, which is an amplified temporal variability. Simultaneously, mean uncertainty values increase and become closer to the magnitudes of the input SSM dataset. Moreover, they no longer diminish with increasing T values (i.e. depth), as is the case in the baseline formulation. This is presumably more realistic since the progressive decoupling between the surface and deeper soil layers can be expected to cause uncertainties to increase rather than to decrease.

Figure 5c shows the impact of accounting for σ(EF) in the total uncertainty budget when using 10 % of T_opt as the T-parameter uncertainty (σ(T)). Considering this term substantially increases the magnitude of the propagated uncertainties and leads them to increase with depth (as does σ(EF)). However, the uncertainties' temporal variability is reduced substantially as the effect of σ(T) is overshadowed by that of σ(EF). Finally, Fig. 5d shows the combined effect of using the MAD of T_opt as its parameter noise σ(T) and accounting for model structural uncertainty σ(EF). Compared to the baseline (Fig. 5a), this yields an increased overall magnitude of the uncertainties; a more realistic increase in (temporal average) uncertainties with depth; and an amplified temporal variability in all layers during transitions between dry and wet conditions, which is also expected (see Fig. 6).

Figure 5 Evaluation of the impact of changes to the baseline method illustrated using an example 2020 time series from 9.875^∘ N, 1.625^∘ E. C3S SSM uncertainties were propagated with the baseline scheme in (a), while (b) and (c) show the individual impacts of increasing the noise of T from 10 % of T_opt to MAD(T_opt) and adding the term σ²(EF), respectively. Combined effects of both changes are shown in (d). The dashed grey line indicates the uncertainty level defined by GCOS (2022) as an accuracy goal for RZSM products.

4.5 Assessment of uncertainty estimates

Similar to De Santis and Biondi (2018), we assess the use of the proposed MAD estimates for σ(T) by computing Pearson's r and root-mean-square differences (RMSDs) with respect to in situ data before and after removing a fixed percentage of the data (5 %, 10 %, 15 %, and 20 %) with the highest uncertainty estimates. In the case of effective correspondence between high values of both the estimated RZSM uncertainties and the observed RZSM deviations from reference in situ measurements, it is expected that the skill metrics will improve due to the masking. This hypothesized correspondence holds well as long as the difference between in situ and satellite-based RZSM values is mainly due to the random errors in the latter. Note that this analysis is only sensitive to the impact of using different values for σ²(T) ($(T_{\mathrm{opt}}/\mathrm{10}{)}^{\mathrm{2}}$ versus MAD(T_opt)²) since the estimated structural uncertainty σ(EF) is constant in time and therefore cannot change the ranking of the total uncertainties.

Figure 6a and d indicate (in magenta shading) 20 % of RZSM layer 2 data with the highest uncertainties masked out in the experiment described above, based on uncertainties estimated with the baseline (b) and our method (d), respectively. Overall, despite the differences in magnitude and amplitude, both our and the baseline method assign the highest uncertainty values to timestamps corresponding to significant soil wetting or drying events. However, with the baseline method the average magnitude of SSM input uncertainty appears to have a greater influence on the calculated RZSM uncertainty estimates. This is most evident when comparing values before and after the inclusion of Metop-A ASCAT into the C3S product in January 2007 (indicated by the dashed vertical line in Fig. 6), which substantially improved data quality thereafter. Specifically, the mean C3S SSM uncertainty dropped from 0.029 to 0.018 m³ m⁻³. Such a clear shift is also visible in the uncertainty values propagated with the baseline method (from 0.008 m³ m⁻³ before to 0.004 m³ m⁻³ after the introduction of Metop-A ASCAT). This causes the baseline method to predict that the majority of the 20 % most uncertain SM values will occur in the pre-ASCAT period. In contrast, in our approach, average uncertainties remain stable (at 0.036 m³ m⁻³) over the entire time period. This suggests that the use of MAD(T_opt)² as an estimate for T-parameter uncertainty reduces the sensitivity to structural breaks, i.e. large variations between the uncertainties in the C3S SSM input sensors, and improves the method's capability to predict day-to-day uncertainty variations. Lastly, after the introduction of ASCAT, both schemes consistently assign higher uncertainties to timestamps characterized by large SM changes. Taken together, while the use of $T_{\mathrm{opt}}/\mathrm{10}$ as the T-parameter uncertainty seems to yield realistic estimates for uncertainty variations due to the use of different C3S SSM input sensors, using MAD(T_opt) as the T-parameter uncertainty seems to better predict day-to-day uncertainty variations in the RZSM estimates.

Figure 6 Differences in uncertainty variations in the baseline (a–b) and our proposed uncertainty estimation approach (c–d). Illustrated using the example of RZSM layer 2 at an arbitrary location in Benin (9.875^∘ N, 1.625^∘ E).

Figure 7 shows the results of the data removal experiment described above, summarized for all considered ISMN stations. To compare the performance with and without the effect of C3S structural breaks on the uncertainty values (see above), results are shown for both the full product period (2002–2020; Fig. 7a–d) and a sub-period without breaks, i.e. from the inclusion of SMAP data onward (1 April 2015–2020; Fig. 7e–h). In both cases, correlation coefficients obtained for the complete time series were compared to those obtained after removing 5 %, 10 %, 15 %, and 20 % of data with the highest associated uncertainties.

Figure 7 Correlations with in situ measurements (y axis) before and after removing a fixed percentage of data with the highest uncertainty (x axis) for the period 2002–2020 (a–d) and 2015–2020 (e–h). Uncertainties were calculated using either $\mathit{\sigma }\left(T\right)=T_{\mathrm{opt}}/\mathrm{10}$ (olive colour) or σ(T)=MAD(T_opt) (orchid colour).

In the case of the full product period (Fig. 7a–d), using $\mathit{\sigma }\left(T\right)=T_{\mathrm{opt}}/\mathrm{10}$ as the T-parameter uncertainty seems to consistently yield more realistic estimates of temporal uncertainty variations than using σ(T)=MAD(T_opt). This is true for all four soil layers. Masking out more uncertain data indicated by either method consistently improves agreement with in situ reference data in the first two product layers. This improvement increases the more data are masked out, as is expected. In the absence of such breaks (Fig. 7e–h) RZSM uncertainty variations seem to be better predicted when using σ(T)=MAD(T_opt) as the T-parameter uncertainty in almost all cases. Notably, in layers 3 and 4, data removal according to either method degraded the agreement with field measurements.

In summary, the propagation of C3S SSM input uncertainties yields accurate predictions of temporal uncertainty variations in RZSM estimates obtained with the EF method for the first two layers (0–10 and 10–40 cm). This is no longer the case for deeper layers (40–100 and 100-200 cm). Note, however, that the RZSM estimates in these layers themselves still exhibit reasonable skill when evaluated against E5L (see Fig. 3).

5 Summary and conclusions

In this study, we computed root-zone soil moisture (RZSM) globally in four depth layers (0–10, 10–40, 40–100, and 100-200 cm) from merged satellite surface soil moisture (SSM) retrievals of the Copernicus Climate Change Service (C3S) COMBINED product v202012 using the exponential filter (EF) method. The EF model parameter T has been optimized at point scale by maximizing the correlation with globally distributed in situ SM measurements from the International Soil Moisture Network (ISMN). The medians of the optimized T values at each layer have been used to compute the global product. A global product intercomparison with ERA5-Land (E5L) reanalysis SM data has shown a satisfactory level of agreement in all layers (global median correlations of the four above-mentioned product layers against E5L reference layers 0–7, 7–28, 28-100, and 100-289 cm were 0.54, 0.47, 0.41, and 0.28, respectively).

Uncertainties in the RZSM estimates obtained with the EF method were calculated using the law for the propagation of uncertainties. Uncertainties in the input SSM data were available in the C3S product and have been calculated by the data producers using triple collocation analysis (TCA). We tested the use of the median absolute deviation of optimized T parameters at the available ISMN locations (MAD(T_opt)) as a proxy for T-parameter noise. Results obtained using MAD(T_opt) in uncertainty propagation were compared with results obtained using 10 % of the optimized T parameter itself ($T_{\mathrm{opt}}/\mathrm{10}$), as done in earlier studies. While the use of $T_{\mathrm{opt}}/\mathrm{10}$ as the T-parameter uncertainty seems to yield realistic estimates for uncertainty variations due to the use of different C3S SSM input sensors, using MAD(T_opt) as the T-parameter uncertainty seems to better predict day-to-day uncertainty variations in the RZSM estimates.

Even though propagating SSM input and model parameter uncertainties yields credible predictions of temporal uncertainty variations, absolute uncertainty magnitudes appear unrealistically small (below 0.01 m³ m⁻³). This is because the propagation of uncertainty only accounts for uncertainties in the data and parameters input to the EF method, but not for limitations of the EF method itself (e.g. the progressive inability of the method to model deeper-layer RZSM due to vanishing surface–root-zone coupling). We proposed to estimate these EF model structural uncertainties as the unbiased root-mean-square differences between RZSM estimates for each of our four product depth layers obtained by applying the EF method to in situ SSM measurements and actual in situ RZSM measurements taken at the same location and depth. This was done at all available ISMN sites, and the median of these estimates was used as a global proxy for EF structural uncertainty for each of the four product depth layers, respectively. Combined, propagated SSM input and model parameter uncertainties and EF structural uncertainties were considered to yield realistic estimates of the total RZSM product uncertainty budget in all layers (global mean uncertainties in the four product layers are 0.031, 0.035, 0.04, and 0.04 m³ m⁻³). Note, however, that a quantitative validation of uncertainty magnitudes is still pending due to the lack of reliable uncertainty reference data on a global scale and for different RZSM depth layers.

The EF parameter uncertainty was estimated on a global scale and can be expected to differ for smaller scales, especially where the variability in environmental conditions is lower. Similarly, estimates of the EF model structural uncertainty are likely to differ on local to regional scales. Also, the structural uncertainty in the EF, here assumed to be constant in time, could in fact vary on a sub-seasonal scale given the phenomena that regulate the process of water transfer in the soil. Moreover, random errors in the in situ measurements were assumed to be negligible and were not accounted for in estimating the structural uncertainty in the model. Nonetheless, it is plausible that the EF structural uncertainty is much greater than the random uncertainty in the in situ sensors. Estimates of the random uncertainty in the in situ sensors could allow for a more accurate estimation of the EF structural uncertainty in the future.

Further insights could also be gained by evaluating the behaviour of the proposed method in propagating uncertainties in different SSM input data, e.g. single-sensor products without structural breaks and non-static input SSM uncertainties obtained by means other than TCA. Nonetheless, this study is an important step towards understanding and describing the uncertainties in EF-based RZSM products.

Appendix A: ISMN references

Mougin et al. (2009)Cappelaere et al. (2009)de Rosnay et al. (2009)Lebel et al. (2009)Galle et al. (2015)Cook (2016a)Cook (2016b)Cook (2018)Musial et al. (2016)Van Cleve et al. (2015)Brocca et al. (2011)Brocca et al. (2011)Zreda et al. (2008)Zreda et al. (2012)Yang et al. (2013)Tagesson et al. (2015)Ikonen et al. (2016)Ikonen et al. (2018)Al-Yaari et al. (2018)Wigneron et al. (2018)Xaver et al. (2020)Zappa et al. (2019)Zappa et al. (2020)Kang et al. (2014)Jin et al. (2014)Blöschl et al. (2016)Vreugdenhil et al. (2013)Jensen and Refsgaard (2018)Bircher et al. (2012)Flammini et al. (2018a)Flammini et al. (2018b)Morbidelli et al. (2011)Morbidelli et al. (2014)Morbidelli et al. (2017)Hollinger and Isard (1994)Biddoccu et al. (2016)Raffelli et al. (2017)Capello et al. (2019)Alday et al. (2020)Osenga et al. (2019)Osenga et al. (2021)Mattar et al. (2014)Mattar et al. (2016)Su et al. (2011)Dente et al. (2012)Beyrich and Adam (2007)Kang et al. (2019)Su et al. (2011)Dente et al. (2012)Su et al. (2011)Dente et al. (2012)Young et al. (2008)Smith et al. (2012)Larson et al. (2008)Hajdu et al. (2019)González-Zamora et al. (2019)Canisius (2011)L'Heureux (2011)Ojo et al. (2015)Rüdiger et al. (2007)Schaefer et al. (2007)Nguyen et al. (2017)Zhao et al. (2020)Zheng et al. (2022)Calvet et al. (2007)Albergel et al. (2008)Calvet et al. (2016)Leavesley et al. (2008)Moghaddam et al. (2010)Moghaddam et al. (2016)Shuman et al. (2010)Marczewski et al. (2010)Zacharias et al. (2011)Bogena et al. (2012)Bogena et al. (2018)Bogena (2016)Loew et al. (2009)Schlenz et al. (2012)Brocca et al. (2008)Brocca et al. (2009)Brocca et al. (2011)Bell et al. (2013)Jackson et al. (2010)Kirchengast et al. (2014)Fuchsberger et al. (2021)Petropoulos and McCalmont (2017)

Table A1 ISMN networks used in this study. A list of all available ISMN networks can be found at https://ismn.earth/en/networks/, last access: 28 August 2023.

Figure A1 Location map of the ISMN in situ stations used in this study and listed in Table A1.