Journal cover Journal topic
Hydrology and Earth System Sciences An interactive open-access journal of the European Geosciences Union
Journal topic
Hydrol. Earth Syst. Sci., 23, 465-491, 2019
https://doi.org/10.5194/hess-23-465-2019
Hydrol. Earth Syst. Sci., 23, 465-491, 2019
https://doi.org/10.5194/hess-23-465-2019

Research article 25 Jan 2019

Research article | 25 Jan 2019

# Integrating multiple satellite observations into a coherent dataset to monitor the full water cycle – application to the Mediterranean region

Integrating multiple satellite observations into a coherent dataset to monitor the full water cycle
Victor Pellet1,2, Filipe Aires1,2, Simon Munier3, Diego Fernández Prieto5, Gabriel Jordá4, Wouter Arnoud Dorigo6, Jan Polcher7, and Luca Brocca8 Victor Pellet et al.
• 1Laboratoire d'Études du Rayonnement et de la Matière en Astrophysique et Atmosphères, Observatoire de Paris, Paris, France
• 2Estellus, Paris, France
• 3Centre National de Recherches Météorologiques, Toulouse, France
• 4Institut Mediterrani d'Estudis Avançats, University of the Balearic Islands, Mallorca, Spain
• 5European Space Agency, Frascati, Italy
• 6Department of Geodesy and Geoinformation, Technical University, Vienna, Austria
• 7Laboratoire de Météorologie Dynamique, Paris, France
• 8Istituto di Ricerca per la Protezione Idrogeologica, Perugia, Italy
Abstract

The Mediterranean region is one of the climate hotspots where the climate change impacts are both pronounced and documented. The HyMeX (Hydrometeorological Mediterranean eXperiment) aims to improve our understanding of the water cycle from the meteorological to climate scales. However, monitoring the water cycle with Earth observations (EO) is still a challenge: EO products are multiple, and their utility is degraded by large uncertainties and incoherences among the products. Over the Mediterranean region, these difficulties are exacerbated by the coastal/mountainous regions and the small size of the hydrological basins. Therefore, merging/integration techniques have been developed to reduce these issues. We introduce here an improved methodology that closes not only the terrestrial but also the atmospheric and ocean budgets. The new scheme allows us to impose a spatial and temporal multi-scale budget closure constraint. A new approach is also proposed to downscale the results from the basin to pixel scales (at the resolution of 0.25). The provided Mediterranean WC budget is, for the first time, based mostly on observations such as the GRACE water storage or the netflow at the Gibraltar Strait. The integrated dataset is in better agreement with in situ measurements, and we are now able to estimate the Bosporus Strait annual mean netflow.

1 Introduction

The Mediterranean region is one of the main climate change hotspots (IPCC2014): its sensitivity to global change is high and its evolution remains uncertain. Its role in the evolution of the global ocean (i.e., mainly salinization and warming), as well as the socio-economic consequences it has for surrounding countries, stress the need for monitoring its water resource. Analyzing the water cycle (WC) and the exchanges among its terrestrial, atmospheric and oceanic branches is critical for estimating the availability of the water in the Mediterranean region. Most previous studies use model outputs or reanalysis , and in situ data networks are too sparse and irregular. A recent paper reviewed the literature on the analysis and quantification of the Mediterranean water budget using observation, model outputs and reanalyses. The WC components are estimated, but their uncertainties remain high. Recommendations are made to increase our use of Earth observations (EO) data, in a coordinated way. EO allow for the monitoring of the WC over long time records, in particular in regions with a low number of in situ stations. But the use of EO data for WC monitoring remains a challenge due to (1) the multiplicity of datasets for the same geophysical parameter, (2) the EO uncertainties (systematic and random errors), and (3) the inconsistency between datasets (for the same component or among the components of the WC). In , EO are used to monitor the WC over the Mediterranean region: it is shown that the WC budget is not closed and that some integration technique should be used to optimize them.

Several approaches have been considered in order to optimize EO datasets at the global scale, for the WC analysis. The features of some “integration” methods presented in the following are synthesized in Table 1.

The “Princeton” approach presented first a work in which they aimed to close the water balance using EO products. In this work, EO datasets such as precipitation were assimilated into a land surface model (the Variable Infiltration Capacity, VIC) using the combination of a Kalman filter and a closure constraint (see Table 1). The resulting “analysis” dataset is not a pure EO product since the VIC model is largely used. In fact, the authors show that the Kalman filtering plus the closure constraint is equivalent to a traditional Kalman filtering, and then to the application of an independent post-filtering that constrains the closure . This post-filtering acts by redistributing the budget residuals within each water component based on the uncertainties of each EO source. Several papers have been published based on this approach . For instance, in , two different precipitation datasets were used over the Mississippi basin. Evapotranspiration was calculated using a revised Penman–Monteith formulation and changes in water storage were estimated from GRACE. For comparison, land surface model outputs, reanalyses data and in situ discharge measurements were used too. The authors concluded that a positive bias of the precipitation datasets leads to an overestimation of the discharge component when the estimation relies on EO data. Meanwhile, the land surface model shows a high degree of agreement with in situ data. The analysis also highlights the importance of error characterization in the individual WC components. relaxed the closure constraint during the assimilation. This is an important feature because tight closure constraint can result in high-frequency oscillations in the resulting combined dataset. A large constraint is used in our approach (see Table 1).

The NASA-NEWS project – the project aims at a better characterization of the WC using EO data. The first step was to improve the coherency of the satellite retrievals, and then to gather the EO dataset and calibrate them. Some information about the uncertainties of the EO datasets was gathered from the data producers, but this information cannot be straightforwardly used further in the integration process since its evaluation is not homogeneous but product-dependent. The WC budget can be closed using the satellite datasets . However, this closure is obtained at the global and annual scales only, and residuals are still significant at regional and monthly scales. then use interpolation for a monthly closure. Closing the budget at the global scale was a first step, and closure must now be obtained at finer spatial and temporal scales to monitor more precisely the distribution of the water components as the EO data are designed to do. In , the storage terms (e.g., groundwater storage) had no significant change when considering annual and global means. This hypothesis was then used at the monthly scale with an optimized interpolation scheme to relax the storage change at the monthly scale. This approach translates into a quadratic quality criterion where storage and fluxes terms are minimized for annual means, at the global scale (see Table 1). One interesting feature in this approach is that both the water and energy cycles were considered simultaneously in the assimilation, taking into account the physical link between the two cycles through the latent heat flux.

Table 1The three main initiatives for budget closure constraint and their technical differences. (In the third column, bold font indicates the new features of the methodology presented in this article.) Subscripts are “l” for land and “o” for ocean; both include the atmosphere. All notations are summarized in Table A1. 1  recently developed a WC-VIC assimilation scheme at the 0.5 pixel scale; 2  used a two-step integration method with annual closure simply downscaled at the monthly scale, plus a Lagrange interpolation for closure relaxation; 3 see for details and hypotheses.

The ESA water cycle initiative – in the context of the ESA WATCHFUL project on water budget closure, described several methodologies (Table 1) to integrate different hydrological datasets with a budget closure constraint. No surface or atmospheric models were used in these integration methods, making the obtained product interesting for model calibration and validation. One of the proposed methods, the so-called simple weighting + post-filtering (SW + PF), was applied by over the Mississippi basin, using satellite datasets for precipitation, evapotranspiration and water storage, and gauge observations for river discharge. After applying a budget closure constraint at the basin scale, the integrated components were compared to various in situ observations, showing good performances of the method. One of the main limitations is the dataset availability of the in situ river discharges. Another concern was the downscaling of the basin closure constraint to the pixel scale. A closure correction model (CCM) is a calibration of the EO that was developed based on the integrated product as the reference . It allows correction of each dataset independently to greatly reduce the budget residuals. This calibration was applied over the basins where river discharges are available and extended to the global scale using an index characterizing the various surface types . This type of post-processing step is anchored in the integration approach, but it can be applied to long time records, at any time or spatial resolution. It can even allow for the reconstruction of missing estimates.

In this paper, we propose several improvements of this line of research. In particular, we propose to close the WC budget not only over land, but also over ocean and in the atmosphere. Furthermore, the budget closure constraint is used simultaneously at different spatial (basin and sub-basin) and temporal (monthly and annual) scales. A new spatial interpolation scheme is proposed to downscale the basin-scale closure constraint to the pixel scale. This new framework is applied to the Mediterranean basin to provide an updated WC budget.

Section 2 presents the study domain and introduces the datasets used in the following. The integration approaches are described together with the other combination techniques in Sect. 3. Section 4 presents the evaluation metrics for the integrated product: its ability to close the WC and its validation with in situ data at the sub-basin or pixel scale. Section 5 presents the WC analysis for the period 2004–2009 using our resulting integrated dataset. Finally, Section 6 concludes the analysis and presents some perspectives. All notations used in the following are summarized in Table A1 in the Appendix.

2 Case study and datasets

This section presents the spatial domain and the datasets used in this study. Table A2 in the Appendix summarizes the main characteristics of these products and more details can be found in . All products have different temporal extents but share a common coverage period 2004–2009.

## 2.1 Mediterranean region

The study domain is represented in Fig. 1. It is the catchment basin of the whole Mediterranean Sea drainage area, computed from each coastal pixel, including all rivers that flow into the sea. Basins have been computed using a hydrographic model at a spatial resolution of 0.25. The resolution of the hydrographic model used to compute the land–sea mask or catchment basin may have an impact on the spatial-average estimates and then on the WC budget residual. This area uncertainty is taken into account in the relaxation of the closure constraint at sub-basin scale (see Table 1). The Mediterranean Sea area (including the Black Sea) is 3.0 million km2, and its drainage area is more than 5 million km2.

Sub-basins were introduced in . They facilitate the analysis of local climate and specific hydrological features. The Mediterranean Sea and the terrestrial sub-basins used in the following are defined as

• the western Maghreb mainly based on the Atlas mountain discharge (MA-DZ-TN);

• the Nile basin and Libyan coast characterizing a Saharan and sub-Saharan climate (LY-EG);

• the Spanish coasts and Pyrénées (ES-Pyr);

• the French, Italian and Adriatic Sea coasts, carrying freshwater from the Alps and the Balkan mountains (Alp-IT-ADR);

• the eastern part of the Mediterranean Sea, Greece, Turkey and Israel (GR-TR-IL); and

• the whole Black Sea drainage catchment, Bulgaria, Georgia, Romania, Russia, Turkey, Ukraine, Slovakia, Hungary, Austria, Slovenia, Bosnia and Serbia (BLS).

In the current study, even if the closure method (PF) is applied over the LY-EG sub-basin, the high uncertainty of the Nile discharge and its particular climate (African monsoon) as well as anthropogenic conditions (most of its water is used for irrigation) make this sub-basin really different from the other sub-basins . Therefore the closure is ensured for the Nile sub-basin, but no spatial extension will be extrapolated over the LY-EG or the south (see Sect. 3).

Figure 1Region of interest. Sub-basins have been computed using a hydrological model , and rivers are from HydroShed (http://www.hydrosheds.org/, last access: 24 January 2019). See text for the definition of the sub-basins.

## 2.2 Original EO datasets

The datasets presented in this section will be used in the integration process. Most of them are satellite products and are commonly used for studying the WC. In order to integrate them, the datasets have been projected onto a common 0.25 spatial resolution grid and re-sampled at the monthly scale.

Precipitation (P) – four satellite-based datasets have been selected. Two are gauge-calibrated products: the Tropical Rainfall Measuring Mission Multi-satellite Precipitation Analysis (TMPA, 3B42 V7) presented in and the Global Precipitation Climatology Project (GPCP, v2) introduced by . Two are uncalibrated products: have unveiled the NOAA CPC Morphing Technique (CMORPH, v1) and developed the Precipitation Estimation from Remote Sensing Information using Artificial Neural Network (PERSIANN, v1). In this study, we use a mix of gauged/ungauged-calibrated precipitation datasets. This choice is motivated by the goal of preserving the original EO spatial pattern where limited gauge density in some areas may corrupt the signal during the gauge-calibration process (in TMPA and GPCP products).

Evapotranspiration (E) – three satellite-based products were chosen to describe evapotranspiration over land: the Global Land Evaporation Amsterdam Model (GLEAM-V3B, Martens et al.2017; Miralles et al.2011); the MODIS Global Evapotranspiration Project (MOD16, Mu et al.2011); and the Numerical Terradynamic Simulation Group product (NSTG, Zhang et al.2010).

Two products were chosen for the evaporation over the sea: the Objectively Analyzed air-sea Fluxes for Global Oceans (OAflux, Sun et al.2003); and The Global Energy and Water Cycle Exchanges Project product (GEWEX-Seaflux, Curry et al.2004).

Water storage change (ΔS) – the terrestrial and seawater storage datasets are all derived from the GRACE mission. The estimates of water storage implicitly include the underground water. Four satellite datasets are based on the spherical decomposition of GRACE measurement: the Jet Propulsion Laboratory (JPL, Watkins and Yuan2014) product; the Centre for Space Research (CSR, Bettadpur2012) product; the German Research Centre for Geoscience (GFZ, Dahle et al.2013) product; and the land-only product from the Groupe de Recherche de Géodésie Spatiale (GRGS, Biancale et al.2005). One extra solution based on the JPL-MASCONS decomposition of GRACE measurement is also used in this work. In order to compute the monthly change in water storage, we applied a centered derivative smoothing filter: [5∕24 3∕8 $-\mathrm{3}/\mathrm{8}$ $-\mathrm{5}/\mathrm{24}$] . The filter is a slightly smoother version of the filter [1∕8 1∕4 $-\mathrm{1}/\mathrm{4}$ $-\mathrm{1}/\mathrm{8}$] presented by . It has been compared with several other filters (results not shown). The chosen filter is a good compromise between its smoothing (that suppresses information) and its ability to de-noise the time series.

Discharge (Rl) – no satellite-based product is available for the discharge with sufficient temporal extent and only a few rivers are still monitored by public or private networks for the Global Runoff Data Centre (GRDC) that collects discharge data at the global scale. The two discharge datasets used in the following are described in . Groundwater discharge is neglected and considered as an uncertainty source.

The CEFREM-V2 dataset of coastal annual discharge into the Mediterranean Sea is based on in situ observations and some indirect estimates using the Pike formula (Pike1964). In addition, the Laboratoire de Météorologie Dynamique developed the Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE, Polcher et al.1998; Ducharne et al.2003) land surface model chosen here to describe the monthly dynamics of the discharge. Two coastal discharges are available from its routing scheme with two different precipitation forcings: GPCC and Climatic Research Unit (CRU) products. We therefore projected the monthly dynamical patterns from ORCHIDEE towards the CEFREM grid. We then scaled the monthly values of ORCHIDEE to match the CEFREM annual values. For comparison purposes, CEFREM total freshwater inflow into the Mediterranean (without the Black Sea) is 400 km3 yr−1, while ORCHIDEE gives a value of 380 km3 yr−1. The scaling is then a simple way to take into account the anthropogenic impact that is not modeled at the annual scale and the 0.5 resolution. The final product then has the spatial resolution and the annual cumulative value of CEFREM, but with the monthly dynamics of the ORCHIDEE model.

Precipitable water change (ΔW) – we considered two datasets for the precipitable water: the ESA Globvapor dataset and the 6 h reanalysis product from the ECMWF reanalyses (ERA-I, Dee et al.2011). The ERA-I reanalysis has been considered here because precipitable water, although model-based, is largely constrained by satellite observations. In order to compute changes in precipitable water, we also applied the derivative filter: [5∕24 3∕8 $-\mathrm{3}/\mathrm{8}$ $-\mathrm{5}/\mathrm{24}$].

Moisture divergence (Div) – due to the limited temporal extent of the satellite-based data, we used the 6-hourly ERA-I reanalysis product . Among the various re-analyses, ERA-I was chosen here in view of previous results demonstrating advantages in the representation of long term wind variability which plays a key role in the representation of moisture divergence. Nevertheless, have shown the limitation of the reanalysis that do not catch moisture divergence events shorter than at the 6 h. This limitation must be taken into account when closing the WC.

Gibraltar netflow (Gib) – the only multiannual estimate of the Gibraltar netflow based on observations is presented in . A monthly reconstruction of the net transport is used with the effects of the atmospheric pressure removed. This is done for consistency with the oceanic water storage from GRACE. The reconstruction technique used to generate that estimate has proven to be effective to simulate the variability but the uncertainties in the mean value stay large. In , an expert-based assessment of the mean transport is presented. Therefore, in this work we substituted the 2004–2016 mean value of the estimate by the estimate proposed by .

The Mediterranean Sea is also connected to the Red Sea with the Suez channel and to the Black Sea with the Bosporus Strait. The netflow at the Suez channel is neglected . Since no in situ reference is available on the Bosporus netflow, the current work gathers the Mediterranean and the Black Sea into a single reservoir for the integration process. An a posteriori estimate of the Bosporus netflow will be given in Sect. 5 using the water budget equations and the integrated estimates for the other water components.

Table 2Comparison of the uncertainty specifications for terrestrial water components. The weights associated with a variable (computed as the ratio between the particular variable uncertainty with respect to the sum of all the uncertainties ${\mathit{\sigma }}_{i}^{\mathrm{2}}/\sum _{i}{\mathit{\sigma }}^{\mathrm{2}}$) are expressed in percentage.

## 2.3 Validation datasets

The ENSEMBLES observation dataset (EOBS) – in order to validate the precipitation, an additional dataset is used: the EOBS dataset developed by EU-FP6 project ENSEMBLES . It was a regional, well-documented and validated in situ gridded daily dataset at the 0.25 spatial resolution, covering the 1950–2007 period.

FLUXNET – ground-based FLUXNET data were used to validate the evapotranspiration and precipitation over several sites in Europe1. These flux measurements were based on a eddy covariance technique. All stations available in Europe for the 2004–2009 period have been selected. In order to avoid coastal contamination, the three seaside towers IT-Ro2, IT-Noe and ES-Amo have been suppressed.

Total and thermosteric sea-level databases – to validate the sea-level output from the integration technique, we use an independent estimate of the Mediterranean water content. The water content can be estimated as total sea level minus the thermosteric variations (i.e., changes in sea level due to thermal expansion/contraction) . Total sea level is obtained from the Ssalto/Duacs altimeter data that are produced and distributed by the Copernicus Marine and Environment Monitoring Service2. The thermosteric sea-level variations are estimated using two ocean regional reanalyses (MEDRYS, Hamon et al.2016; Bahurel et al.2012, MyOcean,) and two global products that include the Mediterranean Sea (the Met Office Hadley Centre EN-v4 Good et al.2013; Ishii et al.2003, ISHII).

## 2.4 EO uncertainty assumptions

Some studies aimed to characterize the uncertainty of satellite-retrieved products: estimating the relative uncertainty of numerous datasets by the distance to the average product or using non-satellite datasets . Nevertheless, such characterizations are generally product- and site-specific, and for some products used in this work, no uncertainty characterization can be found in the literature. For these reasons we considered the same uncertainty as in .

Table 2 summarizes the uncertainty used in the various integration techniques. The uncertainty is associated with a weight which is the ratio of the sum of all the uncertainties in the WC equation and the uncertainty of the considered variable (computed as ${\mathit{\sigma }}_{i}^{\mathrm{2}}/\sum _{i}{\mathit{\sigma }}^{\mathrm{2}}$ and expressed in percentage). Note that uncertainties in Table 2 stand for the merged product and not for a particular satellite dataset (see Eq. 6). The uncertainties are prescribed by the literature, but they are slightly modified from to handle the special case of the Mediterranean region. used uncertainty values of 10 mm month−1 for each of the four P products and the three E products (leading to 5 and 5.8 mm month−1 for the merged P and E estimate), 5 mm month−1 for each of the three ΔS products (leading to 2.9 mm month−1 for the merged product) and 1 mm month−1 for only one R. The choice of these values was motivated by results of the studies cited in Sect. 1. In order to be closer to , on the one hand, we decide to reduce P uncertainty to 4 mm month−1. This is justified since the de-biasing was done toward the gauge-calibrated TMPA dataset (see , for details). On the other hand, we increased E uncertainty up to 6 mm month−1. The uncertainty of the merged ΔS is estimated to be broadly the same since it is mainly driven by the large pixel resolution of GRACE. Finally, the uncertainty of the discharge R has been increased since the product is partially based on model simulations and the groundwater discharge is not included in the analysis (see Sect. 2). For the atmospheric variables, we consider an uncertainty proportional to the range of variability for the precipitable water change: 1 mm month−1. Following the suggestion from , the reanalysis moisture divergence uncertainty has been set to 6 mm month−1 due to its large range of variability and timescale.

3 EO integration methodologies

## 3.1 Closing the water cycle budget

In this section, the notations are introduced, but additional details can be found in . The WC can be described by the following time-varying budget equations:

$\begin{array}{ll}& \frac{\mathit{\delta }{S}_{\mathrm{l}}}{\mathit{\delta }t}={P}_{\mathrm{l}}-{E}_{\mathrm{l}}-{R}_{\mathrm{l}}\phantom{\rule{0.25em}{0ex}}\left(\mathrm{terrestrial}\right),\\ & \frac{\mathit{\delta }{S}_{\mathrm{o}}}{\mathit{\delta }t}={P}_{\mathrm{o}}-{E}_{\mathrm{o}}+{R}_{\mathrm{l}}^{*}-\mathrm{Gib}\phantom{\rule{0.25em}{0ex}}\left(\mathrm{oceanic}\right),\\ \text{(1)}& & \frac{\mathit{\delta }W}{\mathit{\delta }t}={E}_{\mathrm{l}/\mathrm{o}}-{P}_{\mathrm{l}/\mathrm{o}}-\mathrm{Div}\phantom{\rule{0.25em}{0ex}}\left(\mathrm{atmospheric}\right),\end{array}$

where “l” stands for land and “o” for ocean. If all the components in Eq. (1) are expressed in mm month−1 (area-normalized), then a fourth equality is defined: ${R}_{\mathrm{l}}^{*}=\frac{{\mathbf{A}}_{\mathrm{land}}}{{\mathbf{A}}_{\mathrm{Sea}}}\cdot {R}_{\mathrm{l}}$ for total freshwater input/output with Aland is the total drainage area of the Mediterranean (including the Black Sea), and ASea is the total area of both seas.

We first consider the six terrestrial water components ${\mathbit{X}}_{\mathrm{l}}^{\mathrm{t}}=\left({P}_{\mathrm{l}}$, El, Rl, ΔSl, ΔWl, Divl) and the six oceanic water components ${\mathbit{X}}_{\mathrm{o}}^{\mathrm{t}}=\left({P}_{\mathrm{o}}$, Eo, ΔSo, ΔWo, Divo, Gib). We then define ${\mathbit{X}}_{\mathrm{lo}}^{\mathrm{t}}=\left[{\mathbit{X}}_{\mathrm{l}}$, Xo]t. The closure of the water budget can be relaxed using a centered Gaussian random variable r and ${\mathbit{X}}^{\mathrm{t}}\cdot {\mathbf{G}}_{\mathrm{lo}}^{\mathrm{t}}=\mathbit{r}$, with r∼𝒩(O, ) where

$\begin{array}{ll}& {\mathbf{G}}_{\mathrm{lo}}=\\ \text{(2)}& & \left[\begin{array}{cccccccccccc}\mathrm{1}& -\mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{{\mathbf{A}}_{\mathrm{land}}}{{\mathbf{A}}_{\mathrm{Sea}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& -\mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}\end{array}\right]\end{array}$

which is equivalent to the water budget in Eq. (1) and

$\mathrm{\Sigma }=\left(\begin{array}{cc}{\mathbit{\sigma }}_{\mathrm{l}}^{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& {\mathbit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}\end{array}\right),$

with ${\mathbit{\sigma }}_{\mathrm{l}}=\left(\begin{array}{cc}\mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}\end{array}\right)$ representing the standard deviation of the constrained terrestrial and atmospheric water budget residual over land, and ${\mathit{\sigma }}_{\mathrm{o}}=\left(\begin{array}{cc}\mathrm{2}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}\end{array}\right)$ representing the standard deviations of the constrained oceanic and atmospheric water budget residual over sea. Σ assumes no correlation in the imbalance of the three WCs at monthly and annual scales, at sub-basin scales or entire basin scales.

Let

$\begin{array}{}\text{(3)}& {\mathbit{Y}}_{\mathrm{l}}^{\mathrm{t}}=\left(\begin{array}{lll}\left({P}_{\mathrm{1}},\mathrm{\dots },{P}_{\mathrm{p}},& {E}_{\mathrm{1}},\mathrm{\dots },{E}_{q},& {R}_{\mathrm{1}},\mathrm{\dots },{R}_{r},\\ \mathrm{\Delta }{S}_{\mathrm{1}},\mathrm{\dots },\mathrm{\Delta }{S}_{\mathrm{s}},& \mathrm{\Delta }{W}_{\mathrm{1}},\mathrm{\dots },\mathrm{\Delta }{W}_{v},& {\mathrm{Div}}_{\mathrm{1}},\mathrm{\dots },{\mathrm{Div}}_{\mathrm{d}}\right)\end{array}\right)\end{array}$

be the vector of dimension ${n}_{\mathrm{l}}=p+q+m+s+v+d$ gathering the multiple observations available for each water component over land (similarly Yo of dimension no is defined over sea):

• (P1, P2, …, Pp), the “p” precipitation estimates;

• (E1, E2, …, Eq), the q sources of information for evapotranspiration;

• (R1, R2, …, Rm), the m discharge estimates;

• (ΔS1, ΔS2, …, ΔSs), the “s” sources of information for the water storage change;

• (ΔW1, ΔW2, …, ΔWv), the v precipitable water change estimates;

• (Div1, Div2, …, Divd), the “d” moisture divergence.

The aim of this approach is to obtain a linear filter Kan used to obtain an estimate Xan (“an” stands for analysis) of Xlo based on the observations Ylo:

$\begin{array}{}\text{(4)}& {\mathbit{X}}_{\mathrm{an}}={\mathbf{K}}_{\mathrm{an}}\cdot {\mathbit{Y}}_{\mathrm{lo}}\phantom{\rule{0.25em}{0ex}}\mathrm{with}\phantom{\rule{0.25em}{0ex}}{\mathbit{Y}}_{\mathrm{lo}}=\left[{\mathbit{Y}}_{\mathrm{l}},{\mathbit{Y}}_{\mathrm{o}}\right],\end{array}$

where Kan is a $\mathrm{12}×\left({n}_{\mathrm{l}}+{n}_{\mathrm{o}}\right)$ matrix.

## 3.2 Simple weighting (SW)

A general approach to deal with EO datasets in the analysis of the WC is to choose the best individual dataset for each water component. This is the approach developed, for example, in the GEWEX project. In , an optimal selection (OS) was based on the minimization of the water budget residuals to select the best combination of individual datasets. Using the OS principle facilitates finding datasets coherent to each other and with independent errors . But this kind of strategy limits the use of several sources of information to reduce the uncertainties.

On the other hand, the SW approach benefits from the multiplicity of the observations. EO products, and more generally any estimation of a variable via observations, present two types of errors. (1) Systematic errors related, for instance, to the absolute calibration of the sensor. These can be represented by a bias and/or a scaling factor. (2) Random errors related to retrieval algorithm uncertainties or to missing or inaccurate auxiliary information (e.g., cloud mask) or to the sensor itself. These are often characterized by a standard deviation using a Gaussian hypothesis. From a statistical point of view, using the average of several estimates reduces the random errors of the estimation if no bias errors are present in the estimates. The merging process such as in Eq. (4) requires unbiased estimates (Aires2014). The difficulty is that, as for uncertainties (Sect. 2.4), it is rather difficult to obtain bias estimates from the literature for every dataset used in this approach. A pragmatic strategy is to set the reference as the mean state for each component. Then, all the sources of information for this component are bias-corrected toward this reference . A slightly modified version of the bias correction is to choose one reference among the datasets and apply the bias correction. We opted for the modified version and de-biased the EO using the climatological season of the TMPA product as a reference . Therefore, our SW methodology is first based on a seasonal bias correction to reduce the systematic biases and is then followed by a weighted average of the corrected estimates to reduce the random errors. After the seasonal de-biasing, all the precipitation products will have a similar seasonality, but their inter-annual trend or monthly variations will still be different. In particular, the seasonal de-biasing will not change the trend of the EO products.

The SW methodology uses the diversity of datasets to reduce the random errors. Let us consider the p precipitation observations Pi associated with Gaussian errors ϵi∼𝒩(O, σi). σi is the standard deviation of the ith estimate. The SW precipitation estimate PSW is given by the weighted average:

$\begin{array}{}\text{(5)}& {P}_{\mathrm{SW}}=\frac{\mathrm{1}}{p-\mathrm{1}}\sum _{i=\mathrm{1}}^{p}\frac{\sum _{k\ne i}{\left({\mathit{\sigma }}_{k}\right)}^{\mathrm{2}}}{\sum _{k}{\left({\mathit{\sigma }}_{k}\right)}^{\mathrm{2}}}{P}_{i}.\end{array}$

This equation is valid when there is no bias error in the Pis (thanks to the preliminary bias correction) and is optimal when the errors ϵi are statistically independent of each other. This expression is valid for the other water components. The variance of the PSW uncertainties is then given by

$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{{P}_{\mathrm{SW}}}=\frac{\mathrm{1}}{\left(p-\mathrm{1}{\right)}^{\mathrm{2}}}\sum _{i=\mathrm{1}}^{p}{\left(\frac{\sum _{k\ne i}{\left({\mathit{\sigma }}_{k}\right)}^{\mathrm{2}}}{\sum _{k}{\left({\mathit{\sigma }}_{k}\right)}^{\mathrm{2}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{i}^{\mathrm{2}}.\end{array}$

This is important information because it gives the uncertainty of the estimates in Eq. (5). It shows that the PSW errors can be significantly reduced by increasing the number p of observations.

Following Eq. (5) the SW state vector XSW can be defined as

$\begin{array}{}\text{(7)}& {\mathbit{X}}_{\mathrm{SW}}={\mathbf{K}}_{\mathrm{SW}}\cdot {\mathbit{Y}}_{\mathrm{lo}},\end{array}$

where KSW is a $\mathrm{12}×\left({n}_{\mathrm{l}}+{n}_{\mathrm{o}}\right)$ matrix in which each line represents Eq. (5) for 1 of the 12 water components (the first one for the precipitation estimate, the second for the evapotranspiration, etc.) and based on the (nl+no) observations. Since no specific uncertainty specifications were available in the literature for the Mediterranean basin, the uncertainties are assumed to have the same standard deviation σi in the following.

## 3.3 Post-filtering (PF)

In the SW approach, each water component is weighted (see Eqs. 67) based on its a priori uncertainty (Sect. 2), but no closure constraint is imposed on the solution XSW. Several methods were considered in to introduce a WC budget closure constraint on the SW solution. However, Monte Carlo simulations have shown that the SW solution associated with a so-called post-filtering (PF) provides results as good as more complex techniques such as variational assimilation.

The PF approach has been introduced to impose the closure constraint on a previously obtained solution. Here we use XSW as the first guess on the state vector Xlo. In , the PF was used and tested without any model, as a simple post-processing step after the SW. Following , the current study implements the PF filter with a relaxed closure constraint characterized by its uncertainty covariance :

$\begin{array}{}\text{(8)}& {\mathbit{X}}_{\mathrm{PF}}=\left(I-{\mathbf{K}}_{\mathrm{PF}}\cdot {\mathbf{G}}_{\mathrm{lo}}{\mathrm{\Sigma }}^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{lo}}^{\mathrm{t}}\right)\cdot {\mathbit{X}}_{\mathrm{SW}},\end{array}$

where ${K}_{\mathrm{PF}}=\left({\mathbf{B}}_{\mathrm{lo}}^{-\mathrm{1}}+{\mathbf{G}}_{\mathrm{lo}}\sum ^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{lo}}^{\mathrm{t}}{\right)}^{-\mathrm{1}}$ and Blo is the error covariance matrix of the first estimate on Xlo.

We can make XSW explicit to obtain the linear operator Kan of Eq. (4):

$\begin{array}{}\text{(9)}& {\mathbit{X}}_{\mathrm{an}}={\mathbit{X}}_{\mathrm{PF}}=\left(I-{\mathbf{K}}_{\mathrm{PF}}\cdot {\mathbf{G}}_{\mathrm{lo}}{\mathrm{\Sigma }}^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{lo}}^{\mathrm{t}}\right)\cdot {K}_{\mathrm{SW}}\cdot {\mathbit{Y}}_{\mathrm{lo}},\end{array}$

so that ${K}_{\mathrm{an}}=\left(\mathbf{I}-{\mathbf{K}}_{\mathrm{PF}}\cdot {\mathbf{G}}_{\mathrm{lo}}{\mathrm{\Sigma }}^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{lo}}^{\mathrm{t}}\right)\cdot {K}_{\mathrm{SW}}$. The PF step (budget closure) consists in partitioning the budget residual among the 12 components at each time step, independently. This technique allows a satisfactory WC budget closure for each basin to be obtained. The residual term r could be reduced in the SW + PF approach by decreasing the variance Σ in Eq. (8). Nevertheless, an excessively tight closure constraint is in contradiction with the large inherent uncertainties in original observations.

Following we enforced the budget closure by frequency range to avoid high-frequency errors impacting the low-frequency variables such as the evapotranspiration (mainly driven by annual vegetation growth; ). We first decomposed each parameter into high- and low-frequency components considering a cut-off frequency of 6 months (using a FFT decomposition). The budget is then applied independently at low and high frequencies. The high-frequency component of E is then not included in the high budget closure. The linearity of PF and FFT ensures the budget closure of the re-composed final product. In the following temporal multi-scaling, the annual constraint is applied only to the low-frequency budget closure.

Spatial multi-scaling – it is possible to impose a WC budget closure simultaneously over the six sub-basins, the full basin and the ocean (i.e., Mediterranean and Black Sea). Let us consider the total WC state vector

$\begin{array}{}\text{(10)}& {\mathbit{X}}^{\mathrm{t}}={\left[{\mathbit{X}}_{\mathrm{l}}^{\left(\mathrm{1}\right)},{\mathbit{X}}_{\mathrm{l}}^{\left(\mathrm{2}\right)},{\mathbit{X}}_{\mathrm{l}}^{\left(\mathrm{3}\right)},{\mathbit{X}}_{\mathrm{l}}^{\left(\mathrm{4}\right)},{\mathbit{X}}_{\mathrm{l}}^{\left(\mathrm{5}\right)},{\mathbit{X}}_{\mathrm{l}}^{\left(\mathrm{6}\right)},{\mathbit{X}}_{\mathrm{o}}\right]}^{\mathrm{t}}\end{array}$

that includes the six water components ${\mathbit{X}}_{\mathrm{l}}^{i}$ over each sub-basin i of area ${\mathbf{A}}_{\mathrm{l}}^{\left(i\right)}$ and ocean. The “closure” matrix becomes

$\begin{array}{}\text{(11)}& {\mathbf{G}}_{\mathrm{lo}}=\left(\begin{array}{ccccc}{\mathbf{G}}_{\mathrm{l}}^{\left(\mathrm{1}\right)}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathbf{G}}_{\mathrm{l}}^{\left(\mathrm{2}\right)}& \mathrm{\cdots }& \mathrm{0}& \mathrm{0}\\ \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& {\mathbf{G}}_{\mathrm{l}}^{\left(\mathrm{6}\right)}& \mathrm{0}\\ {\mathbf{L}}_{\mathrm{lo}}^{\left(\mathrm{1}\right)}& {\mathbf{L}}_{\mathrm{lo}}^{\left(\mathrm{2}\right)}& \mathrm{\cdots }& {\mathbf{L}}_{\mathrm{lo}}^{\left(\mathrm{6}\right)}& {\mathbf{G}}_{\mathrm{o}}\end{array}\right)\end{array}$

with

$\begin{array}{ll}& {\mathbf{G}}_{\mathrm{l}}^{\left(i\right)}=\left[\begin{array}{cccccc}\mathrm{1}& -\mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& -\mathrm{1}\end{array}\right],\\ & {\mathbf{L}}_{\mathrm{lo}}^{\left(i\right)}=\left[\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& \frac{{\mathbf{A}}_{\mathrm{l}}^{\left(i\right)}}{{\mathbf{A}}_{\mathrm{Sea}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right],\\ \text{(12)}& & {\mathbf{G}}_{\mathrm{o}}=\left[\begin{array}{cccccc}\mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ -\mathrm{1}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& -\mathrm{1}& \mathrm{0}\end{array}\right].\end{array}$

The last row of Glo represents the overall budget closure, including all the sub-basins and the ocean. The dimensions of the covariance matrices Blo and  are increased following the new size of the state vector Xlo. No cross terms in Blo and  are included, meaning that there is no dependency of the first guess and closure errors among the sub-basins.

Temporal multi-scaling – it is also possible to impose a WC budget closure simultaneously at monthly and annual scales. With monthly closure, the annual closure should automatically be obtained, but due to the relaxation of the closure constraint, the annual closure would be relaxed too. We control here the yearly closure constraint with an uncertainty of 1 mm. Furthermore, we impose a yearly closure assuming no groundwater storage change at the annual scale over land (representing an additional constraint on ΔSl to ensure that no bias is introduced for this variable during the PF process). In this framework, monthly closures are now interdependent in the given year and the new state vector is

$\begin{array}{}\text{(13)}& {\mathbit{X}}_{\mathrm{year}}^{\mathrm{t}}={\left[{\mathbit{X}}^{\mathrm{Jan}},\mathrm{\cdots },{\mathbit{X}}^{\mathrm{Dec}}\right]}^{\mathrm{t}},\end{array}$

with Xm the total state vector X defined in Eq. (10), for month m. The closure is applied independently for the 4 years of the 2004–2009 period, but the 12 months of each year are closed independently.

The closure matrix GAlo that includes closure for the 12 months of the year and the full year is derived from the monthly constraint of Eq. (11) and defined as

$\begin{array}{}\text{(14)}& {\mathbf{GA}}_{\mathrm{lo}}=\left(\begin{array}{cccc}{\mathbf{G}}_{\mathrm{lo}}& \mathrm{0}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{0}& {\mathbf{G}}_{\mathrm{lo}}& \mathrm{\cdots }& \mathrm{0}\\ \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }\\ \mathrm{0}& \mathrm{0}& \mathrm{\cdots }& {\mathbf{G}}_{\mathrm{lo}}\\ {\mathbf{N}}_{\mathrm{lo}}& {\mathbf{N}}_{\mathrm{lo}}& \mathrm{\cdots }& {\mathbf{N}}_{\mathrm{lo}}\end{array}\right),\end{array}$

where Nlo is the modified closure matrix Glo in which the matrix ${\mathbf{G}}_{\mathrm{l}}^{\left(i\right)}$ is rewritten in ${\mathbf{N}}_{\mathrm{l}}^{\left(i\right)}$ by imposing ΔSl=0:

$\begin{array}{}\text{(15)}& {\mathbf{N}}_{\mathrm{l}}^{\left(i\right)}=\left[\begin{array}{cccccc}\mathrm{1}& -\mathrm{1}& -\mathrm{1}& \mathbf{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& -\mathrm{1}\end{array}\right].\end{array}$

The last row of GAlo represents the annual budget closure considering no storage change at the annual scale over land, including all the sub-basins and the ocean. The dimensions of the covariance matrices Blo and  are increased once again following the new size of the state vector Xyear. No cross terms in Blo and  are included, meaning that there is no dependency of the first guess and closure errors between the months.

This SW + PF technique is able to deal only with time series (average value over the considered sub-basins), not with maps (pixel), since the discharge is not available at this resolution. Therefore, in order to obtain a multi-component dataset that closes the WC budget and has spatial patterns at the pixel level, another technique needs to be used.

## 3.4 INTegration (INT)

The INT methodology allows one to extrapolate the results obtained with the previous SW + PF, from the sub-basin to pixel scales. To obtain a pixel-wise closure, assimilate satellite data into the VIC model at the pixel scale (0.5) using the VIC pixel water storage and runoff information. extrapolated at the global scale the results of the WC closure of several large river basins around the globe, by using surface classes that intend to discriminate between EO error types, preserving as much as possible the hydrological coherency.

The INT approach proposed here uses the WC closure over the Mediterranean sub-basins to extrapolate the closure correction to the surrounding area. The methodology is presented in its various steps in Fig. 2 for precipitation and evaporation, for a particular month. In this analysis, we consider only the Mediterranean sub-basins and their close neighborhood, so a simple spatial interpolation of the closure correction is assumed to be sufficient.

The SW + PF method (Fig. 2, second row) provides a WC budget closure over the six sub-basins, for each month $m=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{\cdots },\phantom{\rule{0.125em}{0ex}}\mathrm{72}$ of the 2004–2009 period.

Figure 2Steps of the spatialization of the budget closure for the INT solution, from the SW to the INT solutions: precipitation (a) and evapotranspiration (b) for July 2008. Units are in mm month−1.

The INT method requires a scaling factor to go from the SW to the SW + PF solution at the sub-basin scale. We define ${\mathit{\beta }}^{\left(i\right)}\left(m\right)={P}_{\mathrm{PF}}^{\left(i\right)}\left(m\right)/{P}_{\mathrm{SW}}^{\left(i\right)}\left(m\right)$ (for precipitation here), the ratio between the SW and the SW + PF solution, for each sub-basin i and month m. This ratio can be used to scale the SW dataset towards the SW + PF solution at the basin scale, for a particular month m, in the following way:

$\begin{array}{}\text{(16)}& {P}_{\mathrm{INT}}^{\left(i\right)}\left(m\right)={\mathit{\beta }}^{\left(i\right)}\left(m\right)\cdot {P}_{\mathrm{SW}}^{\left(i\right)}\left(m\right)\left(={P}_{\mathrm{PF}}^{\left(i\right)}\left(m\right)\right).\end{array}$

For water storage change or moisture divergence, this β could become negative. In this case, the bias correction ${\mathit{\gamma }}^{\left(i\right)}\left(m\right)={P}_{\mathrm{PF}}^{\left(i\right)}\left(m\right)-{P}_{\mathrm{SW}}^{\left(i\right)}\left(m\right)$ is used instead:

$\begin{array}{}\text{(17)}& \mathrm{\Delta }{S}_{\mathrm{INT}}^{\left(i\right)}\left(m\right)=\mathrm{\Delta }{S}_{\mathrm{SW}}^{\left(i\right)}\left(m\right)+{\mathit{\gamma }}^{\left(i\right)}\left(m\right)\left(=\mathrm{\Delta }{S}_{\mathrm{PF}}^{\left(i\right)}\left(m\right)\right).\end{array}$

The β scaling is defined at the sub-basin scale, but if interpolated spatially, it could be used at the pixel scale to obtain a truly spatialized solution.

Let us define a scaling map at the pixel level α such that for each pixel j in sub-basin i, for each month m, α(j, m)=β(i)(m) (or γ(i)(m)). When used as it is, the convolution of SW and α maps allows for the spatialization of the sub-basins' closure (Fig. 2, third row) with

$\begin{array}{}\text{(18)}& \underset{j\in {\mathbf{A}}_{\mathrm{l}}^{\left(i\right)}}{\iint }{P}_{\mathrm{SW}}\left(j,m\right)\stackrel{\mathrm{˙}}{\mathit{\alpha }}\left(j,m\right)={\mathit{\beta }}^{\left(i\right)}\left(m\right)\cdot {P}_{\mathrm{SW}}^{\left(i\right)}\left(m\right)={P}_{\mathrm{INT}}^{\left(i\right)}\left(m\right).\end{array}$

However, this product presents not only a discontinuity across the sub-basins (where different scaling factors β are defined) but also no value can be provided outside of the sub-basins.

To solve these two issues, the α scaling maps are interpolated/extrapolated.

• Interpolation – a region of 200 km on either side of the frontier between two sub-basins i1 and i2 is defined, and a smooth interpolation is performed between the two scaling factors ${\mathit{\beta }}^{\left({i}_{\mathrm{1}}\right)}\left(m\right)$ and ${\mathit{\beta }}^{\left({i}_{\mathrm{2}}\right)}\left(m\right)$ based on the distance to the frontier. This interpolation of the scaling factors α between two sub-basins can introduce errors (closure residuals can slightly increase), but it will be shown that this effect is limited and that the bottom equations (in parenthesis) in Eqs. (16)–(17) stand overall.

• Extrapolation – an extrapolation of the α maps is then performed to have a scaling factor α outside of the sub-basin domain. This extrapolation is weighted according to the respective distances to the two closest sub-basins.

The INT product is the convolution between the SW dataset with the resulting scaling map α that constrains the WC budget closure. INT is then an optimized version of SW in which the WC budget closure correction has been extended at the pixel scale. The fourth row of Fig. 2 shows the resulting INT product and its spatial coverage. The continuity issues between the sub-basins have been solved, and the extrapolation allows for a spatial coverage over the entire domain.

The extrapolation of a closure constraint is interesting at the technical level because for other regions, or when working at the global scale, some form of interpolation/extrapolation between the monitored sub-basins is required . The extrapolation outside of the Mediterranean region will also allow for the use of more in situ observations for the evaluation; this will help the testing of the generalization ability of the extrapolation scheme. The justification of this interpolation/extrapolation is based on the assumption that most of the WC imbalance is due to satellite errors (this assumption is used for the CAL methodology too). The closure constraint is supposed to improve the satellite estimate by reducing the bias and random errors. If no other information is used (such as surface type; see ), the EO errors can be considered spatially continuous, and it then makes sense to extrapolate results based on this spatial continuity. We perform the main analysis over the Mediterranean basin and test the extrapolation scheme over well-monitored locations.

The difference between the SW and INT estimates, represented in the last row of Fig. 2, is then directly related to the pixel-wise interpolated scaling factor α. Discontinuity between the sub-basins is smoothed. The north of Europe excluding France is mainly driven by the scaling factor on the BLS region. That is consistent with the updated Köppen climate classification . Since the SW + PF solution is available over the 2004–2009 period only, INT can be obtained only over this period.

Table 3Main characteristics of the five merging methods in this study: EO stands for Earth observation satellite datasets, and * means not considering the GRACE period coverage. The last column shows the capability of the methodology to close the WC budget. “ ” means bad closure, “” means quite bad closure, “+” means quite good closure and “+ +” means good closure.

## 3.5 CALibration (CAL)

To obtain the INT solution, many EO datasets were combined: multiple datasets for each water component (the SW part), and for the various WC components (the PF part). However, if one of the datasets is missing, the INT solution cannot be estimated and this will result in a gap in the time record, and shorter time series of the integrated database.

In , a “Closure Correction Model” (CCM) was introduced to correct each dataset independently, based on the results of the SW + PF integration. The CCM is defined as a simple linear transformation with a scaling factor a and an offset b, such that ${\mathbit{X}}_{\mathrm{cal}}=a\cdot {\mathbit{X}}_{\mathrm{obs}}+b$. The CCM parameters a and b were calibrated by computing a linear regression between the original observation datasets and the SW + PF components.

A similar approach can be used, with the INT solution as a reference instead of the SW + PF. Instead of calibrating the original EO datasets using basin-scale data, we propose here to calibrate the SW solution towards the INT solution at the pixel scale. This calibration of the SW allows a long-term dataset at the pixel scale to be obtained (see Table 3) with WC budget closure statistics closer to the INT solution. In our tests (not shown), the linear regression is quite satisfactory for the calibration, and it is not necessary to use a more complex statistical regression tool such as a neural network.

The merging/integration techniques used in this study are described in Table 3.

Table 4rms of the WC budget residual (in mm month−1) over the sub-basin using OS, SW, SW + PF, INT and CAL solutions and for the period 2004–2009. Percentage of improvement of the rms of the residuals from the SW solution to constrained methods is also shown. For comparison purposes, the result using the ERA-I dataset is also depicted.

4 Evaluation of the integrated datasets

In this section, the obtained integrated datasets are first evaluated in terms of WC budget closure. Our EO dataset integration technique is based on the closure of the WC budget. This is a physical constraint, but in some cases (e.g., missing important water component), this constraint could result in a degraded estimation of the components. Therefore, available in situ data (precipitation, evapotranspiration and seawater level) are used to validate some of the water components of the integrated dataset. This evaluation is performed at two different spatial scales: the sub-basin scale and the pixel scale.

## 4.1 Water cycle budget closure

The residuals of the surface and atmospheric WC budgets for the Mediterranean region are computed at the monthly scale, over the 2004–2009 period. The root mean square (rms) statistics of these residuals are summarized in Table 4 for the six considered products (ERA-I, OS, SW, SW + PF, INT and CAL). Percentages of improvement of the rms of the residuals with respect to the SW solution are also shown for comparison purposes.

ERA-I provides the reanalysis product for all the variables except for the water storage and the discharge, to keep the comparison consistent. It should be noted that ERA-I does not have any water conservation constraint. The optimal selection is given by TMPA precipitation; GLEAM evapotranspiration and OAFlux evaporation; GRGS water storage change over land and JPL water storage change over sea; GPCC-forced ORCHIDEE-CEFREM discharge; and the derivate Globvapor for atmospheric water vapor change. Only one dataset is available for the moisture divergence . As shown in , and , the SW merging procedure reduces the WC budget residuals at the sub-basin scale, by reducing the random errors of the EO data. The SW product outperforms the ERA-I reanalysis and the OS product. However, the full closure is generally still not satisfactory with this technique. The SW + PF procedure closes the water budget over all the sub-basins, and over the surface and in the atmosphere, with an rms of the residual of about 4 mm month−1. The surface budget residuals are drastically reduced: from 72 % over the GR-TR-IL sub-basin and up to 94 % for the Mediterranean Sea. This shows the necessity to use a WC budget closure constraint that links the six water components.

The INT product provides satisfactory budget closure results (from 61 % to 94 %), even if they are slightly degraded compared to the SW + PF (due to the interpolation process between sub-basins). Since no interpolation has been applied over the Mediterranean Sea, the statistics are equal to the SW + PF.

The CAL product improves the WC budget residuals less compared to INT. Nevertheless, the rms of the residuals for these products is reduced over all sub-basins compared to the SW solution.

Figure A1 gives, in the Appendix, the 2004–2009 time series of all the water component estimates for the various methodologies (SW, SW + PF, INT and CAL) over the various sub-basins as well as the probability density function of the residuals. This figure shows how the WC closure impacts the time series.

## 4.2 Evaluation at the sub-basin scale

Since the WC budget closure constraint was imposed at the sub-basin scale (see Sect. 3), the evaluation of the integrated product is done at this scale too. Two metrics are used here, the rms of the difference (rmsd) with in situ measurements and the CORRelation (CORR). Only multiple-EO integrated datasets are compared in the two following sections.

Terrestrial precipitation – Table 5 provides the comparison of the EOBS gridded gauge precipitation dataset (Sect. 2.2) with the SW, SW + PF, INT and CAL solutions, in terms of temporal correlation (at the monthly and sub-basin scales), and rmsd, for each sub-basin and for the continental scale (land included in Fig. 1). Since the SW + PF product is defined only on the Mediterranean drainage sub-basins, no statistics are shown for this approach over the continental region (last column). For the rmsd error statistics, results are also provided as improvements compared to the SW solution.

Table 5Comparison of the SW, SW + PF, INT and CAL precipitation solutions with the EOBS dataset, in terms of correlations, rmsd, and percentage of improvement of the rmsd compared to the SW solution.

Over all the sub-basins, the SW + PF methodology improves results compared to the unconstrained SW method. Even if the correlation of SW with EOBS is already good, the closure constraint improves this correlation to 0.84 (from 0.81) over the MA-DZ-TN sub-basin. This is true even over the complex sub-basins such as Alp-IT-ADR. SW + PF also reduces the rmsd with EOBS (by up to 20 %). These results show the positive effect of the closure constraint on precipitation. Without explicitly constraining satellite precipitation products towards the in situ data, SW + PF statistics are still improved.

The INT product shows similar CORR and rmsd statistics as SW + PF over the Mediterranean sub-basins, with a slight decrease in the CORR with EOBS over the ES-Pyr sub-basin. Over the continental region, INT improves the correlation compared to SW (from 0.78 to 0.80) while reducing by 17 % the rmsd. Therefore, the interpolation process between the sub-basins (see the spatialization in Sect. 3.4) does not degrade the solution inside the sub-basins, while the extrapolation outside of them improves the unconstrained SW statistics over the whole continent. This is a true benefit since INT presents comparable performances with the SW + PF in terms of closure capability and closeness to in situ measurements, with the advantage of the spatial variability at the pixel scale.

Finally, the CAL precipitation product shows results as good as SW (slightly better for the whole continental region) for the CORR, and smaller rmsd with EOBS. The CAL product does not close the WC budget as well as the INT solution, but it has the advantage of being available over a longer time record (1980–2012) compared to the 2004–2009 INT period.

Seawater-level change – the seawater storage (related to the seawater level) change over the Mediterranean Sea (excluding the Black Sea) is tested using altimetry and thermal datasets over the 2004–2009 period. First, the thermal content estimates of the four datasets presented in Sect. 2.2 are merged into one single estimate. The merged thermal content estimate is then subtracted from the AVISO altimetry seawater level. The monthly change is then computed using the same derivative filter as the one used for GRACE: [5∕24 3∕8 $-\mathrm{3}/\mathrm{8}$ $-\mathrm{5}/\mathrm{24}$].

Figure 3 shows the altimetry estimate and the various methodologies estimates. Since the Mediterranean Sea is considered without the Black Sea for this evaluation, there is no SW + PF estimate (that added the Mediterranean Sea and the Black Sea). While the SW solution has a 0.52 CORR and a 12.2 mm month−1 rmsd with respect to the altimetry estimate, INT statistic are 0.58 for the CORR and 11.8 mm month−1 for the rmsd, and CAL 0.56 for the CORR and 11.8 mm month−1 for the rmsd. Here again, the INT estimate outperforms the unconstrained SW methodology in both CORR and rmsd. CAL presents also better results than SW but the CORR with altimetry is slightly reduced compared to INT. No inter/extrapolation has been used in INT for the “Mediterranean Sea plus the Black Sea” sub-basin and the improvement of INT versus SW is due only to the impact of the closure constraint. Nevertheless, the SW + PF approach closes the WC over the Mediterranean Sea within the Black Sea (no information about the Bosporus netflow) and the spatial downscaling in INT is needed to discriminate the closure correction above the two seas. Using the WC closure over the Mediterranean and Black Sea improves the water storage change estimates.

Figure 3Seawater-level evaluation of SW, INT and CAL estimates compared to altimetry minus thermal content. The correlation difference is statistically significant at the 70 % level based on the T-test.

## 4.3 Evaluation at the pixel scale

The INT and CAL estimates are here evaluated at the pixel scale, for precipitation and evapotranspiration. Improvements of SW by INT and CAL are measured using in situ measurements of precipitation and evapotranspiration from the FLUXnet database, available over the Mediterranean region, for the 2004–2009 period (Sect. 2.2).

Figure 4 presents the scatterplots of the rmsd between the SW estimate (ESW) and INT or CAL (Ecor for “corrected”) datasets with the FLUXnet evapotranspiration data (EFLUX), for each station. The 1:1 line is also shown. Each dot under the 1:1 line represents an improvement at the corresponding station from SW solution to INT and/or CAL. INT and CAL improve evapotranspiration estimates for more than 53 % of the sites. The distribution of the differences in the encapsulated figure is slightly narrowed by the INT and CAL compared to the SW solution. Location of the station where the closure improves the rmsd with the flux measurement is shown in green if INT and CAL both improve the estimate, blue when only CAL improves, and magenta when only INT improves. Red dots represent stations where there is a degradation in both INT and CAL. The evaluation of the EO estimate at 0.25 spatial resolution using tower sites should be made with caution. The poor performance of satellite estimates over particularly complex topography (mountainous rainfall) or coastal pixels with land–sea contamination could explain the difference between the INT estimate and the FLUXNET measurement on these particular locations.

Figure 5 presents the scatterplots of the rmsd between the SW estimate (PSW) and the INT or CAL (Pcor) datasets with the FLUXnet precipitation data (PFLUX), for each station. Over most stations (82 %), the INT and CAL solutions improve the precipitation estimate compared to SW. Locations of improved sites are shown with the same color code as in Fig. 4. It can be seen in Fig. 5 that red dots are located mainly in mountainous or coastal regions. These two types of landscape are really challenging for precipitation retrieval due to snow precipitation on the one hand and coastal sea–land contamination on the other.

5 A coherent multi-component dataset for the water cycle monitoring

In this section, the integrated datasets are used to deliver updated estimates of the Mediterranean WC budget. The impact of hydrological constraint (PF) as well as the INTegration (INT) and CALibration (CAL) processes on the spatial averaging of the water component estimates and the WC budget residuals, over the several Mediterranean sub-basins, is summarized in Fig. A1 of the Appendix.

Figure 4Top panel: scatterplot of the rmsd between FLUXnet station and the SW, INT and CAL products, for evapotranspiration. Dots under the 1:1 line (green) show improvement, and dots over the line (red) show degradation. INT and CAL results are superposed for some locations, meaning that the linear approximation in CAL is enough to mimic the INT scaling factors at these location. The encapsulated figure shows the distribution of the differences with the Fluxnet estimates. Bottom panel: location of the FLUXnet stations used for validation: green dots show an improvement for INT and CAL compared to SW (INT+ and CAL+), blue dots show improvement only for CAL (INT and CAL+), and magenta only for INT (INT+ and CAL). Red dots is where no improvement is observed (INT and CAL). Blue line limits the total basin area.

## 5.1 Analysis of the Mediterranean WC

The mean fluxes of the Mediterranean WC and associated variability over the 2004–2009 period are depicted in Fig. 6. The WC is analyzed over its natural sub-basin boundaries. The variability is computed as the standard deviation of the annual values over the period. These values have been computed over the respective terrestrial or oceanic sub-basins, considering all the drainage area in western Europe and BLS or in Africa within Turkish coastal region but without considering the Nile River basin (for which just its discharge is represented), Black Sea or Mediterranean Sea. The large font corresponds to the INT, and the little font is for SW. The two values for the netflow estimate at the Bosporus Strait are estimated as the deficit term of the water budget equation, computed over the Mediterranean and the Black Sea independently. Using INT estimates (i.e., closure of the two seas at once), the two values are in better agreement with each other than with the two SW estimates. In the following, only the INT values are described.

Figure 5Same as Fig. 4 but for precipitation.

Figure 6 shows the uneven water contribution between the European (314±57 km3 yr−1 for the total discharge) and the African (21±30 km3 yr−1) drainage area to the Mediterranean Sea budget. Furthermore, it shows the role of the Black Sea in the global Mediterranean WC. Most of the European freshwater flows to the BLS (398±70 km3 yr−1; it represents more than 50 % of the European discharge), where the EP balance allows for an equal contribution to the Mediterranean Sea budget though the Bosporus Strait input. Considering the Nile discharge, the closure optimization increase the discharge value (from 19±6 to 76±30 km3 yr−1). Recent discussions on the Nile discharge can be found in . Our new discharge estimates include the groundwater discharge passing through the aquifers.

Table 6Comparison in the literature for the Mediterranean Sea (without the Black Sea) average annual mean fluxes and their associated variability (in mm yr−1). While the variability of the real product is computed as the standard deviation of annual values, the uncertainty associated with the regional climate models' mean is the inter-model spread. The periods of analysis for the various studies are recalled.

After closure optimization, the annual precipitation, evapotranspiration and moisture divergence over the European drainage area are estimated to be 770±40, 510±10 and 77±60 km3 yr−1, respectively. Europe accumulates most of the moisture coming from the Mediterranean Sea (1787±200 km3 yr−1), while the Black Sea poorly evacuates its moisture towards land (91±60 km3 yr−1). Over land the contribution of the African part to the global moisture divergence is quite high (274±43 km3 yr−1 mainly due to the presence of the Atlas mountains). The two netflow estimates at the Bosporus Strait are very close, with a difference lower than its associated uncertainty in Fig. 6. Freshwater inputs at the two Mediterranean straits (Bosporus and Gibraltar) compensate for the very large evaporation loss (3372±88 km3 yr−1) of the Mediterranean Sea. This process represents more than twice the precipitation (1499±102 km3 yr−1).

Figure 6Mean annual fluxes (km3 yr−1) of the Mediterranean WC and associated uncertainties in SW (small font) and INT (large font) during the 2004–2009 period.

Figure 6 represents the whole WC over the considered region with its main features: the role of the Mediterranean Sea as the moisture and energy reservoir for the surrounding land; the poor contribution of the African coasts in terms of water resource, and the role of the Black Sea as the buffer process for the freshwater input. This quasi-triangular process emphasizes the hydrological link between the surrounding land and the two seas.

## 5.2 Comparison of the Mediterranean flux estimates with the literature

Table 6 summarizes the comparison of the various estimates of the water fluxes in the current analysis with what can be found in the literature. The various annual mean estimates are based on different time periods and comparison must be taken with caution since some interannual variability is likely to be due to the change in hydrologic regime. focused on the Mediterranean Sea heat and water budget using an ensemble of ERA-40-driven high resolution Regional Climate Models (RCMs) from the FP6-EU ENSEMBLE database. The atmospheric budget was not considered in and no moisture divergence estimate was provided. For comparison purposes, we decided to select the RCM ensemble-mean estimate and two particular models: the Danish HIRHAM and the Canadian CRCM . These two models have been selected since their EP estimates are the extremes of the RCMs ensemble. In , the netflow at Gibraltar was estimated as the deficit term of the WC: $\mathrm{Gib}=E-P-R-\mathrm{Bos}$.

analyzed the WC over the Mediterranean region in the context of the NAO teleconnection over the 1979–1993 period using two reanalyses (ERA-40 and NCEP-NCAR) for precipitation, evaporation and moisture divergence. They used the discharge data from the monitored rivers through the Mediterranean Hydrological Cycle Observing System (MED-HYCOS) and GRDC. Their estimate includes a total Mediterranean input of 100 mm yr−1 from MED-HYCOS and the Bosporus input of 75 mm yr−1 from the literature (Lacombe and Tchernia, 1972). estimated the netflow at Gibraltar as the balance of the Mediterranean water deficit using the equation $\mathrm{Gib}=\mathrm{Div}-R-\mathrm{Bos}$ coming from the oceanic and atmospheric budgets and the null assumptions about the storage change. used old versions of the reanalyses and some remarks have already been made on the precipitation and evapotranspiration estimates for these versions. Nevertheless, to our knowledge, was the last effort to estimate the atmosphere WC over the Mediterranean.

reviewed the state-of-the-art in the quantification of the various water component estimates. Their estimates presented in Table 6 are the best consensual values among the scientific community. They are based on several studies and take into account the results of and for example. In particular, the mean Gibraltar netflow estimate from has been commented on and the new mean is provided in .

Table 6 also shows the results from before and after their satellite data optimization based on a variational assimilation at the annual scale. The constraints of the fluxes over the Mediterranean Sea and the Black Sea were made independently (considering no netflow at the Bosporus Strait). The Mediterranean Sea was closed with no exchange to the Atlantic at Gibraltar (no netflow). provided no explicit discharge for the Mediterranean Sea, but only for the Eurasian continent.

For the four mentioned articles, only the Mediterranean Sea without the Black Sea is considered. No estimate from SW + PF methodology is provided in Table 6. Our integrated dataset is the only one to use direct observations for the Gibraltar netflow and to compute the Bosporus netflow via a WC budget. For all estimates, Table 6 presents their associated variability. While the variability of real products is computed as the standard deviation of the annual values, the variability associated with the RCM mean is the inter-model spread (i.e., proxy of the uncertainty).

Evaporation – the RCM ensemble mean for the annual evaporation is 1254 mm yr−1 with an inter-model spread of 164 mm.yr−1. Some RCMs evaluated higher annual evaporation, as the HIRHAM model that estimated 1377±55 mm yr−1. On the contrary, found a lower evaporation with the reanalyses (1113 and 934 mm yr−1 for NCEP and ERA). estimated much higher evaporation and higher annual variability with an mean annual value of 1391±157 mm yr−1 using only OAFlux and 1420±109 mm yr−1 after optimization. Our unconstrained SW solution gives an annual value of 1300±34 mm yr−1 and our constrained INT product gives 1295±33 mm yr−1. The CAL estimate is close to the INT solution.

Precipitation – the RCM ensemble mean for the annual precipitation was 442±84 mm yr−1 which is quite close to the NCEP reanalyses value in . Satellite estimates in both and the current study indicate higher precipitation: 576 and 571 mm yr−1 in and from 573 to 577 mm yr−1 after the closure constraint in this work. SW, INT and CAL products present similar precipitation estimates at the annual scale due to the quite low uncertainty associated with the precipitation in the integration. Even if the spread among the RCMs was lower than for the evaporation, some RCMs such as CRCM did compute even larger precipitation than what has been retrieved from satellites (606±80 mm yr−1). had already noted that gauges-calibrated satellite datasets such as GPCP and TMPA tend to have higher precipitation values than what was simulated in the RCMs. Precipitation over the Sea is a sensitive variable and its validation is difficult due to the lack of buoys. The ERA reanalyses value in was low compared with the NCEP estimate.

Evaporation minus precipitation focused on EP to assess the physical consistency in the RCMs. They assumed that a model with high evaporation tends to have higher precipitation. The averaged EP budget among the RCMs was 812±180 mm yr−1 and the range was between 952±80 (HIRHAM model) and 602±107 mm yr−1 (CRCM model). The inter-model spread was high for the EP budget, stressing the difficulties in providing realistic water budget evaluation. found a similar EP budget, but the associated variability was high too due to the uncertainty in evaporation. Our EP estimates are, respectively, 726±57 and 719±60 mm yr−1 before and after the closure constraint. These values are lower but still in the RCM ensemble range. They are closer to what found with NCEP reanalyses. consider the net surface flux to be 900±200 mm yr−1, which is in good agreement with the CRCM model estimate. found a similar EP budget but with far higher evaporation estimates, which seemed quite unrealistic. Furthermore, their closure constraint tends to increase the evaporation value and then the EP budget.

Discharge – only the RCMs providing the runoff have been used to compute the annual value of R (124±46 mm yr−1) in . found comparable values for the discharge, considering only the monitored rivers. did not include explicit discharge into the Mediterranean Sea since the closure was done at the global scale (Eurasian continent) and no value was provided for the Mediterranean freshwater input. Our discharge estimate is increased from 144±21 in SW to 155±15 mm yr−1 in INT after the optimization. This increase is mainly driven by the re-evaluation of the Nile discharge that present larger discharge (76 km3 yr−1) after closure. All these discharge estimates are lower than the value prescribed in (200±10 mm yr−1).

Black Sea discharge – the RCM ensemble-mean value for the freshwater input through the Bosporus Strait was 87±60 mm yr−1, stressing the high discrepancies among the RCMs. closed independently the Mediterranean and Black seas, with no exchange between the two oceanic basins (i.e., the netflow equals zero). In the current approach, the Black Sea discharge is computed as the deficit in the water budget for the Mediterranean Sea, by considering the netflow at Gibraltar (Gib) corrected from : $\mathrm{Bos}=E-P-R-\mathrm{Gib}$. The SW product presents an unrealistic value for the Black Sea discharge (2.0±615 mm yr−1); this is mainly due to the high uncertainty associated with the netflow at Gibraltar. On the other hand, the closure constraint improves the Bosporus netflow estimate, which equals 129±60 mm yr−1 with INT after optimization. The value is close to the deficit of the Black Sea water budget (computed after optimization), 132±60 mm yr−1 (not shown in Table 6), stressing the consistency between the two seas' water budget. The value is still higher than the estimate of 75 mm yr−1 in .

Gibraltar netflow considered no flow at Gibraltar when closing the Mediterranean WC and then provided no estimate for this variable. Both and evaluated the netflow by closing the WC over the Mediterranean region but they used different assumptions and equations. The estimate in is based on the oceanic closure while it is based on both the oceanic and atmospheric closure in . The RCM ensemble mean was 540±150 mm yr−1 in , while found lower value with the reanalyses (493 and 370 mm yr−1 with NCEP and ERA). give two values for the netflow at Gibraltar: one from direct observations but suffering from large uncertainties (850±400 mm yr−1), and the other as the deficit of the water budget (600±200 mm yr−1). The value in INT and CAL estimate are impacted by the closure constraint. The netflow estimate after optimization (428±124 mm yr−1) is lower than what can be found in but in the range of the RCMs water budget deficit.

textitMoisture divergence – no moisture divergence was provided by the RCMs in . found moisture divergence to be 659 mm yr−1 in NCEP and 488 mm yr−1 in ERA. estimated the divergence to be 848±105 mm yr−1 after optimization. The difference between and estimates and what is found in the current study is mainly driven by the discrepancy between the three reanalyses: Modern-Era Retrospective Analysis for ReSearch and Applications (MERRA) used by , NCEP and ERA-40 by , and ERA-I used in the current analysis. Recent works focusing on atmospheric reanalyses comparisons have demonstrated the ERA-I quality. have stressed the ERA-I performances in the representation of long term wind variability, critical for the representation of moisture divergence. have pointed out that satellite derived EP (SeaFlux-GPCP) correlates well with ERA-Interim atmospheric moisture divergence. have assessed the performance of ERA-I reanalysis for atmospheric moisture budgets consideration.

6 Conclusions

The main goal of this work was to build a multi-component dataset describing the Mediterranean hydrology by constraining the WC closure. Various methodologies have been presented and particular attention has been placed on the INTegration method. This approach fulfills the previously stated objectives: being a pixel-wise dataset but in which the WC closure is controlled. INT is an integrated dataset that shows several benefits compared to previous studies. The INT product allows reduction of the rms of the WC budget residual down to 3.55 mm month−1 over land and 5.27 mm month−1 in the atmosphere. These reductions represent an improvement of, respectively, 78 % and 80 % compared to the best unconstrained satellite combination dataset. The temporal coverage of INT is limited by the common coverage period 2004–2009 of all the satellite estimates used in this study (see Table A2).

The INT dataset has been evaluated at various scales. Even if the evaluation is a difficult task and the presented work is not exhaustive, our results show that the consideration of the WC closure allows one to reduce differences with the available in situ measurements. At the sub-basin scale, the overall precipitation is closer to the in situ gridded EOBS dataset after being constrained. The seawater-level estimate is also improved compared to the altimetry estimate. At the pixel scale, the INT estimate shows a better agreement with in situ tower measurements from the FLUXnet2015 database.

The WC has been analyzed in terms of long-term means over the 2004–2009 period and compared with the previous literature. The INT methodology has improved estimates of the Mediterranean water components. The INT product provides more realistic values for both the Bosporus and Gibraltar netflows by constraining them with the satellite observations. Note that the Bosporus estimate is mainly driven by the Gibraltar estimate and can then be improved if the Gibraltar netflow evaluation would become more accurate.

This study conducted on the Mediterranean Sea is innovative from previous work. The Mediterranean WC has already been well investigated by and relying on models and reanalyses. At the global scale, close independently the Mediterranean and Black seas using satellite observations, while close the Mediterranean Sea WC in estimating the Gibraltar netflow as the WC budget deficit. This study aims to provide a full description of the WC, based on fewer hypotheses. It is the first effort to close the WC at the surface and in the atmosphere over the whole Mediterranean region, using satellite observations and in situ measurement for the Gibraltar netflow.

There are still large uncertainties in the WC components, but the INT methodology appears to be a valuable approach, in particular to include coherency among these components. The current work has also introduced the CAL product which is a calibration of the satellite products that can be used to extrapolate in time the closure constraint. The CAL product is less efficient in closing the WC, but presents the advantage of having longer temporal coverage. Several improvements will be considered in the future: (1) more accurate in situ observations (e.g., Bosporus netflow estimate or coastal discharges) should lead to improved estimates. (2) New WC inputs could be considered (e.g., groundwater exchange or horizontal exchange at the oceanic sub-basin scale) to better characterize the flux and stock terms in the WC. (3) The use of other sources of EO estimates should be considered. For example, the evapotranspiration estimate based on the closure of the energy cycle could be tested. This dataset could be an opportunity to (4) close simultaneously the water and the energy cycles and should lead to a better estimate of the evapotranspiration over land. Multiple-component dataset INT shows promising aspects for forcing, calibrating or constraining regional models with a water conservation constraint. Some developments and evaluations are still required before the production of a climate data record can be started. The two databases (INT and CAL) can however be obtained under request to the corresponding author or via the HyMeX data server (http://mistrals.sedoo.fr/HyMeX/, last access: 24 January 2019).

Data availability
Data availability.

Data sets are available upon request by contacting the corresponding author.

Appendix A

Figure A1 compares the six water component estimates and the pdf of the two WC budget errors. The estimates are for the six terrestrial sub-basins, the oceanic part and the total land (in column) through the various methodologies presented in the study: SW, SW + PF, INT and CAL. The figures show by how much the water budget residuals are reduced and how the water components are impacted.

Table A1 gathers the notation used in the study.

Table A2 lists the datasets used in the study.

Figure A1Comparison of the six water component estimates and the pdf of the two WC budget errors (in row). The estimates are for the six terrestrial sub-basins, the oceanic part and the total land (in column) through the various methodologies presented in the study: SW, SW + PF, INT and CAL.

Table A1Notation used in this study.

Table A2Overview of the various datasets used in this study. Their common coverage period, on which the WC budget is estimated, is 2004–2009.

Author contributions
Author contributions.

VP and FA developed the methodology, analysed the results and wrote the paper. VP obtained the datasets and performed the experiments. SM and DF participated in discussions about the methodology. GJ, WD, JP and LB participated in interesting general discussions.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement
Special issue statement.

Acknowledgements
Acknowledgements.

We would like to thank the ES (European Space Agency) and its Support To Science Element program for funding the “Water Cycle Observation Multi-mission Strategy For Mediterranean region” project (ESRIN contract no. 4000114770/15/I-SBO; wacmosmed.estellus.fr). We are grateful for the E-OBS dataset from EU-FP6 project ENSEMBLES (http://ensembles-eu.metoffice.com, last access: 24 January 2019), and the data providers of the ECA & D project (http://www.ecad.eu, last access: 24 January 2019). We would like to thank the WACMOS-Med partners for the interesting related discussions and Philippe Drobinski and Véronique Ducrocq for their support to WACMOS-Med and his role in the HYMEX project.

Edited by: Eric Martin
Reviewed by: Bob Su and two anonymous referees

References

Adler, R. F., Huffman, G. J., Chang, A., Ferraro, R., Xie, P.-P., Janowiak, J., Rudolf, B., Schneider, U., Curtis, S., Bolvin, D., Gruber, A., Susskind, J., Arkin, P., and Nelkin, E.: The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979–Present), J. Hydrometeorol., 4, 1147–1167, https://doi.org/10.1175/1525-7541(2003)004<1147:TVGPCP>2.0.CO;2, 2003. a, b

Aires, F.: Combining Datasets of Satellite-Retrieved Products. Part I: Methodology and Water Budget Closure, J. Hydrometeorol., 15, 1677–1691, https://doi.org/10.1175/JHM-D-13-0148.1, 2014. a, b, c, d, e, f, g, h, i

Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: Crop evapotranspiration: Guidelines for computing crop water requirements, FAO irrigation and drainage paper 56, FAO, Rome, 300, 6541, 1998. a

Ashouri, H., Hsu, K.-L., Sorooshian, S., Braithwaite, D. K., Knapp, K. R., Cecil, L. D., Nelson, B. R., Prat, O. P., Ashouri, H., Hsu, K.-L., Sorooshian, S., Braithwaite, D. K., Knapp, K. R., Cecil, L. D., Nelson, B. R., and Prat, O. P.: PERSIANN-CDR: Daily Precipitation Climate Data Record from Multisatellite Observations for Hydrological and Climate Studies, B. Am. Meteorol. Soc., 96, 69–83, https://doi.org/10.1175/BAMS-D-13-00068.1, 2015. a, b

Bahurel, P., Ocean, M., Bell, M.-O., and Le Traon, P.: Ocean Monitoring and Forecasting in the European MyOcean GMES Marine Initiative, available at: http://marine.copernicus.eu/wp-content/uploads/2016/06/r652_9_pierre_bahurel_talk.pdf (last access: 24 January 2019), 2012. a

Bettadpur, S.: GRACE 327-742 (CSR-GR-12-xx) Gravity Recovery and Climate Experiment UTCSR Level-2 Processing Standards Document) (For Level-2 Product Release 0005), Technical report, Center for Space Research, The University of Texas, Austin, 2012. a, b

Biancale, R., Lemoine, J.-M., Bruinsma, S., and Bourgogne, S.: 10-day time series of the geoid from GRACE and LAGEOS data Evaluation of solutions Processing standards Mean field and 10-day solutions Noise estimation over global oceans, CNES, Toulouse, France, available at: ftp://ftp.csr.utexas.edu/pub/grace/Proceedings/GSTM_2009_Extras/Posters/Biancale_Poster.pdf (last access: 24 January 2019), 2005. a, b

Brown, P. J. and Kummerow, C. D.: An Assessment of Atmospheric Water Budget Components over Tropical Oceans, J. Climate, 27, 2054–2071, https://doi.org/10.1175/JCLI-D-13-00385.1, 2013. a

Chen, X., Su, Z., Ma, Y., Yang, K., Wen, J., and Zhang, Y.: An improvement of roughness height parameterization of the Surface Energy Balance System (SEBS) over the Tibetan plateau, J. Appl. Meteorol. Clim., 52, 607–622, https://doi.org/10.1175/JAMC-D-12-056.1, 2013. a

Curry, J. A., Bentamy, A., Bourassa, M. A., Bourras, D., Bradley, E. F., Brunke, M., Castro, S., Chou, S. H., Clayson, C. A., Emery, W. J., Eymard, L., Fairall, C. W., Kubota, M., Lin, B., Perrie, W., Reeder, R. A., Renfrew, I. A., Rossow, W. B., Schulz, J., Smith, S. R., Webster, P. J., Wick, G. A., Zeng, X., Curry, J. A., Bentamy, A., Bourassa, M. A., Bourras, D., Bradley, E. F., Brunke, M., Castro, S., Chou, S. H., Clayson, C. A., Emery, W. J., Eymard, L., Fairall, C. W., Kubota, M., Lin, B., Perrie, W., Reeder, R. A., Renfrew, I. A., Rossow, W. B., Schulz, J., Smith, S. R., Webster, P. J., Wick, G. A., and Zeng, X.: SEAFLUX, B. Am. Meteorol. Soc., 85, 409–424, https://doi.org/10.1175/BAMS-85-3-409, 2004. a, b

Dahle, C., Flechtner, F., Gruber, C., König, D., König, R., Michalak, G., and Neumayer, K.-H.: GFZ GRACE Level-2 Processing Standards Document for Level-2 Product Release 0005, Scientific Technical Report – Data; 12/02, Deutsches GeoForschungsZentrum GFZ, Potsdam, p. 20, https://doi.org/10.2312/GFZ.b103-12020, 2013. a, b

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteorol. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011. a, b, c, d, e, f, g

De Geeter, J., Van Brussel, H., De Schutter, J., and Decreton, M.: A smoothly constrained Kalman filter, IEEE Trans. Pattern Anal. Mach. Intell., 19, 1171–1177, https://doi.org/10.1109/34.625129, 1997. a

Ducharne, A., Golaz, C., Leblois, E., Laval, K., Polcher, J., Ledoux, E., and de Marsily, G.: Development of a high resolution runoff routing model, calibration and application to assess runoff from the LMD GCM, J. Hydrol., 280, 207–228, https://doi.org/10.1016/S0022-1694(03)00230-0, 2003. a

Eicker, A., Forootan, E., Springer, A., and Longuevergne, L.: Does GRACE see the terrestrial water cycle `intensifying'?, J. Geophys. Res., 121, 733–745, 2016. a

Falge, E., Aubinet, M., Bakwin, P., Baldocchi, D., Berbigier, P., Bernhofer, C., Black, T., Ceulemans, R., Davis, K., Dolman, A., Goldstein, A., Goulden, M., Granier, A., Hollinger, D., Jarvis, P., Jensen, N., Pilegaard, K., Katul, G., Kyaw Tha Paw, P., Law, B., Lindroth, A., Loustau, D., Mahli, Y., Monson, R., Moncrieff, P., Moors, E., Munger, J., Meyers, T., Oechel, W., Schulze, E.-D., Thorgeirsson, H., Tenhunen, J., Valentini, R., Verma, S., Vesala, T., and Wofsy, S.: FLUXNET Research Network Site Characteristics, Investigators, and Bibliography, 2016, https://doi.org/10.3334/ORNLDAAC/1530, 2017. a, b, c

Fenoglio-Marc, L., Kusche, J., and Becker, M.: Mass variation in the Mediterranean Sea from GRACE and its validation by altimetry, steric and hydrologic fields, Geophys. Res. Lett., 33, 1–5, https://doi.org/10.1029/2006GL026851, 2006. a

Good, S. A., Martin, M. J., and Rayner, N. A.: EN4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates, J. Geophys. Res.-Oceans, 118, 6704–6716, https://doi.org/10.1002/2013JC009067, 2013. a

Hamon, M., Beuvier, J., Somot, S., Lellouche, J.-M., Greiner, E., Jordà, G., Bouin, M.-N., Arsouze, T., Béranger, K., Sevault, F., Dubois, C., Drevillon, M., and Drillet, Y.: Design and validation of MEDRYS, a Mediterranean Sea reanalysis over the period 1992–2013, Ocean Sci., 12, 577–599, https://doi.org/10.5194/os-12-577-2016, 2016. a

Harzallah, A., Jordà, G., Dubois, C., Sannino, G., Carillo, A., Li, L., Arsouze, T., Cavicchia, L., Beuvier, J., and Akhtar, N.: Long term evolution of heat budget in the Mediterranean Sea from Med-CORDEX forced and coupled simulations, Clim. Dynam., 51, 1145, https://doi.org/10.1007/s00382-016-3363-5, 2018. a

Haylock, M. R., Hofstra, N., Klein Tank, A. M. G., Klok, E. J., Jones, P. D., and New, M.: A European daily high-resolution gridded data set of surface temperature and precipitation for 1950–2006, J. Geophys. Res., 113, D20119, https://doi.org/10.1029/2008JD010201, 2008. a, b

Hesselbjerg Christensen, J. and Meteorologisk Institut: The HIRHAM4 regional atmospheric climate model, DMI, Copenhagen, 1996. a

Huffman, G. J., Bolvin, D. T., Nelkin, E. J., Wolff, D. B., Adler, R. F., Gu, G., Hong, Y., Bowman, K. P., and Stocker, E. F.: The TRMM Multisatellite Precipitation Analysis (TMPA): Quasi-Global, Multiyear, Combined-Sensor Precipitation Estimates at Fine Scales, J. Hydrometeorol., 8, 38–55, https://doi.org/10.1175/JHM560.1, 2007. a, b

IPCC: IPCC, 2014: Climate Change 2014: Synthesis Report, in: Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Tech. rep., edited by: Core Writing Team, Pachauri, R. K., and Meyer, L. A., IPCC, Geneva, Switzerland, 2014. a

Ishii, M., Kimoto, M., and Kachi, M.: Historical Ocean Subsurface Temperature Analysis with Error Estimates, Mon. Weather Rev., 131, 51–73, https://doi.org/10.1175/1520-0493(2003)131<0051:HOSTAW>2.0.CO;2, 2003. a

Jordà, G. and Gomis, D.: On the interpretation of the steric and mass components of sea level variability: The case of the Mediterranean basin, J. Geophys. Res.-Oceans, 118, 953–963, https://doi.org/10.1002/jgrc.20060, 2006. a

Jordà, G., Sánchez-Román, A., and Gomis, D.: Reconstruction of transports through the Strait of Gibraltar from limited observations, Clim. Dynam., 48, 851–865, https://doi.org/10.1007/s00382-016-3113-8, 2016. a, b, c, d, e

Jordà, G., Von Schuckmann, K., Josey, S. A., Caniaux, G., García-Lafuente, J., Sammartino, S., Özsoy, E., Polcher, J., Notarstefano, G., Poulain, P. M., Adloff, F., Salat, J., Naranjo, C., Schroeder, K., Chiggiato, J., Sannino, G., and Macías, D.: The Mediterranean Sea heat and mass budgets: Estimates, uncertainties and perspectives, Prog. Oceanogr., 156, 174–208, https://doi.org/10.1016/j.pocean.2017.07.001, 2017. a, b, c, d, e, f, g, h, i, j, k

Joyce, R. J., Janowiak, J. E., Arkin, P. A., and Xie, P.: CMORPH: A Method that Produces Global Precipitation Estimates from Passive Microwave and Infrared Data at High Spatial and Temporal Resolution, J. Hydrometeorol., 5, 487–503, 2004. a, b

Kottek, M., Grieser, J., Beck, C., Rudolf, B., and Rubel, F.: World Map of the Köppen–Geiger climate classification updated, Meteorol. Z., 15, 259–263, https://doi.org/10.1127/0941-2948/2006/0130, 2006. a

L'Ecuyer, T. S., Beaudoing, H. K., Rodell, M., Olson, W., Lin, B., Kato, S., Clayson, C. A., Wood, E., Sheffield, J., Adler, R., Huffman, G., Bosilovich, M., Gu, G., Robertson, F., Houser, P. R., Chambers, D., Famiglietti, J. S., Fetzer, E., Liu, W. T., Gao, X., Schlosser, C. A., Clark, E., Lettenmaier, D. P., and Hilburn, K.: The Observed State of the Energy Budget in the Early Twenty-First Century, J. Climate, 28, 8319–8346, https://doi.org/10.1175/JCLI-D-14-00556.1, 2015. a, b

Ludwig, W., Dumont, E., Meybeck, M., and Heussner, S.: River discharges of water and nutrients to the Mediterranean and Black Sea: Major drivers for ecosystem changes during past and future decades?, Prog. Oceanogr., 80, 199–217, https://doi.org/10.1016/j.pocean.2009.02.001, 2009. a, b

Margat, J.: Programme des Nations Unies pour l'environnement plan d'action pour la Méditerranée: situation et perspectives, MAP Technical Report Series No. 158 PNUE/PAM, 2004 Programme des Nations Unies pour l'environnement/Plan d'action pour la Méditerranée (PNUE/PAM), Athènes, Grèce, 2004. a

Mariotti, A., Struglia, M. V., Zeng, N., and Lau, K.-M.: The Hydrological Cycle in the Mediterranean Region and Implications for the Water Budget of the Mediterranean Sea, J. Climate, 15, 1674–1690, https://doi.org/10.1175/1520-0442(2002)015<1674:THCITM>2.0.CO;2, 2002. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

Martens, B., Miralles, D. G., Lievens, H., van der Schalie, R., de Jeu, R. A. M., Fernández-Prieto, D., Beck, H. E., Dorigo, W. A., and Verhoest, N. E. C.: GLEAM v3: satellite-based land evaporation and root-zone soil moisture, Geosci. Model Dev., 10, 1903–1925, https://doi.org/10.5194/gmd-10-1903-2017, 2017. a, b

McCabe, M., Wood, E., Wójcik, R., Pan, M., Sheffield, J., Gao, H., and Su, H.: Hydrological consistency using multi-sensor remote sensing data for water and energy cycle studies, Remote Sens. Environ., 112, 430–444, https://doi.org/10.1016/j.rse.2007.03.027, 2008. a

Miralles, D. G., De Jeu, R. A. M., Gash, J. H., Holmes, T. R. H., and Dolman, A. J.: Magnitude and variability of land evaporation and its components at the global scale, Hydrol. Earth Syst. Sci., 15, 967–981, https://doi.org/10.5194/hess-15-967-2011, 2011. a

Mu, Q., Zhao, M., and Running, S. W.: Improvements to a MODIS global terrestrial evapotranspiration algorithm, Remote Sens. Environ., 115, 1781–1800, https://doi.org/10.1016/j.rse.2011.02.019, 2011. a, b

Munier, S. and Aires, F.: A new global method of satellite dataset merging and quality characterization constrained by the terrestrial water cycle budget, Remote Sens. Environ., 205, 119–130, https://doi.org/10.1016/j.rse.2017.11.008, 2017. a, b, c, d, e, f

Munier, S., Aires, F., Schlaffer, S., Prigent, C., Papa, F., Maisongrande, P., and Pan, M.: Combining data sets of satellite-retrieved products for basin-scale water balance study: 2. Evaluation on the Mississippi Basin and closure correction model, J. Geophys. Res., 119, 12100–12116, https://doi.org/10.1002/2014JD021953, 2014. a, b, c, d, e, f, g, h, i

Pan, M. and Wood, E. F.: Data Assimilation for Estimating the Terrestrial Water Budget Using a Constrained Ensemble Kalman Filter, J. Hydrometeorol., 7, 534–547, https://doi.org/10.1175/JHM495.1, 2006. a, b, c

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

Pellet, V., Aires, F., Mariotti, A., and Fernandez, D.: Analyzing the Mediterranean water cycle via satellite data integration, Pure Appl. Geophys., 175, 3909–3937, 2017. a, b, c, d, e, f, g, h, i, j

Pike, J.: The estimation of annual run-off from meteorological data in a tropical climate, J. Hydrol., 2, 116–123, https://doi.org/10.1016/0022-1694(64)90022-8, 1964. a

Plummer, D. A., Caya, D., Frigon, A., Côté, H., Giguère, M., Paquin, D., Biner, S., Harvey, R., de Elia, R., Plummer, D. A., Caya, D., Frigon, A., Côté, H., Giguère, M., Paquin, D., Biner, S., Harvey, R., and de Elia, R.: Climate and Climate Change over North America as Simulated by the Canadian RCM, J. Climate, 19, 3112–3132, https://doi.org/10.1175/JCLI3769.1, 2006. a

Polcher, J., Mcavaney, B., Viterbo, P., Gaertner, M.-A., Hahmann, A., Mahfouf, J.-F., Noilhan, J., Phillips, T., Pitman, A., Schlosser, C. A., Schulz, J.-P., Timbal, B., Verseghy, D., and Xue, Y.: A proposal for a general interface between land surface schemes and general circulation models, Global Planet. Change, 19, 261–276, 1998. a, b

Rodell, M., Beaudoing, H., L'Ecuyer, T., Olson, W., Famiglietti, J., Houser, P., Adler, R., Bosilovich, M., Clayson, C., Chambers, D., Clark, E., Fetzer, E., Gao, X., Gu, G., Hilburn, K., Huffman, G., Lettenmaier, D., Liu, W., Robertson, F., Schlosser, C., Sheffield, J., and Wood, E.: The Observed State of the Water Cycle in the Early 21st Century, J. Climate, 28, 8289–8318, https://doi.org/10.1175/JCLI-D-14-00555.1, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y

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

Sanchez-Gomez, E., Somot, S., Josey, S. A., Dubois, C., Elguindi, N., and Déqué, M.: Evaluation of Mediterranean Sea water and heat budgets simulated by an ensemble of high resolution regional climate models, Clim. Dynam., 37, 2067–2086, https://doi.org/10.1007/s00382-011-1012-6, 2011. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Schneider, N., Schröder, M., Lindstrot, R., Preusker, R., Stengel, M., and ESA DUE GlobVapour Consortium: ESA DUE GlobVapour water vapor products: Validation, AIP Conference Proceedings, 1531, 484–487, https://doi.org/10.1063/1.4804812, 2013. a, b

Seager, R. and Henderson, N.: Diagnostic computation of moisture budgets in the ERA-interim reanalysis with reference to analysis of CMIP-archived atmospheric model data, J. Climate, 26, 7876–7901, https://doi.org/10.1175/JCLI-D-13-00018.1, 2013. a, b

Sheffield, J., Ferguson, C. R., Troy, T. J., Wood, E. F., and McCabe, M. F.: Closing the terrestrial water budget from satellite remote sensing, Geophys. Res. Lett., 36, 1–5, https://doi.org/10.1029/2009GL037338, 2009. a

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

Stopa, J. E. and Cheung, K. F.: Intercomparison of wind and wave data from the ECMWF Reanalysis Interim and the NCEP Climate Forecast System Reanalysis, Ocean Model., 75, 65–83, https://doi.org/10.1016/j.ocemod.2013.12.006, 2014. a, b

Su, Z.: The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes, Hydrol. Earth Syst. Sci., 6, 85–100, https://doi.org/10.5194/hess-6-85-2002, 2002. a

Su, Z., Timmermans, W., Zeng, Y., Schulz, J., John, V. O., Roebeling, R. A., Poli, P., Tan, D., Kaspar, F., Kaiser-Weiss, A. K., Swinnen, E., Toté, C., Gregow, H., Manninen, T., Riihelä, A., Calvet, J. C., Ma, Y., and Wen, J.: An overview of european efforts in generating climate data records, B. Am. Meteorol. Soc., 99, 349–359, https://doi.org/10.1175/BAMS-D-16-0074.1, 2018. a

Sun, B., Yu, L., and Weller, R. A.: Comparisons of Surface Meteorology and Turbulent Heat Fluxes over the Atlantic: NWP Model Analyses versus Moored Buoy Observations, J. Climate, 16, 679–695, https://doi.org/10.1175/1520-0442(2003)016<0679:COSMAT>2.0.CO;2, 2003. a, b

Trenberth, K. E., Fasullo, J. T., and Mackaro, J.: Atmospheric Moisture Transports from Ocean to Land and Global Energy Flows in Reanalyses, J. Climate, 24, 4907–4924, https://doi.org/10.1175/2011JCLI4171.1, 2011. a

Troy, T. J. and Wood, E. F.: Comparison and evaluation of gridded radiation products across northern Eurasia, Environ. Res. Lett., 4, 045008, https://doi.org/10.1088/1748-9326/4/4/045008, 2009. a

Troy, T. J., Sheffield, J., Wood, E. F., Troy, T. J., Sheffield, J., and Wood, E. F.: Estimation of the Terrestrial Water Budget over Northern Eurasia through the Use of Multiple Data Sources, J. Climate, 24, 3272–3293, https://doi.org/10.1175/2011JCLI3936.1, 2011. a

Watkins, M. M. and Yuan, D.-N.: GRACE Gravity Recovery and Climate Experiment JPL Level-2 Processing Standards Document For Level-2 Product Release 05.1, Technical report, Jet Propulsion Laboratory, California Institute of Technology, 2014. a, b

Watkins, M. M., Wiese, D. N., Yuan, D.-N., Boening, C., and Landerer, F. W.: Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons, J. Geophys. Res.-Solid, 120, 2648–2671, https://doi.org/10.1002/2014JB011547, 2015. a, b

Wu, H., Kimball, J. S., Mantua, N., and Stanford, J.: Automated upscaling of river networks for macroscale hydrological modeling, Water Resour. Res., 47, W03517, https://doi.org/10.1029/2009WR008871, 2011. a, b

Yilmaz, M. T., DelSole, T., and Houser, P. R.: Improving Land Data Assimilation Performance with a Water Budget Constraint, J. Hydrometeorol., 12, 1040–1055, https://doi.org/10.1175/2011JHM1346.1, 2011. a, b

Zhang, K., Kimball, J. S., Nemani, R. R., and Running, S. W.: A continuous satellite-derived global record of land surface evapotranspiration from 1983 to 2006, Water Resour. Res., 46, W09522, https://doi.org/10.1029/2009WR008800, 2010.  a, b

Zhang, Y., Pan, M., and Wood, E. F.: On Creating Global Gridded Terrestrial Water Budget Estimates from Satellite Remote Sensing, Surv. Geophys., 37, 1–20, https://doi.org/10.1007/s10712-015-9354-y, 2016. a, b

Zhang, Y., Pan, M., Sheffield, J., Siemann, A. L., Fisher, C. K., Liang, M., Beck, H. E., Wanders, N., MacCracken, R. F., Houser, P. R., Zhou, T., Lettenmaier, D. P., Pinker, R. T., Bytheway, J., Kummerow, C. D., and Wood, E. F.: A Climate Data Record (CDR) for the global terrestrial water budget: 1984–2010, Hydrol. Earth Syst. Sci., 22, 241–263, https://doi.org/10.5194/hess-22-241-2018, 2018. a, b

FLUXNET2015 datasets; https://fluxnet.fluxdata.org (last access: 24 January 2019).

CMEMS http://www.marine.copernicus.eu (last access: 24 January 2019).

Special issue