<figure><img src="image_1.png" /></figure>

HAL
open science

# Combining Datasets of Satellite-Retrieved Products.
## Part I: Methodology and Water Budget Closure

Filipe Aires

► To cite this version:

Filipe Aires. Combining Datasets of Satellite-Retrieved Products. Part I: Methodology and Water Budget Closure. Journal of Hydrometeorology, 2014, 15 (4), pp.1677-1691. 10.1175/JHM-D-13-0148.1. hal-02559427

HAL Id: hal-02559427

https://hal.science/hal-02559427

Submitted on 12 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Combining Datasets of Satellite-Retrieved Products. Part I:
Methodology and Water Budget Closure

FILIPE AIRES

Estellus, France, and Laboratoire de l'Etude du Rayonnement et de la Matière en Astrophysique, CNRS, Observatoire de Paris, Paris,
France, and Earth Institute and Department of Earth and Environmental Engineering, Columbia University, New York, New York

(Manuscript received 29 August 2013, in final form 5 March 2014)

## ABSTRACT

This study addresses in general terms the problem of the optimal combination of multiple observation datasets. Only satellite-retrieved geophysical parameter datasets are considered here (not the raw satellite observations). This study focuses on the terrestrial water cycle and presents methodologies to obtain a coherent dataset of four water cycle key components: precipitation, evapotranspiration, runoff, and terrestrial water storage. Various innovative “integration” methodologies are introduced: simple weighting (SW), constrained linear (CL), optimal interpolation (OI), and neural networks (NN). The term “integration” will be used here, not “assimilation,” as no model will be included in the data fusion process. A simple post-processing filtering (PF) step can be used to impose the water cycle budget closure after the integration method. It is shown that this constraint actually improves the estimation of the water cycle components. The integration techniques are tested using real observation data over the Mississippi and Niger basins from satellite and in situ measurements. A Monte Carlo experiment with a synthetic uncertainty perturbation model is used to measure the ability of the SW, OI, and NN, with or without the PF step, to retrieve the four water cycle components. Once the PF closure constraint is added, the methodologies have equivalent accuracies. The need for these types of methodologies should increase in the future since multiple observation datasets are now available and the climate community needs to combine them into a unique, optimal, and coherent dataset of multiple parameters. A companion paper will test these methodologies on satellite observation datasets at the basin and global scales.

## 1. Introduction

The number of satellites observing Earth's surface and atmosphere has been constantly increasing over the last three decades (Fernandez-Prieto et al. 2012). In addition, the exponential increase of computational power has allowed more and more groups (i.e., operational centers, space agencies, universities, public laboratories, or private companies) to develop their own satellite retrievals from these observations. As a consequence, multiple observation datasets exist for each geophysical variable, and the scientific community is now confronted with the problem of choosing a single dataset for each variable. However, each dataset source has its own advantages and drawbacks, so this choice is

difficult. Another difficulty is that the available datasets for different geophysical variables are rarely coherent with each other because each one generally uses different assumptions, auxiliary information, or models. These incoherencies in multivariable datasets complicate their use in studies focusing on variable dependencies, such as process-oriented analysis or sophisticated validation metrics (Aires and Rossow 2003; Foley et al. 2013). In this context, international programs such as the Global Energy and Water Cycle Experiment (GEWEX; www.gewex.org/GDAP.html) are focusing on the integration of multiple datasets with a particular emphasis on coherency constraints.

As a consequence, there is a strong need for methodologies able to optimally produce, in a coherent way, multivariable datasets from the fusion of various observation datasets. These types of methods should be able to assimilate very different sources of information, such as 1) satellite observations, 2) in situ measurements, and 3) physical constraints, when they are available. Variational

<footnote>Corresponding author address: F. Aires, Estellus/LERMA, Observatoire de Paris, 61 avenue de l'Observatoire, 75014 Paris, France.
E-mail: filipe.aires@estellus.fr</footnote>
assimilation (Rodgers 2000; Kalnay et al. 1996) could be a good candidate, but in order to obtain a purely observational dataset, it is important not to use any model (hydrological, surface, or climate). This is important because one of the main objectives for building such a multivariable dataset is to obtain a reference that is as independent as possible from any model, so that it can be used to validate and calibrate models.

In this paper, various methodologies will be presented and tested by studying the terrestrial water cycle budget. The continental water budget is a key diagnostic of the effects of climate change and variability. It is also the most important part of the water cycle for impact studies on the populations, and its change remains uncertain. The closure of the terrestrial water cycle has been studied using the following water budget equation:

$$ \frac{\partial S}{\partial t} = P - E - R, \qquad (1) $$

where S represents the terrestrial water storage, R is the global runoff, P is the precipitation, and E is the evaporation (evapotranspiration).

Experiments have been conducted to close this budget (Sheffield et al. 2009; Sahoo et al. 2011; Pan et al. 2012). Most studies have targeted the runoff (freshwater discharge) as a diagnostic. It is at the interface between land and ocean, and it can be measured by in situ observations (even if there are now attempts to measure it with the use of satellite observations). In Sheffield et al. (2009), precipitation, evapotranspiration, and water storage are quantified using a collection of satellite measurements while the streamflow is used as the closure target and estimated with gauge measurements. Large uncertainties in the individual components make the closure difficult. For 10 basins over 3 years, Sahoo et al. (2011) merge different products for precipitation and evapotranspiration, using weighted values based on their errors, to mitigate errors in the individual satellite products. A constrained ensemble Kalman filter is then used to close the water budget and provide a constrained estimate of the water budget. Pan et al. (2012) adopt a similar strategy to estimate the water budget over 32 basins for the period of 1984–2006. However, they used in situ observations, land surface model simulations, and global reanalyses together with the remote sensing products. In Azarderakhsh et al. (2011), the Amazon basin is studied using water storage, precipitation, and evaporation satellite observations. The water discharge is estimated using the water balance equation and results are compared to gauge-based measurements. This study stresses the need for subbasin studies and the addition of temporal lags in the analysis.

From the previous experiments and literature, it is clear that each individual dataset is not accurate enough to directly close the water budget. In this paper, the objective is then to estimate new optimized states $X_a = (P, E, Q, \Delta S)$ from various estimates of the different components, plus additional constraints, in particular, the water budget closure. Pan and Wood (2006), Sahoo et al. (2011), and Pan et al. (2012) have performed a similar task, but they employed a Kalman filtering solution, and a land surface model was used. In this study, the emphasis is to obtain purely observation-driven estimates, limiting the intervention of any model. The qualities and limitations of various approaches will be investigated. The methods should utilize as many sources of information as possible for each water cycle component (e.g., various satellite precipitation estimates), but no model. They should combine these estimates by using a priori information on their respective uncertainties; the combination of sources is intended to reduce these uncertainties. They should use the water budget closure as independent a priori information to help characterize uncertainties and to obtain a new optimized dataset of terrestrial water components. Furthermore, the methodology should be flexible enough to ensure that as much a priori information as possible about these water components can be used. The weighting and closure constraints should be performed simultaneously to obtain an optimal solution. The approach should be flexible enough to ensure that it can evolve with time, when more information on the water cycle can be integrated. A thorough analysis of the results should be performed over space and time (by analyzing the space-time variability of the weighting) to ensure the best quality for a long-term and globally integrated dataset that can be used for climate studies.

Section 2 presents various integration techniques. The uncertainties of these methods will be assessed using Monte Carlo experiments, and section 3 describes the approach used to perform this task. The results are presented and analyzed in section 4. Section 5 provides the main conclusions of this study and presents various perspectives for this work.

## 2. Integration methods

### a. General context

The goal of this study is to introduce different methods able to combine various datasets of multiple variables, in such a way that a priori constraints among the state variables are imposed to the state estimate.

The major developments for the introduction of state constraints have been accomplished in the framework of
Kalman filtering. A first approach to use such state constraint consists of projecting the state variables in a space where the constraint is satisfied. However, this approach has some limitations. For example, the interpretability of the results is more difficult since both the state variables and the physical equations are modified (Simon and Chia 2002). Another solution uses the particular constraint as a “perfect measurement,” meaning that there is no uncertainty on this constraint. The drawback of this method is that the introduction of such “hard” or “strong” constraints into the system increases the problem dimension and can introduce numerical instabilities (Hayward 1998; Porrill 1988; De Geeter et al. 1997; Kalman 1960). Simon and Chia (2002) proposed to postprocess the unconstrained Kalman solution to impose a linear state constraint: at each time step, the unconstrained solution is projected into a constrained state space. This avoids the numerical instabilities associated with the perfect measurement approach. Furthermore, it does not increase the dimension of the problem since the state variables are kept in the same space (this facilitates the interpretation of results). In addition, it is shown that this projected solution is better than the unconstrained one since the elements of the error covariance matrix are smaller in magnitude. This postprocessing can be defined using maximum a posteriori, least squares, or linear projection approaches.<fnref n="1" /> An even more advanced solution is proposed by Ko and Bitmead (2007) that projects the full dynamical system into a constrained space. This system-projection approach is shown to be even better than the postprocessing solution.

However, these developments have been made in the framework of Kalman filtering where a state variable dynamical system is considered. This often means that a model is used so that the state variables at time $t$ are dependent on the state variables at earlier times (i.e., the dynamical system). As previously mentioned, the goal of this paper is to produce state estimates that are independent from any model. Therefore, the Kalman filtering cannot be considered here; only direct state estimation will be proposed. However, the technical developments made for the Kalman filtering will be useful for the integration tools proposed in the following.

## b. Notations

Let $\mathbf{X}^T = (P, E, Q, \Delta S)$ be the state variable quantifying the four water cycle components of interest in this paper (a superscript T is the transpose symbol). The goal of this study is to estimate $\mathbf{X}$. The closure of the water

cycle is obtained when $\mathbf{X}^T \cdot \mathbf{G} = 0$, where $\mathbf{G}^T = (1, -1, -1, -1)$, from Eq. (1).

Multiple observations can be obtained for each component:

* $(P_1, P_2, \dots, P_p)$, the $p$ precipitation estimates;

* $(E_1, E_2, \dots, E_q)$, the $q$ sources of information for evapotranspiration;

* $(Q_1, Q_2, \dots, Q_r)$, the $r$ runoff estimates; and

* $(\Delta S_1, \Delta S_2, \dots, \Delta S_s)$, the $s$ sources of information for the groundwater storage change.

Except for the runoff $Q$ that comes from in situ measurements, in the following, “observation” will be referring to a satellite observation. Note also that, in the following, observations will refer to water cycle components (for $P$, $E$, $Q$, or $\Delta S$) and not to actual raw satellite observations such as radiances measured by satellite instruments. The aim of this study is to optimally combine observation products, not to develop a remote sensing algorithm.

Let

$$ \mathbf{Y}_e^T = (P_1, \dots, P_p, E_1, \dots, E_q, Q_1, \dots, Q_r, \Delta S_1, \dots, \Delta S_s) \quad (2) $$

be the vector of dimension $n = p + q + r + s$ gathering all the initial water cycle observations. The observing system can be represented using<fnref n="2" />

$$ \mathbf{Y}_e = \mathbf{A} \cdot \mathbf{X} + \varepsilon; \quad (3) $$

no bias error is considered here. The observation operator $\mathbf{A}$ has the following structure:

$$ \mathbf{A} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 1 \end{pmatrix}. \quad (4) $$

<footnote n="1">¹ As often seen, the various solutions are equivalent when the same assumptions are used (i.e., linearity of the operators and Gaussian character of the considered random variables).</footnote>

<footnote n="2">² In this case, observations are equal to the actual water cycle component plus an additive noise. More complex observing systems could be considered here and, in this case, the linear function $\mathbf{A}$ would be replaced by a more complex relation.</footnote>
Off-diagonal terms could be introduced if this information was available.

In this matrix, row i includes a “1” in the first column if observation $\mathbf{Y}_i$ measures $P$ (or $E$, $Q$, or $\Delta S$ for the second, third, and fourth columns, respectively). The $n$-vector $\boldsymbol{\varepsilon}$ represents the uncertainties of each observation in $\mathbf{Y}_\varepsilon$. The observation uncertainty is characterized by the covariance matrix $\mathbf{R} = \text{cov}(\boldsymbol{\varepsilon})$, often chosen to be diagonal:

$$ \mathbf{R} = \begin{pmatrix} \sigma_1^2 & 0 & 0 \\ 0 & \ddots & 0 \\ 0 & 0 & \sigma_n^2 \end{pmatrix}. \qquad (5) $$

In this study, the goal is to estimate a relationship (i.e., filter) that finds an estimate $\mathbf{X}_a$ (a stands for analysis) of $\mathbf{X}$ based on the observations $\mathbf{Y}_\varepsilon$. Therefore, an inverse problem needs to be solved where Eq. (3) is inverted from observation $\mathbf{Y}_\varepsilon$ to $\mathbf{X}_a$ (Tarantola 1987; Rodgers 2000).

A simple inversion of matrix $\mathbf{A}$ is not good enough. Since $\mathbf{A}$ is not square, a pseudoinverse or a generalized inversion of matrix $\mathbf{A}$ could be a candidate, but these simple algebra techniques suffer from numerical instabilities. In addition, other sources of information are available, so more sophisticated methods should be considered to exploit them. For example, the information on the observation noise (i.e., matrix $\mathbf{R}$) is very valuable information that allows for weighing the observations $\mathbf{Y}_\varepsilon$ based on their reliability. Furthermore, in addition to the observations $\mathbf{Y}_\varepsilon$, a priori information $\mathbf{X}_b$ on $\mathbf{X}$ can also be used (b stands for “background,” the terminology used in numerical weather prediction centers). An important question concerns of course the origin of this first guess (FG) information since it is not always easy to obtain such a priori information, especially when models are excluded. The covariance matrix $\mathbf{B}$ of the FG errors needs to be specified based on the source of information for $\mathbf{X}_b$.

### c. Simple weighting

Let us first consider the simple case where two observations $P_1$ and $P_2$, with Gaussian uncertainties $\varepsilon_1 \sim N(0, \sigma_1)$ and $\varepsilon_2 \sim N(0, \sigma_2)$, are available for the estimation of precipitation $P$. In this case, the optimal linear estimator (Rodgers 2000) is given by

$$ P_a^{\text{sw}} = \frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2} P_1 + \frac{\sigma_1^2}{\sigma_1^2 + \sigma_2^2} P_2. \qquad (6) $$

This type of equation is valid when there is no bias error in the sources of information $P_1$ and $P_2$. Note that to be an optimal solution, the errors in $P_1$ and $P_2$ must be independent from each other. This type of estimator is

more frequently used when merging raw satellite observations because the bias correction is easier in this case (Thome 2004; Aires et al. 2010). However, the variables that are combined here are geophysical products originating from a remote sensing algorithm and, most of the time biases do exist in this case, unless a calibration procedure has been used to reduce them. When information on the bias is available (e.g., by comparison to in situ measurements), this information can be used to suppress the bias. However, in many cases, bias information is difficult to obtain, and a pragmatic strategy is to estimate the mean state for each variable as the average of each of its observations: $\bar{P} = \text{mean}(P_1, P_2)$ (the overbar refers to the mean). Then all the sources of information $P_i$ on this variable are bias corrected toward this average: $\bar{P}_i = P_i - \bar{P}_i - \bar{P}$, for $i = 1, 2$.

It is possible to extend the result of Eq. (6) to the case where more than two observations are available for a variable and where more than one variable is retrieved. Let us find the optimal linear estimator, such that

$$ \mathbf{X}_a^{\text{sw}} = \mathbf{K} \cdot \mathbf{Y}_\varepsilon, \qquad (7) $$

where $\mathbf{K}$ is a $4 \times n$ matrix. To be optimal, $\mathbf{K}$ needs to take into account the observation uncertainties of the observation $\mathbf{Y}_\varepsilon$ described by $\mathbf{R}$.

Let us consider the first line of $\mathbf{K}$ used to retrieve the first water cycle component in $\mathbf{X}$, that is, $P$. As stated earlier, there are $p$ available observations: ($P_1, P_2, \dots, P_p$). The first line of $\mathbf{K}$ is composed of zeroes, except for the first $p$ columns that are given by

$$ K_{1j} = \frac{1}{p-1} \left( \sum_{\substack{1 \le k \le p \\ k \ne j}} \sigma_k^2 \right) / \left( \sum_{1 \le k \le p} \sigma_k^2 \right) \cdot P_k, \quad (8) $$

for $j = 1, \dots, p$. This expression extends easily to the second (for $E$), third (for $Q$), and fourth (for $\Delta S$) rows. Again, this expression is the optimal solution only if all the errors are independent from each other. If a covariance matrix describing the dependencies of the input datasets errors was available, an eigendecomposition of this covariance matrix would create a new input dataset with uncorrelated errors. Unfortunately, this type of covariance matrix for the input dataset errors is generally not available.

### d. A constrained linear technique

This approach is based on a simple linear model for the weighting of the satellite observations:
$$ \mathbf{X}_{a}^{\mathrm{CL}} = \mathbf{K} \cdot \mathbf{Y}_{\mathrm{g}}, \qquad (9) $$

with additional constraints on gain matrix **K** and on solution **X**<sub>a</sub><sup>CL</sup>.

We note that **K** is the (4, n = p + q + r + s) matrix that linearly links the new water cycle components to the satellite estimates:

$$ \mathbf{K}^{\mathrm{T}} = \begin{pmatrix} k_{1,1} & 0 & 0 & 0 \\ \vdots & 0 & 0 & 0 \\ k_{1,p} & 0 & 0 & 0 \\ 0 & k_{2,p+1} & 0 & 0 \\ 0 & \vdots & 0 & 0 \\ 0 & k_{2,p+q} & 0 & 0 \\ 0 & 0 & k_{2,p+q+1} & 0 \\ 0 & 0 & \vdots & 0 \\ 0 & 0 & k_{2,p+q+r} & 0 \\ 0 & 0 & 0 & k_{4,p+q+r} \\ 0 & 0 & 0 & \vdots \\ 0 & 0 & 0 & k_{4,p+q+r+s} \end{pmatrix} \quad (10) $$

Note that the simple weighting (SW) filter follows such structural constraint.

The **K** coefficients in Eq. (10) are constrained by

$$ L_P = \sum_{i=1}^{p} k_{1,i} = 1, \quad L_E = \sum_{i=p+1}^{p+q} k_{2,i} = 1, $$

$$ L_{\Delta Q} = \sum_{i=p+q+1}^{p+q+r} k_{3,i} = 1, \quad \text{and} \quad L_{\Delta S} = \sum_{i=p+q+r+1}^{p+q+r+s} k_{4,i} = 1 \qquad (11) $$

(i.e., each line L in **K** adds up to 1). This ensures that estimated water components in **X** are a weighted sum of the observations for that particular component. However, if the estimation of **K** is a well-posed problem, it might be possible to use simpler constraints, for example, the coordinates in **K** can be constrained to be 0 ≤ k<sub>ij</sub> ≤ 1, or their amplitude can be limited. This later and simple procedure is a Tikhonov regularization often used in regression problems. For instance, in the neural network (NN) theory, it is called “weight decay” (Bishop 1996).

Note that if only constraints of Eqs. (11) were applied to the filter of Eq. (10), the constrained linear (CL) solution would be very close to the SW one. The SW exploits the uncertainty characteristics to obtain the good coefficients in filter **K**. The specificity of the CL

approach is that, in addition to the structural constraint in **K**, a water closure constraint is also added. This closure can be written as a constraint on the solution **X**:

$$ \mathbf{C}_{\mathrm{closure}} = \sum_{i=1}^{M} (\mathbf{G} \cdot \mathbf{W} \cdot \mathbf{Y}_{\mathrm{g}}^{i})^{\mathrm{T}} \nu^{-1} (\mathbf{G} \cdot \mathbf{W} \cdot \mathbf{Y}_{\mathrm{g}}^{i}), \quad (12) $$

where *M* is the number of samples (i.e., observations **Y**<sup>i</sup>) used in this criterion. The term *v* is a scalar; it represents the variance of closure errors. This number can be defined a priori from previous experiments or as an empirical regularization term. In this case, the smallest value for *v* that still ensures water budget closure in the estimate **X**<sub>a</sub> is chosen.

To find the optimal solution **K**, an optimization algorithm needs to be used to minimize **C**<sub>closure</sub> under constraints of Eqs. (11):

$\hat{\mathbf{K}} = \min_{\mathbf{K}}(\mathbf{C}_{\mathrm{closure}})$, under:

$$ L_P = 1, L_E = 1, L_Q = 1, \quad \text{and} \quad L_{\Delta S} = 1. \quad (13) $$

This optimization needs to be performed simultaneously on all the *N* available samples. The samples can be pooled into a general dataset, and the optimization be used to obtain a general filtering matrix **K**. As previously mentioned, matrix **K** has a similar structure to the SW solution. As a consequence, SW can be a good first guess solution for CL, and this can speed up the optimization process and the quality of the final results.

It should be noted that in this CL solution, **R** and **B** are not used, though they are also very valuable constraints. It is an intermediate solution between the SW and the following optimal interpolation (OI). This method will not be tested in the following experiments because it is more difficult to implement and the application that will be presented in this paper does not require such complexity. However, it is a very interesting approach, especially for problems where complex physical constraints need to be imposed on the solution.

### e. The optimal interpolation method

The goal of this approach is to simultaneously exploit all the available a priori information and to find the best compromise for **X**<sub>a</sub>. For this purpose, a cost function (i.e., quality criterion) is first defined as

$$ \begin{aligned} 2J(\mathbf{X}) = & (\mathbf{X} - \mathbf{X}_b)^{\mathrm{T}} \mathbf{B}^{-1} (\mathbf{X} - \mathbf{X}_b) \\ & + (\mathbf{Y}_g - \mathbf{A} \cdot \mathbf{X})^{\mathrm{T}} \mathbf{R}^{-1} (\mathbf{Y}_g - \mathbf{A} \cdot \mathbf{X}) \\ & + (\mathbf{G} \cdot \mathbf{X})^{\mathrm{T}} \mathbf{C}^{-1} (\mathbf{G} \cdot \mathbf{X}). \end{aligned} \quad (14) $$
The first right-hand term is related to the information provided by an FG variable $X_b$. Since the focus of this paper is to obtain a model-independent estimate of the water cycle components, this FG needs to originate from observations such as available climatology, simple averaging of the observations for each component, or the result of an SW filtering. The second term is related to the observations $Y_e$. The third term represents the water budget closure constraint. This OI filter can be used without the FG information $X_b$, but, of course, this would be detrimental to the quality of the solution.

The optimal solution $X_a^{OI}$ is chosen to minimize the quadratic cost function $J(X)$. This is obtained when equating to zero the gradient $\nabla J(X) = 0$. By posing $Y_e - A \cdot X = Y_e - A \cdot X_b - A \cdot (X - X_b)$, and neglecting the constant terms, we search for the gradient of

$$ (X - X_b)^T (B^{-1} + A^T R^{-1} A) (X - X_b) - 2 (X - X_b)^T A^T R^{-1} (Y_e - A X_b). \quad (15) $$

This gradient is given by

$$ \nabla J(\mathbf{X}) = (\mathbf{X} - \mathbf{X}_b)^T (\mathbf{B}^{-1} + \mathbf{A}^T \mathbf{R}^{-1} \mathbf{A}) + \mathbf{G}^T \mathbf{C}^{-1} \mathbf{G}) \cdot (\mathbf{X} - \mathbf{X}_b) - \mathbf{A}^T \mathbf{R}^{-1} (\mathbf{Y}_e - \mathbf{A} \mathbf{X}_b) + \mathbf{G}^T \mathbf{C}^{-1} \mathbf{G} \mathbf{X}_b. \quad (16) $$

Therefore, the solution follows:

$$ \mathbf{X}_a^{OI} = \mathbf{X}_b + \mathbf{K} [\mathbf{A}^T \mathbf{R}^{-1} (\mathbf{Y}_e - \mathbf{A} \mathbf{X}_b) - \mathbf{G}^T \mathbf{C}^{-1} \mathbf{G} \mathbf{X}_b], \quad (17) $$

where

$$ \mathbf{K} = (\mathbf{B}^{-1} + \mathbf{A}^T \mathbf{R}^{-1} \mathbf{A} + \mathbf{G}^T \mathbf{C}^{-1} \mathbf{G})^{-1}. \quad (18) $$

It can be seen that in this approach, the linear matrix $K$ is built using both the $R$ and $B$ matrices, and this is a valuable feature. No constraints of the form $L_p, L_e, L_Q$, and $L_{\Delta S}$ or no structural constraint such as in the filter $K$ of Eq. (9) are included in this solution, so there is a flexibility in the shape of the filter. Note that this OI filter does not change with the input data, except if the a priori information error statistics in matrix $B$ or the observations errors in $R$ change from one input datum to another.

In the following experiments, three tests will be conducted. The OI1 and OI2 solutions refer to the OI estimate with and without the closure term in Eq. (14) (i.e., last term in the equation with $G$). The solution OI3 will be the OI1 solution plus the postprocessing filtering (PF) step that will be introduced in the following to impose closure in a postprocessing step.

### f. A nonlinear neural network weighting technique

In this approach, the estimate $X_a^{NN}$ of the water components is provided by a nonlinear NN model $g$:

$$ \mathbf{X}_a^{NN} = g(\mathbf{Y}_e). \quad (19) $$

This is an extension of the SW, CL, or OI solutions that use a linear filtering instead of the nonlinear model $g$. This nonlinear character of $g$ means that the weighting of the information in $Y_e$ is state dependent; this is generally a useful feature.

The details of the NN methodology will not be discussed here. The interested reader is referred to, for example, Bishop (1996) and Aires et al. (2004b). The NN model $g$ is a parametric model, and a training procedure needs to be used to estimate these parameters. The training of $g$ first requires an optimization procedure. The usual back-propagation algorithm (Rumelhart et al. 1986) will be used here. A training dataset $\mathcal{D}$ is required too: it gathers a set of input-output couples $(\mathbf{Y}_e, \mathbf{X})$. In this study, these couples will be built with an ensemble of water cycle components $X$ that follow a closure constraint from Pan et al. (2012) (see section 3a). Synthetic observations $Y_e$ are generated using a specified noise model for $e$. The goal of the NN model $g$ is therefore to retrieve the target water cycle components $X$ from the perturbed observations $Y_e$. The noise model for the perturbations in $Y_e$ is not provided to the NN model; only their implicit characteristics in dataset $\mathcal{D}$ are available to the NN.

The goal of this study is to retrieve, from the observations $Y_e$, an estimation $X_a$ such that its components close the water cycle budget. How can this closure constraint be imposed on the NN outputs? The first solution is easy in practice; it consists of using the learning dataset $\mathcal{D}$ outputs $X$ that respect this closure. In general, a constraint in the learning dataset is reproduced by the trained NN, but this is not a strong explicit constraint on the model $g$. This is the solution tested in this paper since the dataset described in section 3a respects this budget closure. Another solution would be to introduce the closure constraint (previously defined as $C_{closure}$) in the loss function that is minimized during the NN training. Such complex constraints have been used, for example, in Aires et al. (1999), but their implementation is complex. In particular, the back-propagation algorithm would need to be redefined in this case. This approach will not be tested here.

In this NN technique, the weights for each water cycle component in $Y_e$ can vary with the state of the system. However, there is no constraint on the summation of the weights, as in the $L_p, L_e, L_Q$, and $L_{\Delta S}$ terms in
section 2d. This can potentially be dangerous because compensation phenomena could occur (Aires et al. 2004a): an error in the weight of one observation in $Y_e$ is compensated by an opposite error in the weight of another observation. This can lead to very high positive and negative weight values in the NN filter. This is a classic symptom of ill-conditioned problems. In the case of such problems, regularization solutions tend to reduce the potential for such high-value weights, for example, the traditional NN weight decay (Bishop 1996).

As previously mentioned, the learning of the NN requires a dataset that includes the four water cycle component targets that respect the closure constraint. As a consequence, the NN learning inference strategy cannot easily be used on real conditions. However, in the experiments of section 4, this technique will be used to test the state dependency of the obtained filters.

### g. The postprocessing filtering

This interesting approach has been presented, for example, in Pan and Wood (2006). It is not an integration method, but rather a postprocessing method to impose the closure constraint on a previously obtained solution $X_b$. Thus, an estimate of the state variables $X_b = K \cdot Y_e$ is first obtained using the observations $Y_e$ and some filter $K$. In Pan and Wood (2006), the FG uses also some constraints from the variable infiltration capacity (VIC) macroscale hydrologic model, but other FG solutions (from the SW, CL, OI, or NN estimates) can be used instead. This will be tested in section 4.

A PF process is then applied on the FG solution $X_b$ in order to introduce the closure constraint:

$$X_a^{\text{PF}} = X_b + K(0 - G X_b), \quad (20)$$

where $K = BG^T(GBG^T)^{-1}$. This can be rewritten as

$$X_a^{\text{PF}} = K_{\text{PF}} \cdot X_b, \quad (21)$$

where $K_{\text{PF}} = [\text{Id} - BG^T(GBT^T)^{-1}\text{G}]$, where $\text{Id}$ is the identity matrix. This expression shows that this post-processing is a filtering of $X_b$ based on $G$ with the goal of obtaining the closure. Note again that this filtering is constant and identical for any situation unless there is enough information to have an a priori error matrix $B$ dependent on the situation.

In this paper, this approach is used, for example, with $X_b = X_a^{\text{SW}}$, that is, the SW solution of section 2c. Note that this postprocessing step does not use the observations $Y_e$ (used for the FG $X_b$ only), so the observations are not used twice to estimate $X$, which would have been a problem (Rodgers 2000).

## 3. Dataset and Monte Carlo scheme

### a. Hydrological dataset

The numerical experiments of this paper use the terrestrial water cycle dataset produced in Pan et al. (2012). This dataset includes the $P, E, Q$, and $\Delta S$ components at the monthly time scale for 23 years (1984–2006) and over 32 major basins. Therefore, the dataset $D$ is a matrix of dimension $276 \times 4$ for each basin. In this dataset, the terrestrial water budget closure is respected since this was the goal of its construction. The method used to impose this closure was based on the PF method, with an FG provided by a combination of observations and modeling constraints. Figure 1 represents the four water components in $D$ for the Mississippi basin. The budget is not represented in this figure, but, by construction, it is equal to zero for all time steps. It can be noted that large interannual variations are present, but no easy long-term trend can be identified here.

The tests conducted in this paper are performed over two basins only: Mississippi and Niger. These basins are major hydrological regions. The mean and standard deviation (STD) of the four water components for these two basins are presented in Table 1. The average values are close for both basins but the variability can differ: evapotranspiration has a higher STD for Mississippi, but precipitation is the more variable component for Niger. In both cases, the runoff $Q$ has the lowest variability.

From this dataset $D$, synthetic observations are built for testing the ability of the integration methods of section 2 to retrieve $X$. The perturbed observations will be built based on the following structure:

$$Y_e^T = (P_1, P_2, P_3, E_1, E_2, E_3, Q_1, Q_2, \Delta S_1, \Delta S_2, \Delta S_3). \quad (22)$$

Each water component will be observed three times, except for runoff $Q$, which is observed twice only. In the following Monte Carlo experiments, each one of these observations will be chosen equal to the true component in $X$ plus a random perturbing noise using Eq. (3).

The a priori information on observation uncertainty is given in Table 1 (right columns). These uncertainties are synthetic and have been chosen based on an educated physical guess. There is no given information on potential biases in the observations (bias will be present in the following Monte Carlo experiments). It will be seen in the following section that the actual observation errors used in the Monte Carlo experiments will be different from this a priori information in order to get closer to real conditions. For example, the observation error covariance matrix $R$ is well determined for the simulation of the experimental datasets, but our
<figure><img src="image_3.png" /><figcaption>FIG. 1. The four water cycle components (P, E, Q, and ΔS) over the Mississippi basin for 1984–2006. Dataset is from Pan et al. (2012).</figcaption></figure>

knowledge of **R** is only partial when testing the integration techniques. Tests have been conducted with a perfect knowledge of these observation errors; obviously, the retrieval results are better, but this is not a realistic situation.

b. Perturbation model

To test the integration methodologies presented in section 2, a Monte Carlo approach is applied to the dataset $\mathcal{D}$ of the previous section. A strategy could be to first perturb the data **X** once to obtain one synthetic observation **Y**<sub>ε</sub> and then test which method retrieves **X** better from **Y**<sub>ε</sub>. However, the results obtained using this strategy would be highly dependent on the chosen perturbation ε. The Monte Carlo approach consists of drawing a large number of perturbations ε in order to obtain reliable statistics on the quality of the integration methods.

The first very important component of a Monte Carlo experiment is therefore the statistical model for the perturbations ε. The perturbation model represents both errors in observation **Y**<sub>ε</sub> and uncertainties of our

knowledge of these observation errors. In this study, the goal is to obtain a model that represents a bias and variance error on each one of the observations in **Y**<sub>ε</sub>. The model also needs to represent the uncertainties on these observation errors, since most of the time, the user of the integration models has only a limited knowledge of these observation errors.

Let $S(t)$ be the time series of one water cycle component. The time series are perturbed using the following model:

$$S'(t) = S(t) + a(t) + b, \quad (23)$$

with $a(t)$ following a normal distribution $a(t) \sim \mathcal{N}(0, \sigma_a)$, where $\sigma_a \sim \mathcal{U}(0, +20)$, that is, a uniform distribution between 0 and +20, and $b \sim \mathcal{U}(-10, +10)$, also a uniform distribution. When the random errors $a(t) + b$ are such that the term $S(t) + a(t) + b$ is negative, the sampling is redone in order to avoid negative P, E, or Q components. The random bias $b$ is independent from the time index $t$, so it represents a random bias applied for the whole time series for each N Monte Carlo members.
<table><caption>TABLE 1. Statistics on the Mississippi and Niger water components and a priori information on observation uncertainty characteristics for observation $Y_e$.</caption><thead><tr><th rowspan="2">Water component (mm)</th><th colspan="2">Mississippi</th><th colspan="2">Niger</th><th colspan="3">Observation errors (mm)</th></tr><tr><th>Mean</th><th>STD</th><th>Mean</th><th>STD</th><th>Set 1</th><th>Set 2</th><th>Set 3</th></tr></thead><tbody><tr><td>P</td><td>64.7</td><td>23.8</td><td>51.3</td><td>54.4</td><td>15</td><td>16</td><td>11</td></tr><tr><td>E</td><td>48.1</td><td>36.2</td><td>35.2</td><td>16.1</td><td>13</td><td>8</td><td>11</td></tr><tr><td>Q</td><td>16.8</td><td>8.5</td><td>16.9</td><td>14.7</td><td>5</td><td>7</td><td></td></tr><tr><td>ΔS</td><td>-0.2</td><td>23.4</td><td>-0.8</td><td>30.0</td><td>13</td><td>18</td><td>11</td></tr></tbody></table>

The term $a(t)$ is a random Gaussian noise with no bias and a variance that is drawn uniformly from -10 and +10 for each time series. Different parameters for the two uniform distributions $U(0, +20)$ and $U(-10, +10)$ could have been chosen for each observation $Y_e$ (from $P_1$ to $\Delta S_3$). For simplicity of presentation, they are chosen here to be equal for each observation type. These noise characteristics constitute the synthetic errors of observations that will be used in the Monte Carlo experiment, but our a priori knowledge of these errors is imperfect and is provided in Table 1.

Other distributions of errors could have been chosen, for instance, a lognormal distribution for precipitation. Also, the error model can be chosen to be state dependent, for example, the random term is not the same for a low and a high runoff. Again, for simplicity, we have chosen to use a Gaussian distribution of the random errors because it is a good assumption when no other information on the uncertainties is available, and because often, errors become Gaussian when a first guess is used. A uniform distribution for parameters $\sigma_a$ and $b$ was chosen because it is the “no information” distribution in Bayesian statistics and because it represents well the type of a priori knowledge that is generally available about uncertainties.

### c. Monte Carlo algorithm

The Monte Carlo scheme for the estimation of the integration error statistics is described below to facilitate the understanding of the approach. It uses $N = 1000$ samples; the stability of the results to this number has been validated. The inputs are the datasets $D_{Mississippi}$ and $D_{Niger}$ (i.e., two $276 \times 4$ matrices for Mississippi and Niger) gathering the time series of the state variables $X$ for both basins. The Monte Carlo algorithm draws random samples $Y_e$ using the perturbation model. The outputs of the Monte Carlo are the retrieval statistics of $(X - X_a)$ for the various methodologies. We will test only six methodologies: SW, OI, NN, and for each one of them, their PF-filtered versions. No off-diagonal error in $R$ for the input datasets is considered in these experiments.

Data: Mississippi and Niger data: $D_{Mississippi}$ and $D_{Niger}$.

Result: Retrieval accuracies of the integration methodologies.

Define $R$ and $B$ (i.e., a priori information on observation and first guess errors);

For $i = 1$ to $N$ (i.e., number of Monte Carlo samples, 1000) do

For $j = 1$ to $J$ (i.e., number of observations, 276) do

Draw bias $b$ from $U(-10, +10)$;

Draw $\sigma_a$ from $U(0, +20)$;

For $t = 1$ to $T$ (number of time steps) do

Draw $a(t)$ from $N(0, \sigma_a)$;

End

Obtain observation noises $\epsilon_j(t) = a(t) + b$ for each time step $t$.

End

Obtain $Y_e = A \cdot X + \epsilon(t)$ for each time step.

Estimate $X_a^{SW}$ using Eq. (7);

Estimate $X_a^{OI}$ using Eq. (17);

Estimate $X_a^{NN}$ using Eq. (19);

Apply the PF filtering to each estimate and obtain:

$X_a^{SW+PF}, X_a^{OI+PF}, X_a^{NN+PF}$;

End

Perform the $(X_a - X)$ root-mean-square (RMS) statistics on the pooled ($N$ samples) dataset for the six methods.

## 4. Comparison of the integration methods

In this section, the integration methods presented in section 2 are tested on a dataset of real data but are artificially perturbed by a random noise following the Monte Carlo experiment described in the previous section.

### a. Filters

The first filter considered here is the PF $K_{PF}$ of Eq. (21). Figure 2 represents this $4 \times 4$ nonsymmetric matrix. It can be seen that $K_{PF}$ is close to the identity matrix, but the two more uncertain variables, $\Delta S$ and even more $P$, have a lower diagonal amplitude, which means that
<figure><img src="image_4.png" /><figcaption>FIG. 2. The values indicate how the weighting is done in the PF filter between observational water components and resulting components with closure.</figcaption></figure>

these variables can be modified more than the two others to reach the closure. The closure is possible thanks to the off-diagonal terms. This closure filter can be applied a posteriori to the SW, CL, OI, and NN filters.

In Fig. 3, the 4 × 11 matrix filters of the SW and OI methods are represented without and with the PF. These filters are identical for Mississippi and Niger since they are not dependent on the state variables **X** or observations **Y**<sub>e</sub>. They would change if our a priori information on observation uncertainty in Table 1 was different for the two basins. It can be seen that the SW filter has no off-diagonal terms: the three observations for precipitation are used only for the precipitation, and similarly for evapotranspiration, runoff, and terrestrial water storage change. This is a structural constraint. Again, the observations with less uncertainties have a larger weight. The application of the PF filter on the SW solution keeps a structure that is relatively close to the SW filter, but off-diagonal terms are introduced to obtain the water budget closure. The application of the PF filter on the OI solution has a relatively small impact because a closure constraint was already introduced in the OI filter [Eq. (17)]. However, the weights on the precipitation data are particularly low, meaning that the determination of the precipitation in **X**<sub>a</sub><sup>OI</sup> is almost as much determined by the precipitation observations in **Y**<sub>e</sub> as by the other water cycle components and the closure constraint.

Figure 4 represents the Jacobians $(\partial g / \partial Y_i)(Y_e)$ of the neural network integration model. They are the analog of the linear filters SW or OI since these Jacobians represent a linearization of the nonlinear NN function g. As such, they can be directly compared to Fig. 3. The

Jacobians are dependent on the actual observations $\mathbf{Y}_e$, so they are represented here for both Mississippi and Niger, and for January and July means (i.e., monthly averages over the 23 years of data). It can be noted that January and July Jacobians are not identical. Mississippi and Niger Jacobians differ too. This means that the NN is able to optimize its weighting of the observations in $\mathbf{Y}_e$ to the particular location and time. However, the overall structure of the filters is similar and close to the filters of Fig. 3, which is sign of robustness of the integration methodologies. The closure constraint introduced with PF has almost no impact on the NN filters (not shown) because the closure is already good with the NN integration solution.

b. Integration results

Figure 5 represents the time series of the precipitation, evapotranspiration, runoff, water storage change, and terrestrial water budget for year 1984 for the Niger basin. There are three observations for each water component (red). These observations are deviations from the target cycle (dashed black line). The SW solution (green) is a weighting of the observations. It can be seen in the time series of the budget that this weighting of the observations is already a good improvement, even if no closure constraint is used, since the residuals are much lower. The NN estimate (blue) is an even better solution: the residuals are, in this case, negligible. It confirms that the NN is able to obtain the closure even if this constraint was only implicitly imposed in the learning dataset $\mathcal{D}$.

The RMS errors of the retrieval of the four water cycle components and the residuals of the terrestrial water budget are represented in Table 2 for the Mississippi and Niger basins, for the whole 1984–2006 record. These estimates result from the Monte Carlo experiments described in section 3. The STD of each variable is represented first for comparison purposes. Typical errors for the observations are also represented: set 1 is for ($P_1, Q_1, \Delta S_1$) in Eq. (22) and set 2 is for ($P_2, Q_2, \Delta S_2$). The RMS errors for set 1 are the same for both basins since no location dependency has been introduced in the observation noise. However, for each water component, errors from the two sets are different.

The results of this table are also represented in Fig. 6 for Mississippi (continuous lines) and for Niger (dashed lines). Only the SW, SW+PF, and NN results are presented since the other integration methods are comparable to one or the other of these three methods. Each integration technique is better than the use of a single set of observations. This is expected, of course, since the combination of multiple observations is better than the use of only one. SW is a very good improvement over
<figure><img src="image_5.png" /><figcaption>FIG. 3. Filters for SW, SW + PF, OI2, and OI3 = OI1 + PF. The values in these matrices indicate how the weighting is performed by each filter on each input dataset, in order to obtain the results integrated dataset of the four water cycle components.</figcaption></figure>

the use of a single set of observations; the RMS errors can be divided by a factor of 2. However, the resulting budget is still not acceptable (RMS about 9.25 mm for Mississippi). The introduction of the postprocessing filtering PF to SW solves this problem: the budget reaches closure. In addition, the introduction of the closure constraint has not degraded the retrievals of SW; it has actually improved it significantly (Ko and Bitmead 2007), for example, from 5.38 to 4.36 for precipitation. This means that the closure constraint is good a priori information and that its use in the integration process is very valuable.

The RMS errors of OI1 are almost identical to SW. This is surprising since more information is used in the OI integration techniques, that is, the observing operator **A**. However, this information is relatively limited in this experiment because of the simplicity of the **A** structure. Its use in the integration might only be useful for more sophisticated models such as when **A** is the linearization of the radiative transfer in numerical weather prediction centers for traditional weather forecasting. The introduction of the closure constraint in the OI2 estimation has a very strong positive impact in the water budget closure that decreases from 9.25 for the OI1 estimate, to 0.01 with OI2, again for the Mississippi case. This demonstrates the ability of the OI technique to utilize external constraints such as the closure. Again, the introduction of the closure constraint also improves the estimate of the water components. Note that results for the SW + PF and the OI2 solution are almost identical. Furthermore, the OI1 + PF solution, which includes the PF filtering, performs as well as the OI2 estimate. This means that the introduction of the closure can be done simultaneously with the observations or as a postprocessing step.

The NN approach appears to be quite a good solution. It outperforms the other methods in terms of retrieval accuracy. Furthermore, it is possible to confirm that the closure constraint imposed only implicitly in the dataset $\mathcal{D}$ used to train the NN is respected in the NN outputs since the NN budget has RMS error close to zero. The use of the PF filter does not improve the solution: the retrieval uncertainties are the same and the closure of the budget is only slightly reduced. The fact that the NN follows the closure well using only an implicit constraint could be surprising. However, the PF processing to explicitly impose the closure (section 2g) is just a 4 × 4 matrix close to identity (Fig. 2), so if the NN “learns” during the training that such postprocessing is useful to obtain optimal retrievals, it is easy for the NN to represent it in its model g. Another remark is that these NN experiments have been performed without any regularization imposing constraints on the weight amplitude (section 2f), so the fact that the results appear robust means that the problem is well posed.

## 5. Conclusions, discussion, and perspectives

A multitude of satellite observations are available to monitor Earth's surface and atmosphere at various wavelengths, from various orbits (geostationary or polar orbiters), and with different space and time samplings. In general, these numerous observations are used independently to derive one particular geophysical variable. These observations have resulted in a multiplicity of datasets for the same geophysical variable. Each dataset has its own advantages and limitations [e.g., one soil moisture dataset can perform better over semiarid regions, another can be better for vegetated areas
<figure><img src="image_6.png" /><figcaption>FIG. 4. NN Jacobians ($\partial g/\partial Y_i$) ($Y_i$) (analog of the linear filters of Fig. 3) for (from top to bottom) Mississippi in January and July and Niger in January and July.</figcaption></figure>

(Liu et al. 2011), so users have difficulties in choosing one dataset among the others.

This situation is close to what happens in the climate modeling community, where many climate models are available and each one is offering valuable information. To solve this issue and exploit optimally all the information available about the future climate, the scientific community has averaged the forecasts from the major climate models to obtain 1) a more reliable estimate of the future climate and 2) a quantification of model dispersion as a surrogate for climate forecast uncertainty. In a similar way, the Earth observation community needs to organize a better framework to combine, as well as possible, all the available satellite observations. International programs such as the GEWEX program or space agencies such as the European Space Agency (ESA) now focus on the production of unified and coherent global satellite datasets.

It has been shown in the literature that, if possible, it is optimal to combine the satellite measurements before the retrieval so that the potential synergy among the observations could be exploited by the retrieval algorithm itself (Aires et al. 2012). However, this is not always possible in practice, so there is a need for methods combining a posteriori multiple geophysical datasets from independent retrievals. The integration solutions

presented in this study intend to optimize this a posteriori combination using a priori information on the observation uncertainty estimates. Four methodologies have been introduced in order to obtain optimized terrestrial water cycle component, in particular, estimates of precipitation, evapotranspiration, runoff, and terrestrial water storage change: simple weighting (SW), the constrained linear (CL) method, the optimal interpolation (OI) method, and the neural network (NN) technique. All of them can be combined, if necessary, with a postprocessing filtering (PF) step to impose the water budget closure (Pan and Wood 2006). These integration methodologies are able to exploit the synergy of multiple observational datasets from satellite observations or in situ measurements (Aires 2011; Aires et al. 2011).

The presented methodologies are also able to use a physical constraint that links the four water cycle components, for example, the closure of the terrestrial water cycle budget. This means that, first, multiple water cycle component estimates are optimized at once, and, second, that these components become more coherent with each other. The Monte Carlo numerical experiments show that the introduction of this additional physical constraint actually helps the retrieval process. This a priori physical constraint is independent from any model, so the final dataset of optimized water cycle
<figure><img src="image_7.png" /><figcaption>FIG. 5. Time series of the *P*, *E*, *Q*, *ΔS*, and terrestrial water budget for year 1984 for the Niger basin. Observations are in red, the SW estimates are in green, the NN estimate in blue, and the target is dashed black.</figcaption></figure>

components is model independent<fnref n="3" /> and it can be used to validate or calibrate hydrological, surface, or climate models. Among the four methods, the solutions seem to have equivalent performances once the closure constraint is used. The introduction of the closure constraint is very important; it actually improves not only the budget but also the estimation of the four components.

The perspectives for this work are many fold. In a companion paper, various global datasets providing multiple estimates of each terrestrial water cycle component will be gathered. The integration approaches will be used to obtain a global optimized and coherent dataset of the four water cycle components. The weighting of the various initial observation datasets will be investigated. Some of the techniques have their weighting dependent on the state variables of the

hydrologic system or on the observation uncertainty. These variations of the weighting will be investigated globally. Furthermore, precipitation, evaporation, and water storage changes can be observed globally, but runoff is only measured in some local stations. To obtain an optimized and global dataset of (*P*, *E*, *ΔS*), the local weights in the runoff measurement locations need to be interpolated to the locations without runoff measurements. Another way to improve the water cycle dataset is to introduce even more constraints. For instance, links to other water cycle components such as soil moisture (Liu et al. 2011; Kolassa et al. 2013) or inundations (Prigent et al. 2007) could be identified. An example of possible links could be the correlation between two water cycle variables. Then, the introduction of these links as new constraints in the integration scheme could be considered. This would require new technical developments.

Another very important topic concerns the estimation of errors for the input datasets, before integration. We have seen in this paper that this information has a strong impact on most integration methods. Unfortunately,

<footnote n="3">³ Note that the techniques presented here could be used together with an FG estimate from a numerical model (like in variational assimilation) but, in this case, their potential use for model validation would vanish.</footnote>
<table><caption>TABLE 2. RMSE statistics from the Monte Carlo experiments for the retrieval of *P*, *E*, *Q*, Δ*S*, and the resulting water cycle budget for the Mississippi and Niger basins. The STD of each variable is represented for comparison purposes. Integration methodologies are SW, SW + PF, OI1, OI2, OI1 + PF, NN, and NN + PF. Statistics are performed for the whole 1984–2006 period.</caption><thead><tr><th></th><th><i>P</i></th><th><i>E</i></th><th><i>Q</i></th><th>Δ<i>S</i></th><th>Budget</th></tr></thead><tbody><tr><td><strong>Mississippi</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>STD</td><td>23.71</td><td>36.11</td><td>8.40</td><td>23.40</td><td>0.00</td></tr><tr><td>Obs. set 1</td><td>8.90</td><td>9.56</td><td>4.07</td><td>8.10</td><td>15.94</td></tr><tr><td>Obs. set 2</td><td>10.69</td><td>7.44</td><td>4.97</td><td>11.04</td><td>17.85</td></tr><tr><td>SW</td><td>5.38</td><td>4.78</td><td>3.18</td><td>4.78</td><td>9.25</td></tr><tr><td>SW + PF</td><td>4.36</td><td>4.10</td><td>2.99</td><td>4.17</td><td>0.00</td></tr><tr><td>OI1</td><td>5.38</td><td>4.80</td><td>3.18</td><td>4.75</td><td>9.25</td></tr><tr><td>OI2</td><td>4.36</td><td>4.12</td><td>2.99</td><td>4.16</td><td>0.01</td></tr><tr><td>OI1 + PF</td><td>4.36</td><td>4.12</td><td>2.99</td><td>4.16</td><td>0.00</td></tr><tr><td>NN</td><td>4.18</td><td>3.76</td><td>2.65</td><td>3.83</td><td>0.05</td></tr><tr><td>NN + PF</td><td>4.18</td><td>3.76</td><td>2.65</td><td>3.83</td><td>0.00</td></tr><tr><td><strong>Niger</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>STD</td><td>54.31</td><td>16.09</td><td>14.65</td><td>29.90</td><td>0.00</td></tr><tr><td>Obs. set 1</td><td>9.03</td><td>9.50</td><td>4.13</td><td>7.93</td><td>15.92</td></tr><tr><td>Obs. set 2</td><td>10.88</td><td>7.42</td><td>4.95</td><td>10.77</td><td>17.71</td></tr><tr><td>SW</td><td>5.40</td><td>4.94</td><td>3.16</td><td>4.82</td><td>9.29</td></tr><tr><td>SW + PF</td><td>4.38</td><td>4.27</td><td>2.97</td><td>4.23</td><td>0.00</td></tr><tr><td>OI1</td><td>5.40</td><td>4.97</td><td>3.16</td><td>4.77</td><td>9.28</td></tr><tr><td>OI2</td><td>4.38</td><td>4.29</td><td>2.97</td><td>4.21</td><td>0.01</td></tr><tr><td>OI1 + PF</td><td>4.38</td><td>4.29</td><td>2.97</td><td>4.21</td><td>0.00</td></tr><tr><td>NN</td><td>3.80</td><td>3.09</td><td>1.92</td><td>3.60</td><td>0.02</td></tr><tr><td>NN + PF</td><td>3.80</td><td>3.09</td><td>1.92</td><td>3.60</td><td>0.00</td></tr></tbody></table>

error information is very difficult to obtain from data producers. In our experiments, Gaussian and uniform distributions have been used, but more complex distributions might be more adequate (such as a lognormal distribution for precipitation errors). The distribution of errors might also be dependent on the state (e.g., bigger relative errors when the signal to measure is smaller). An important point also concerns the correlation of errors among the input datasets. Satellite datasets use often common features (like the same radiative transfer model or a same auxiliary dataset used for the retrieval) so correlation of errors might appear. Even if no common feature is used among the various input datasets, the cases with larger errors for one dataset might be the cases with larger errors for another dataset; this would create error correlations. If such complex information on the errors was available, it could be used in the integration method. It might be possible to use Bayesian statistics to estimate “implicitly” the uncertainties of each input datasets and their relationship during the integration step. This will be the subject of a forthcoming study.

*Acknowledgments.* We thank the ESA (European Space Agency) and in particular Diego Fernandez for funding the project “Watchful: STSE Water Cycle

<figure><img src="image_8.png" /><figcaption>FIG. 6. RMSE of the retrieved water component (*P*, *E*, *Q*, and Δ*S*) using method SW (red), SW + PF (green), and NN (blue). The other integration methods are similar to one of these curves. The resulting water cycle budget is also represented. Results are for the whole 1984–2006 time series. Mississippi basin is in continuous lines and Niger in dashed lines.</figcaption></figure>

Feasibility” (4000107122/12/I-NB; http://watchful.geo.tuwien.ac.at) and Marcela Doubkova for managing this project. We also thank colleagues Wolfgang Wagner, Richard Kidd, Wouter Dorigo, Patrick Matgen, and Paul Bates from the Watchful project for interesting related discussions and in particular Stefan Schlaffer, who co-ordinated this project. We are also grateful to Ming Pan and colleagues from Princeton University for providing the dataset used in this study.

## REFERENCES

Aires, F., 2011: Measure and exploitation of multisensor and multiwavelength synergy for remote sensing: 1. Theoretical considerations. *J. Geophys. Res.*, **116**, D02301, doi:10.1029/2010JD014701.

—, and W. Rossow, 2003: Inferring instantaneous, multivariate and nonlinear sensitivities for the analysis of feedback processes in a dynamical system: Lorenz model case study. *Quart. J. Roy. Meteor. Soc.*, **129**, 239–275, doi:10.1256/qj.01.174.

—, M. Schmitt, N. Scott, and A. Chédin, 1999: The weight smoothing regularization for MLP for resolving the input contribution’s errors in functional interpolations. *IEEE Trans. Neural Networks*, **10**, 1502–1510, doi:10.1109/72.809096.

—, C. Prigent, and W. B. Rossow, 2004a: Neural network uncertainty assessment using Bayesian statistics with application to remote sensing: 3. Network Jacobians. *J. Geophys. Res.*, **109**, D10305, doi:10.1029/2003JD004175.

—, —, and —, 2004b: Neural network uncertainty assessment using Bayesian statistics: A remote sensing application. *Neural Comput.*, **16**, 2415–2458, doi:10.1162/0899766041941925.

—, F. Bernardo, H. Brogniez, and C. Prigent, 2010: Calibration for the inversion of satellite observations. *J. Appl. Meteor. Climatol.*, **49**, 2458–2473, doi:10.1175/2010JAMC2435.1.
—, M. Paul, C. Prigent, B. Rommen, and M. Bouvet, 2011: Measure and exploitation of multisensor and multiwavelength synergy for remote sensing: 2. Application to the retrieval of atmospheric temperature and water vapor from MetOp. *J. Geophys. Res.*, **116**, D02302, doi:10.1029/2010JD014702.

—, O. Aznay, C. Prigent, M. Paul, and F. Bernardo, 2012: Synergistic multi-wavelength remote sensing versus a posteriori combination of retrieved products: Application for the retrieval of atmospheric profiles using MetOp-A. *J. Geophys. Res.*, **117**, D18304, doi:10.1029/2011JD017188.

Azarderakhsh, M., W. B. Rossow, F. Papa, M. Norouzi, and R. Khanbilvardi, 2011: Diagnosing water variations with the Amazon basin using satellite data. *J. Geophys. Res.*, **116**, D24107, doi:10.1029/2011JD015997.

Bishop, C., 1996: *Neural Networks for Pattern Recognition*. Clarendon Press, 482 pp.

De Geeter, J. D., H. V. Bruseel, and J. D. Schutter, 1997: A smoothly constrained Kalman filter. *IEEE Trans. Pattern Anal. Mach. Intell.*, **19**, 1171–1177, doi:10.1109/34.625129.

Fernandez-Prieto, D., P. van Oevelen, Z. Su, and W. Wagner, 2012: Advances in earth observation for water cycle science. *Hydrol. Earth Syst. Sci.*, **16**, 543–549, doi:10.5194/hess-16-543-2012.

Foley, A. M., and Coauthors, 2013: Evaluation of biospheric components in Earth system models using modern and palaeo- observations: The state-of-the-art. *Biogeosciences*, **10**, 8305– 8328, doi:10.5194/bg-10-8305-2013.

Hayward, S., 1998: Constrained Kalman filter for least-squares estimation of time-varying beamforming weights. *Mathematics in Signal Processing IV*, J. McWhirter and I. Proudler, Eds., Oxford University Press, 113–125.

Kalman, R., 1960: A new approach to linear filtering and prediction problems. *Trans. ASME, 82D*, 35–45.

Kalnay, E., and Coauthors, 1996: The NCEP/NCAR 40-Year Re-analysis Project. *Bull. Amer. Meteor. Soc.*, **77**, 437–471, doi:10.1175/1520-0477(1996)077<0437:TNYRP>2.0.CO;2.

Ko, S., and R. R. Bitmead, 2007: State estimation for linear systems with state equality constraints. *Automatica*, **43**, 1363–1368, doi:10.1016/j.automatica.2007.01.017.

Kolassa, J., F. Aires, J. Polcher, C. Prigent, C. Jimenez, and J. Pereira, 2013: Soil moisture retrieval from multi-instrument observations: Information content analysis and retrieval methodology. *J. Geophys. Res. Atmos.*, **118**, 4847–4859, doi:10.1029/2012JD018150.

Liu, Y., and Coauthors, 2011: Developing an improved soil moisture dataset by blending passive and active microwave satellite-based retrievals. *Hydrol. Earth Syst. Sci.*, **15**, 425–436, doi:10.5194/hess-15-425-2011.

Pan, M., and E. F. Wood, 2006: Data assimilation for estimating the terrestrial water budget using a constrained ensemble Kalman filter. *J. Hydrometeor.*, **7**, 534–547, doi:10.1175/JHM495.1.

—, A. K. Sahoo, T. J. Troy, R. K. Vinukollu, J. Sheffield, and E. F. Wood, 2012: Multisource estimation of long-term terrestrial water budget for major global river basins. *J. Climate*, **25**, 3191–3206, doi:10.1175/JCLI-D-11-00300.1.

Porrill, J., 1988: Optimal combination and constraints for geometrical sensor data. *Int. J. Robotics Res.*, **7**, 66–77, doi:10.1177/027836498800700606.

Prigent, C., F. Papa, F. Aires, W. Rossow, and E. Matthews, 2007: Global inundation dynamics inferred from multiple satellite observations, 1993–2000. *J. Geophys. Res.*, **112**, D12107, doi:10.1029/2006JD007847.

Rodgers, C., 2000: *Inverse Methods for Atmospheric Sounding: Theory and Practice*. Series on Atmospheric, Oceanic and Planetary Physics, Vol. 2, World Scientific Publishing, 240 pp.

Rumelhart, D., G. Hinton, and R. Williams, 1986: Learning internal representations by error propagation. *Foundations, Vol. 1, Parallel Distributed Processing: Explorations in the Microstructure of Cognition*, MIT Press, 318–362.

Sahoo, A. K., M. P. Pan, T. J. Troy, R. K. Vinukollu, J. Sheffield, and E. F. Wood, 2011: Reconciling the global terrestrial water budget using satellite remote sensing. *Remote Sens. Environ.*, **115**, 1850–1865, doi:10.1016/j.rse.2011.03.009.

Sheffield, J., C. Ferguson, T. Troy, E. Wood, and M. McCabe, 2009: Closing the terrestrial water budget from satellite remote sensing. *Geophys. Res. Lett.*, **36**, L07403, doi:10.1029/2009GL037338.

Simon, D., and T. L. Chia, 2002: Kalman filtering with state equality constraints. *IEEE Trans. Aerosp. Electron. Syst.*, **38**, 128–136, doi:10.1109/7.993234.

Tarantola, A., 1987: *Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation*. Elsevier, 613 pp.

Thome, K., 2004: In-flight intersensor radiometric calibration using vicarious approaches. *Post-Launch Calibration of Satellite Sensors*, S. A. Morain and A. M. Budge, Eds., ISPRS Book Series, Vol. 2, Taylor and Francis, 93–102.