High-resolution satellite products improve hydrological modeling in northern Italy

Satellite-based Earth observations (EO) are an accurate and reliable data source for atmospheric and environmental science. Their increasing spatial and temporal resolutions, as well as the seamless availability over ungauged regions, make them appealing for hydrological modeling. This work shows recent advances in the use of high-resolution satellite-based EO data in hydrological modeling. In a set of six experiments, the distributed hydrological model Continuum is set up for the Po River basin (Italy) and forced, in turn, by satellite precipitation and evaporation, while satellite-derived soil moisture (SM) and snow depths are ingested into the model structure through a data-assimilation scheme. Further, satellite-based estimates of precipitation, evaporation, and river discharge are used for hydrological model calibration, and results are compared with those based on ground observations. Despite the high density of conventional ground measurements and the strong human influence in the focus region, all satellite products show strong potential for operational hydrological applications, with skillful estimates of river discharge throughout the model domain. Satellite-based evaporation and snow depths marginally improve (by 2 % and 4 %) the mean Kling–Gupta efficiency (KGE) at 27 river gauges, compared to a baseline simulation (KGE_mean= 0.51) forced by high-quality conventional data. Precipitation has the largest impact on the model output, though the satellite data on average shows poorer skills compared to conventional data. Interestingly, a model calibration heavily relying on satellite data, as opposed to conventional data, provides a skillful reconstruction of river discharges, paving the way to fully satellite-driven hydrological applications.

1 Introduction

Remote sensing of the Earth from space is a ripe yet ever growing sector, with countless applications and users worldwide. Hydrological sciences have already benefited enormously from Earth observation (EO) data (see e.g., McCabe et al., 2017; Chen and Wang, 2018; Alfieri et al., 2018), thanks to global and independent datasets for the different components of the water and energy cycles as well as anthropogenic processes such as irrigation (Massari et al., 2021). Hydrological models play a crucial role in monitoring and forecasting, thanks to their ability to reproduce the physical processes that govern the water cycle. Their successful implementation is strongly conditioned by the availability of consistent, accurate, and seamless hydrometeorological datasets for the considered focus region, space/time resolution, and period of interest. Conventional data, including ground observations and weather radars, are traditionally popular sources of dynamic data to force these models. Yet, they are not viable options for the still vast ungauged regions of the world. Satellite and independent products offer a range of alternatives to fill such gaps, thanks to their massive contribution to the atmospheric reanalyses (see e.g., Hersbach et al., 2020) as well as with independent products. Hydrological models can benefit from dynamic data (either ground or satellite-based) in various forms: (1) as forcing datasets, (2) as assimilation datasets, (3) as benchmark data for model calibration and improved parameterization, and (4) to investigate process understanding.

Forcing data are mandatory input for hydrological models. Key variables are precipitation, air temperature, and evaporation or, alternatively, the meteorological variables needed to estimate them. Their influence on hydrological modeling was assessed, for example, by Wu et al. (2017) and Beck et al. (2017) for precipitation datasets, Dembélé et al. (2020a) for evaporation datasets, and Dembélé et al. (2020b) for combinations of temperature and precipitation datasets. The latter found a reduced influence of the choice of temperature datasets on the output discharge, though these can significantly impact evaporation and soil moisture (SM) estimates. Data assimilation methods are designed to merge measurements of any type with estimates from geophysical models (Reichle, 2008), to compensate for errors in the forcing data, model structural deficiencies, and update their state variables at the initial or intermediate simulation steps (Spaaks and Bouten, 2013). Relevant applications of assimilating satellite products in hydrological modeling include SM (Massari et al., 2015; Wanders et al., 2014), water storage (Li et al., 2012), snow cover (Thirel et al., 2013), evaporation (Hartanto et al., 2017), land surface temperature (Silvestro et al., 2013), water levels (Paiva et al., 2013), discharge (Ishitsuka et al., 2020), water extent (Revilla-Romero et al., 2016; Hostache et al., 2018), and multi-variable combinations (Wongchuig-Correa et al., 2020). Hydrometeorological data have also been used as a benchmark to train the model parameters through machine-learning techniques (Mosaffa et al., 2022) or calibration techniques based on minimization of cost functions computed between simulated and observed variables (Pechlivanidis et al., 2011; Demirel et al., 2018). Additionally, satellite estimation of river levels shows promising applications in the field. It has been tested in the calibration of hydrological (Getirana et al., 2013; Dhote et al., 2021) and hydraulic (Domeneghetti et al., 2021) models.

As part of the Green Deal and Digital Strategy, the European Commission recently launched the Destination Earth program (https://digital-strategy.ec.europa.eu/en/policies/destination-earth, last access: 28 July 2022). This is a joint effort involving key European institutions to develop a very high-precision digital model, or “Digital Twin”, of the Earth for monitoring and predicting environmental change and human impacts, to ultimately support sustainable development. The present work strives toward that direction by contributing to the development of a Digital Twin Earth that is focused on the water cycle and hydrological processes. It highlights the potential of high-resolution satellite products in describing the water cycle and monitoring hydrological extremes and water resources. Through various dedicated experiments, we test the influence of five new high-resolution satellite-derived datasets on the performance of CIMA's distributed hydrological model Continuum (Silvestro et al., 2013), set up for the entire Po River basin in northern Italy. These include (1) GPM-SM2RAIN (Massari et al., 2020) precipitation and (2) the Global Land Evaporation Amsterdam Model (GLEAM; Miralles et al., 2011) evaporation as dynamic forcing; data assimilation of (3) C-SNOW (Lievens et al., 2019) snow depth and (4) RT1 (Quast et al., 2019) SM; and model calibration using (5) satellite-based river discharge (Tarpanelli et al., 2020) as a benchmark. By comparing results with observed river discharge over 2017–2019, and with a simulation forced by conventional data, we investigate the relative impact of these high-resolution satellite products. Further, we take the first steps toward hydrological modeling, fully relying on satellite data, by calibrating and subsequently running the model using SM2RAIN satellite precipitation and GLEAM evaporation as forcing, and satellite-based estimates of river discharge as benchmark data for the calibration.

2 Case study and data

2.1 Case study – the Po River basin

The Po River basin has a catchment area of about 74 000 km² that is shared between Italy (95 %) and Switzerland (5 %). It is fed by tributaries from the Alps in the north and west, and by the Apennines in the south. The basin elevation ranges between 4800 m and the sea level, and hence it features a variety of climatic and hydrological regimes, from a glacial and snow-rain type in the mountain area to a pluvial yet drier regime in the lowland section. The region is considered to be highly vulnerable to flooding, both economically and with respect to loss of life (Domeneghetti et al., 2015). The basin plays a significant role in the Italian economy, hosting approximately 25 % of the Italian population, producing 40 % of the national GDP, and consuming 48 % of national produced energy. The Po River flows through the Po Plain, one of the largest contiguous agricultural areas of Europe. This causes more than 30 % of water to be extracted from surface water and used for agricultural purposes. Although water is sufficient for all uses under average climate conditions, recent periods of prolonged drought led to substantial economic losses and threats to water security (Mysiak et al., 2013), thus a comprehensive evaluation of the impacts of human activities on water resources in the area is a far-reaching matter. Given its large socio-economic influence, the Po River basin has already been investigated through a number of modeling approaches forced by in situ data and by Earth system models, especially to predict the impact of inundation and of climate change (e.g., Ravazzani et al., 2015; Vezzoli et al., 2015; Nogherotto et al., 2019) while applications including satellite products remain scarce.

2.2 Static data

In the choice of spatial information, large-scale datasets were deliberately used over more detailed local data, in line with the concept of the Digital Twin Earth and in view of the plan to extend the simulation area for a continental or global application. We used the Digital Elevation Model (DEM) from the global USGS Hydrologic Derivatives for Modeling and Analysis (HDMA; Verdin, 2017) at 3 arcsec spatial resolution (about 90 m at the Equator), which comes with pre-computed and corrected hydrological derivatives including channel network and macro basins. The DEM was upscaled at the chosen model resolution of 1 km through cubic resampling, to define the computational grid and compute the necessary hydrological derivatives (flow accumulation, drainage directions, and channel network). The river network is defined by cells with an upstream area larger than 240 km², following previous applications of Continuum in northern Italy. To improve its spatial representation, DEM was carved with a high-resolution stream network of the main rivers taken from the Italian Institute for Environmental Protection and Research, while dikes were manually placed at specific locations, especially in flat areas.

The curve number map used to model direct runoff and infiltration from excess rainfall was derived from the ESA-CCI 2018 land cover map (ESA, 2017) at 300 m resolution, together with information on the soil characteristics. Hydrologic soil groups were extracted from the HYSOGs250m (Ross et al., 2018), while for soil texture identification, we applied the USDA method (Shirazi and Boersma, 1984) using the global maps of the fractions of sand and clay from the International Soil Reference and Information Centre (ISRIC) SoilGrids (Hengl et al., 2017), combined with the ESA CCI SoilMoisture (Dorigo et al., 2017) global map of soil porosity. Glacial areas used in the cryospheric model S3M (see Sect. 3.1) were taken from the Randolph Glacier Inventory (RGI) v6 (Raup et al., 2007). Vegetation coverage is taken from the global land cover map ECOCLIMAP (Faroux et al., 2013).

Point information for a set of 99 reservoirs and the three major lakes (Maggiore, Como, and Garda) was included in the model setup (Fig 1). Information on the dams and the corresponding reservoirs was provided by the Italian Civil Protection Department (DPC) and from the Global Reservoir and Dam (GRanD) database (Lehner et al., 2011). Data ingested for each dam include the maximum stored volume, initial volume, maximum non-damaging discharge at the outflow gates, weir length, maximum storage level, outflow coefficient, and coordinates of the release point. For lakes, required metadata are the outlet coordinates, minimum volume inducing outflow discharge, initial volume, and emptying coefficient.

Figure 1 Simulated domain (blue line) and river network (dark green) of the Po River basin. Symbols show the point features considered in the hydrological model.

2.3 Dynamic data

The hydrological model used requires input maps of precipitation, air temperature, relative humidity, wind speed and incoming solar radiation. Alternatively, both actual and potential evaporation can be provided as dynamic input, where the latter is used to estimate actual evaporation from lakes and reservoirs. In such cases wind speed maps are not needed by the model. The baseline hydrological simulation uses conventional meteorological data as input. Precipitation fields were estimated with the Modified Conditional Merging (MCM) technique (Bruno et al., 2021), which incorporates precipitation gauges and radar estimates. The MCM is an improvement of the conditional merging (CM) technique proposed by Sinclair and Pegram (2005), which estimates the structure of covariance and the length of spatial correlation at every gauge, taking it from the cumulated radar precipitation fields. For the Po River basin, MCM is based on 1377 precipitation gauges and on the mosaic of the Italian weather radars.

Hourly maps of the weather variables collected for the Po River basin ultimately include 1258 temperature stations, 608 for relative humidity, 460 for wind speed, and 278 for solar radiation. Temperature maps include an altitude correction algorithm with temperature gradients estimated at every time step by linearly interpolating available data at different elevations. They also include an outlier removal algorithm which discards station data with a deviation of more than 20 ^∘C from the corresponding temperature–elevation interpolating line.

Discharge data at 27 river gauging stations with hourly sampling frequency for the years 2016–2019 were provided by DPC and the regional hydrometeorological offices. Twenty-two stations were selected for model calibration, while five were retained for validation only (Fig. 1). Validation stations were chosen to represent different areas of the Po River basin, including a mix of small and large sub-catchments with varying influences on lakes and reservoirs.

2.4 Satellite products and validation

2.4.1 Precipitation

The precipitation dataset used in this work and referred to as SM2RAIN (Fig. 2) merges SM2RAIN-ASCAT (Brocca et al., 2019) and GPM mission IMERG-LR (Huffman et al., 2015) datasets, both available at 10 km spatial resolution. Unlike the work of Massari et al. (2020), where the fusion of the two datasets was based on an optimal interpolation technique, here we relied on a triple collocation (TC)-based merging using the signal-to-noise ratio (SNR), as in Gruber et al. (2017). In particular, to derive the merged dataset we seek the optimality in a least squares sense, so that the variance of residual random errors is minimized. This leads to a weighted average between SM2RAIN-ASCAT and IMERG-LR, i.e.,

$$\begin{array}{}\text{(1)}& P_{\mathrm{DTE}}=w_{\mathrm{1}}{P}_{\text{SM2RAIN-ASCAT}}+w_{\mathrm{2}}{P}_{\text{IMERG-LR}},\end{array}$$

where the weights w₁ and w₂ are calculated as

$$\begin{array}{}\text{(2)}& w_{\mathrm{1}}={\displaystyle \frac{{\mathrm{SNR}}_{\mathrm{1}}}{{\mathrm{SNR}}_{\mathrm{1}}+{\mathrm{SNR}}_{\mathrm{2}}}}{w}_{\mathrm{2}}={\displaystyle \frac{{\mathrm{SNR}}_{\mathrm{2}}}{{\mathrm{SNR}}_{\mathrm{1}}+{\mathrm{SNR}}_{\mathrm{2}}}},\end{array}$$

where SNR is estimated as the ratio between the variance of the true signal and that of the considered satellite product, multiplied by a parameter representing the systematic error (see Gruber et al., 2017). The subscripts 1 and 2 refer to the SM2RAIN-ASCAT and IMERG-LR datasets, respectively. Under the assumption that the two datasets are independent (as also required by TC), the random error of the merged time series is lower than those of the individual input datasets.

Figure 2 Daily Pearson correlation coefficient between SM2RAIN precipitation and the MCM (radar-gauge) precipitation dataset. The median correlation r_SM2RAIN= 0.76 largely improves that of the two individual products, i.e., r_{SM2RAIN-ASCAT}= 0.66 and r_IMERG-LR= 0.67.

The TC was applied to the triplet: SM2RAIN-ASCAT, IMERG-LR, and the MCM radar-gauge precipitation dataset. Note that, unlike the use of random error variances in Crow et al. (2015), weights calculated in Eq. (2) do not require the assumption of null systematic differences between the datasets, thanks to the self-consistency of SNR (see Gruber et al., 2017 for further details). Before the weights can be used to merge the datasets, relative systematic differences (i.e., long-term bias) have to be corrected to make the weights obtained by Eq. (2) converge to the optimal weights in a least squares sense (Crow et al., 2015). Given the nature of the precipitation signal (containing many null values), this rescaling has been done by means of a multiplicative factor to the mean with respect to MCM. The fusion of the two datasets was only done for the time steps where IMERG-LR was greater than zero. Due to the high sensitivity of the GPM mission, values with zero precipitation in IMERG-LR were set to zero. Hourly data were obtained by imposing the sub-daily temporal pattern of IMERG-LR to the merged dataset. The 10 km resolution dataset generated was thus resampled at 1 km resolution through bilinear interpolation for use in the hydrological model.

2.4.2 Evaporation

GLEAM (Miralles et al., 2011) is a state-of-the-art methodology to derive evaporation and its various components (i.e., transpiration, bare soil evaporation, interception loss, snow sublimation, and open-water evaporation). It combines global satellite observations of meteorological (precipitation, near-surface net radiation, and air temperature) and surface (soil and vegetation water content, and snow water equivalent (SWE)) variables that are informative for the evaporation process. The model is based on the equation by Priestley and Taylor (1972) to estimate potential evaporation. Those estimates are then constrained based on root-zone SM, which results from a precipitation-driven running water balance in which satellite-based SM can be assimilated. Interception loss is independently estimated through an adapted Gash analytical model (Miralles et al., 2010). Since its first version, GLEAM has been widely deployed at coarse resolution for climatic studies. In the past few years, it has been further developed to solve higher spatial and temporal resolutions. For example, Martens et al. (2018) obtained accurate results in an implementation over the Netherlands at 100 m resolution. For this work, GLEAM was applied over the entire Po River basin to produce both potential and actual evaporation estimates at 1 km resolution.

Since measurements of evaporation in the focus region are limited, the performance of the 1 km evaporation dataset was inferred on the basis of the FluxNet IT-Tor site, located in the mountainous Val d'Aosta region in the northwestern part of the domain (Fig. 3). While based on one station only, the performance (Pearson correlation coefficient r= 0.83) is in line with results obtained in the high-resolution implementation across the Netherlands, where Martens et al. (2018) found a median temporal correlation coefficient of 0.76 across 29 sites.

Figure 3 Average GLEAM daily actual evaporation for 2017–2019 in the Po River basin at 1 km resolution (a) and comparison with daily evaporation from the FluxNet site IT-Tor for 2018 (b). Pearson correlation coefficient r= 0.83. The location of the FluxNet site is marked with a red star in (a).

2.4.3 Soil moisture (SM)

High-resolution SM was retrieved from incidence angle-dependent Sentinel-1 backscatter measurements at 500 m spatial sampling (∼ 1 km spatial resolution) (Bauer-Marschallinger et al., 2019) by using a modeling approach based on the time series of first-order radiative transfer (RT1; see Quast and Wagner, 2016; Quast et al., 2019). The RT1 model uses auxiliary leaf area index (LAI) time series provided by the ECMWF ERA5-Land reanalysis dataset (Muñoz-Sabater et al., 2021) to correct effects induced by seasonal dynamics of vegetation. The retrieval is then performed via a nonlinear least-squares regression that optimizes static and dynamic model parameters to minimize the difference between measured and modeled backscatter for a set of ∼ 300 000 pixels over a 4-year time period (2016–2019). The resulting SM product represents a percentage measure of the relative moisture saturation of the soil surface. The performance of the obtained SM time series was validated with in situ observations as well as compared to top-layer (0–7 cm) SM estimates from ERA5-Land. In addition, the spatial distribution of the resulting auxiliary model parameters (single-scattering albedo, soil scattering directionality) was analyzed with respect to the Climate Change Initiative (CCI) land cover (ESA, 2017) classifications to assess the physical plausibility of the resulting parametrization. The observed spatial pattern of the parameters indicate a close connection to the associated land cover, following some expected variations, e.g., higher single-scattering albedo over forested areas compared to croplands.

The RT1 high-resolution SM product over the Po River basin shows an overall good performance compared to ERA5-Land SM, with a median Pearson correlation coefficient of 0.55 for croplands and 0.65 over areas primarily covered by natural vegetation (i.e., tree, shrub, herbaceous cover). Validation was performed using in situ SM for the Oltrepo station (Bordoni et al., 2019) located in Canneto Pavese (PV, Italy), which resulted in a correlation coefficient of 0.58 (raw data) and 0.73 (with a 10-daily rolling mean) (Fig. 4). These results highlight the potential of Sentinel-1 observations for high-resolution SM retrievals and their use in applied science.

Figure 4 Time series of RT1 surface SM compared to in situ SM in Oltrepo. A 10-daily rolling mean is applied to the RT1 retrievals to reduce noise. The shading indicates the corresponding standard deviation (SD). Pearson correlation coefficient of 0.58 (raw data) and 0.73 (10-daily rolling mean).

2.4.4 Snow depth

Snow-depth data were obtained from the Sentinel-1-based product proposed by Lievens et al. (2019). The data product has a 1 km spatial resolution and daily granularity, and is available through the public repository of the C-SNOW project (https://ees.kuleuven.be/project/c-snow, last access: 28 July 2022). The mapping algorithm is based on a change-detection approach and has been validated across the mountain regions of the entire Northern Hemisphere.

For the scope of the present study, C-SNOW data during the period September 2016–April 2020 were evaluated with 172 ultrasonic snow-depth sensors across the Po River basin (Fig. 5a). Of the evaluation dataset, 77 % is located in the range 1000–2500 m above sea level (a.s.l.) (Fig. 5b), a frequent condition in the Alps (Avanzi et al., 2021). Observed snow-depth data were processed by (1) setting to missing any negative value, (2) applying climatological thresholds for maximum and minimum snow depth to remove spikes, and (3) using a threshold on the 6 h moving coefficient of variation to detect periods with grass interference (Avanzi et al., 2014). Data were then aggregated at daily resolution, and C-SNOW data were extracted for the same locations and data range. The evaluation confirmed previous results by Lievens et al. (2019), with C-SNOW successfully reproducing the seasonality and magnitude of snow depth as measured by snow-depth sensors (Fig. 5c and d). Root mean square errors (RMSEs) ranged from less than 20 cm below 1000 m a.s.l. to 60 cm or more above 2000 m a.s.l., though with no significant trend in the bias versus the elevation.

Figure 5 Evaluation of satellite based C-SNOW snow-depth estimates. (a) Location of the 172 sensors across the Po River basin, and (b) their elevation distribution. (c) Comparison between the interquartile range of C-SNOW and in situ measurements at the 172 sensors. (d) Comparison between daily C-SNOW estimates and ground-based snow-depth measurements for all sites.

2.4.5 River discharge

River discharge time series from satellite remote sensing are estimated by integrating data from two sensors: altimeter and multispectral. Traditionally defined as the product of cross-sectional river flow area and velocity, river discharge is calculated by assuming that the satellite sensors measure the two quantities (Tarpanelli et al., 2015). Specifically, once the cross-sectional geometry is known, flow area is calculated as a function of the water height derived from satellite altimetry (Abdalla et al., 2021), while flow velocity, usually measured through in situ instruments (current meter, acoustic Doppler current profiler, velocimeter), is linked to the reflectance measured by the near-infrared signal of the multispectral sensor (Tarpanelli et al., 2013), relying on the reflectance ratio between a dry (C) calibration pixel and the corresponding wet (M) measurement pixel.

Multi-mission satellite altimetry data coming from Saral/Altika, Cryosat-2, and Sentinel-3A and 3B are used to derive densified water level time series (Zakharova et al., 2020) at five stations along the main reach of the Po River basin named Piacenza, Cremona, Borgoforte, Sermide, and Pontelagoscuro (i.e., virtual stations in Fig. 1). At these stations, the multi-mission reflectance was extracted from the MODIS (Aqua and Terra), OLCI (Sentinel-3A) and MSI (Sentinel-2) sensors following the methods shown in Tarpanelli et al. (2020). Here, river discharge (Q) is estimated as the product of flow velocity (Tarpanelli et al., 2020) and flow area, as a function of altimetry-derived water height (H) (Tarpanelli et al., 2015):

$$\begin{array}{}\text{(3)}& Q=\mathit{\alpha }{\left(H\right)}^{\mathit{\beta }}{\left(C/M\right)}^{\mathit{\gamma }},\end{array}$$

where the parameters α, β, γ were calibrated using observed discharges at the five stations. The resulting time series for each station are illustrated in Fig. 6 against the in situ observations recorded at the gauged stations. Performance metrics (Supplement, Table S1) show skillful performance of the method in representing the observed daily discharges at the five stations, with average Nash–Sutcliffe (NS) of 0.81, KGE of 0.88 and relative RMSE (rRMSE) of 26 %.

Figure 6 Comparison between discharges simulated by the multi-mission approach versus observations at five gauging stations in the Po River basin in terms of time series (left column) and scatter plot (right column): (a, f) Piacenza, (b, g) Cremona, (c, h) Borgoforte, (d, i) Sermide and (e, l) Pontelagoscuro.

3 Methods

3.1 Hydrological modeling

Continuum (Silvestro et al., 2013) is a distributed hydrological model relying on a morphological approach that is based on the identification of the drainage network components (Giannoni et al., 2000). It is a tradeoff between empirical and physically-based models, reproducing all the main hydrological processes by relying on parameterization. The physical description of the hydrological processes is comparatively simple, resulting in high computational efficiency, yet generally skillful performance (Silvestro et al., 2013). Continuum reproduces the spatiotemporal evolution of runoff, SM, energy fluxes, surface soil temperature, snow accumulation and melting. Evaporation is estimated through a bulk formulation by solving the mass and energy balance as described in Silvestro et al. (2013) and related appendix, though it can also be provided as input variable. Deep flow and water table evolution are based on the Darcy equation, where each cell drains towards the neighboring cells following the 2D water table gradient and their hydraulic head, while a distributed interaction between water table and soil surface is represented through parameterization. A force-restore equation (Dickinson, 1988) is used to model the surface energy balance and enables the estimation of land surface temperature.

To simulate the cryospheric processes, we used S3M version 5 (Avanzi et al., 2022) – a one-layer snow model accounting for precipitation-phase partitioning, snowpack accumulation and melt, snow rheology and hydraulics, as well as glacier melt (Terzago et al., 2020; Avanzi et al., 2022). With its hybrid approach to snowmelt, which decouples the radiation- and temperature-driven contributions, S3M combines a parsimonious formulation with a substantial physical realism. For this work, S3M and Continuum were set up and run across the entire Po River basin (drainage area of 74 000 km²), with a constant grid spacing of 1 km and time resolution of 1 h.

3.2 Model calibration

To improve the representation of the hydrological states, Continuum was calibrated in the focus region using discharge data as benchmark. We deployed a multi-site calibration procedure that iteratively searches the model parameterization that best matches the available discharge observations over the calibration period at the 22 considered calibration stations (Fig. 1), through minimization of a cost function. Hydrological simulations run for the model calibration cover the 2 years starting on 1 January 2018, while the calibration period starts on 1 July 2018, leaving out the initial 6 months for model warm-up. The calibration tool perturbs six scalar parameters related to four physical hydrological features: infiltration velocity at saturation (cf), field capacity (ct), curve number (CN), and water sources (ws).

While the calibrated value of ws is a constant for the entire region of interest, the calibration of ct, cf, and CN consists of a rescaling of their default maps to the best value, thus preserving their spatial pattern, which depends on geographic spatial datasets of soil characteristics and land cover. The cost function, based on the Kling–Gupta efficiency (KGE; Gupta et al., 2009), computes an error between the duration curves at each percentile, weighted with the logarithm of the upstream area, to give higher weight to the downstream stations without neglecting the contribution of the most upstream ones.

The calibration procedure was performed through the implementation of a parallel search algorithm. The algorithm performs an iterative exploration of the 6D parameter space; the exploration starts with N= 20 initial values sampled with a Gaussian Latin hypercube approach. For each of these N parameter sets, a hydrological simulation is performed over the calibration period, and the cost function is computed to map the error hypersurface. The point that minimizes the cost function is used as the center of the following iteration, until the algorithm converges to an optimal solution.

3.3 Data assimilation of satellite snow and SM products

Satellite-derived SM from the Sentinel-1 RT1 product was assimilated into the Continuum model through a nudging technique (Stauffer and Seaman, 1990; Lakshmivarahan and Lewis, 2013). The nudging scheme is a computationally inexpensive approach and is particularly suitable for applications in operational frameworks for flood predictions. The update is performed when the satellite data become available, on average once per day for SM, following the equation

$$\begin{array}{}\text{(4)}& X_{\mathrm{MOD}}^{+}\left(t\right)=X_{\mathrm{MOD}}^{-}\left(t\right)+G\left[X_{\mathrm{OBS}}\left(t\right)-X_{\mathrm{MOD}}^{-}\left(t\right)\right],\end{array}$$

where $X_{\mathrm{MOD}}^{+}$ represents the updated modeled variable, $X_{\mathrm{MOD}}^{-}$ is the prior modeled value, X_OBS is the observation, and G is the kernel function. Thus, the correction term represents the difference between observed (X_OBS) and modeled variable multiplied by G that takes into account the uncertainties of both model and satellite observations. In this application we used a constant value of G= 0.45, following the recommendations by Laiolo et al. (2016), who estimated optimal G values from a test on four different satellite-derived SM products. In addition, we used G= 0 in areas with low Pearson correlation coefficient (r < 0.7) between satellite-derived and modeled SM in the simulation period.

The assimilation of satellite-derived C-SNOW maps into S3M was performed using the same approach and assuming G= 1 to mimic direct insertion. The C-SNOW maps provide snow depths, while S3M supports assimilation in the form of snow water equivalent (SWE), which is a more suitable variable to assimilate for controlling the water balance. Thus, snow depths from C-SNOW were converted into SWE using simulated snow-density values (see Avanzi et al., 2022). Along with snow-depth information, we rely on C-SNOW to determine snow-covered and snow-free areas, and then assimilated this information into S3M to clip modeled snow cover according to the satellite information. More information on the theoretical background of SWE assimilation in S3M can be found in Avanzi et al. (2022).

4 Results

4.1 Baseline run

The hydrological model Continuum was first calibrated using conventional meteorological data and observed discharges at the 22 calibration stations described above. The calibrated setup was then run over the years 2016–2019 to produce a baseline simulation for 2017–2019, excluding 2016 as model warm-up. Average evaporation in 2017–2019 computed by Continuum for the Po River basin is 950 mm yr⁻¹ and results are 21 % less than the GLEAM average of 1200 mm yr⁻¹ in the same time span. A comparison of simulated versus observed hourly river discharges is shown in Fig. 7 for five sample stations, while six performance metrics are shown for all 27 discharge stations in Fig. 8 and in Table S2 (see Supplement). Dimensionless scores, including KGE and its three decomposition terms, i.e., correlation (r), bias rate, and coefficient of variation rate (CV rate), increase on average with the upstream area. Note that all four scores have optimum value at 1. The mean KGE across all the stations (KGE = 0.51) rises to 0.63 and 0.70 for basins larger than 1000 and 10 000 km², respectively. Similar trends versus the same classes of upstream area are found in the mean r (0.75, 0.86, 0.88), while bias rate (0.98, 0.99, 0.94) and CV rate (0.89, 0.94, 1.13) are slightly deteriorated for basins larger than 10 000 km².

Figure 7 Observed versus simulated (baseline) discharge (m³ s⁻¹) for the years 2017–2019 at five river gauging stations.

Differences in the mean KGE, r and CV rate between validation and calibration stations are not statistically significant in a two-sample t-test for the mean. Only the mean bias rate of the two samples is statistically different at 5 % significance level, with validation stations having an average 30 % negative bias in comparison to an average 5 % positive bias of the calibration stations.

Figure 8 Skills of the baseline run versus upstream area at the 27 measurement stations. Dashed lines indicate the optimum value of each score.

4.2 Model runs with satellite input

In a second phase, we performed four hydrological simulations. Each simulation is based on the configuration and input data of the baseline run and by replacing in turn one input dataset with one of four satellite products described in Sect. 2.3: (1) precipitation from SM2RAIN; (2) evaporation from GLEAM; data assimilation of (3) SM from RT1, and (4) snow depths from C-SNOW. Two additional configurations were run including multiple satellite-based data sources: (5) all four satellite Earth observation datasets, hereafter referred to as EO, and (6) a combination of the satellite precipitation and evaporation, referred to as SM2RAIN + GLEAM. The spatial distribution of the performance of the six model simulations at the 27 river gauges is shown as maps of KGE (Fig. 9) and its three decomposition terms (see Supplement). Further, boxplots of KGE of the six experiments and comparison with the baseline run are shown in Fig. 10.

Figure 9 Spatial distribution of the Kling–Gupta efficiency (KGE) of the six model runs driven by the four input satellite products versus observed discharges at the measurement stations. Validation stations are marked with a bold circle. Multi-product experiments are in row 1, while single-product experiments are in row 2 (forcing input) and 3 (data assimilation input).

Figure 10 Boxplots comparing the KGE of simulated river discharges for all the considered experiments versus observations at the calibration and validation stations. The no-skill line at 1–2${}^{\mathrm{1}/\mathrm{2}}$ is indicated with a solid horizontal line. In the three rightmost columns, PE stands for precipitation and evaporation, while Q stands for discharge.

Results denote a generally skillful reconstruction of river discharges for all experiments, with mean KGE at the 27 stations ranging between 0.13 (SM2RAIN + GLEAM) and 0.53 (C-SNOW), all well above the no-skill threshold of KGE₀= 1–2${}^{\mathrm{1}/\mathrm{2}}\cong -\mathrm{0.41}$ (see Knoben et al., 2019). Simulations including C-SNOW and GLEAM perform on average better than the baseline run, with mean improvements in KGE of 0.02 and 0.01 (+4 % and +2 %), respectively. The largest differences in the overall performances are due to the wide range of the mean bias across the six simulations, with the largest bias rates for SM2RAIN + GLEAM (1.58) and EO (0.69), and the lowest bias rate for GLEAM (1.02) and C-SNOW (0.97), both improving that of the baseline run (0.95). On the other hand, average correlations across the six experiments fall in a much narrower interval, ranging between 0.61 for EO and 0.75 for both C-SNOW and the baseline run. Running the model with all EO data produces on average a 28 % deterioration of the mean performance (KGE = 0.37), though it surprisingly generates the best performance at the five validation stations (KGE = 0.54) among all simulations (Fig. 10).

The six simulations forced by satellite products were then compared to the baseline run to detect similarities and deviations in the entire simulation domain, including areas where no observations are available. To reduce the correlation effects along the river network, we consider only one value per simulated river reach, located just upstream of each confluence. The use of RT1 and GLEAM does not result in significant spatial differences with respect to the baseline (Fig. 11). As expected, using C-SNOW results in differences mainly in alpine areas, especially in Ticino (Switzerland), where the MCM dataset used in the baseline run is known to underestimate precipitation rates due to the lack of ground measurements outside the Italian territory. Larger deviations are visible in the runs including SM2RAIN, particularly in the upper Po River basin in the west and the upper Adda River in the north, confirming the stronger sensitivity of river discharge to precipitation dynamics.

Figure 11 Spatial distribution of the KGE of discharges of the six model runs driven by the four input satellite products versus the baseline run, at each modeled river reach. Multi-product experiments are in row 1, while single-product experiments are in rows 2 (forcing input) and 3 (data assimilation input).

Figure 12 Comparison of observed and simulated hydrographs (m³ s⁻¹) for the events between October and December 2019 at two validation stations: Candoglia (a) and Piacenza (b), together with KGE of calculated versus the observed discharges for the same three months.

Figure 12 shows a comparison of the six simulations forced by satellite products, the baseline run, and observed discharges at two validation stations, for a series of moderate to high intensity events which hit a large portion of the Po River basin in the autumn of 2019. The second of the three main events, during the second half of November, caused the exceedance of the maximum alert level and widespread flooding in several river sections in the main reach of the Po River basin across the Lombardy and Emilia-Romagna regions, including the area of Piacenza (Fig. 12, bottom). In Piacenza, all model simulations performed reasonably well, with maximum error on the peak discharge below 20 %. The best performances over the three months are found in the baseline run and in the two runs with data assimilation (RT1 and C-SNOW), all three with KGE = 0.89. Lower performances are produced by the three runs forced by SM2RAIN, mainly due to an overestimation of the first event in late October 2019. At the Candoglia station, results show an opposite pattern, with the best performances by SM2RAIN and SM2RAIN + GLEAM, both with KGE = 0.74 over the three months. This mildly improves on the performance of the baseline run (KGE = 0.71).

4.3 Sensitivity of satellite data to three model parameterizations

A subsequent experiment investigated the performance of the hydrological model in reproducing discharges at the 27 river gauges, by forcing it with the satellite datasets SM2RAIN and GLEAM. We compared the results of three model runs over 2017–2019 in detail, using three different model parameterizations obtained through dedicated calibrations (over 2018–2019), derived by applying the steps described in Sect. 3.2 to different configurations of input and benchmark discharges:

1. The first is the simulation SM2RAIN + GLEAM described in Sect. 4.2, i.e., run with the model parameters obtained by calibrating with conventional ground observations (interpolated measurements and MCM precipitation) and optimizing the objective function using observed discharge at the 22 calibration stations as benchmark (obs PE, obs Q in Fig. 10).

2. The simulation SM2RAIN + GLEAM run on a model calibration forced by the same satellite datasets SM2RAIN and GLEAM as input and optimizing the objective function using observed discharge at the 22 calibration stations as benchmark (EO PE, obs Q in Fig. 10).

3. The simulation SM2RAIN + GLEAM run on a model calibration forced by the satellite datasets SM2RAIN and GLEAM as input and optimizing the objective function using satellite-derived discharge estimates at the five virtual stations (see Sect. 2.3.5) as a benchmark (EO PE, EO Q in Fig. 10).

It is worth noting that SM and snow-depth data were not used in this experiment because they are not model input variables but rather assimilation variables, hence the calibration procedure described in Sect. 3.2 would not be directly applicable. Results from simulation no. 2 forced by the same SM2RAIN and GLEAM used in the calibration shows the lowest performance among the three (mean KGE = 0.07 over all 27 stations). Simulation no. 3 (EO PE, EO Q) gives satisfactory performance (mean KGE = 0.10), relatively close to no. 1 (mean KGE = 0.13), despite relying largely on satellite data. Interestingly, the five validation stations on average outperform the set of calibration stations, with average KGE of 0.38, 0.30, and 0.29 for the three experiments. The performance of the three model runs versus the upstream area at the 27 stations (Fig. 13) shows a general improvement in the correlation with the upstream area, while for the other metrics trends are less clear. Simulation no. 3 shows reduced variability (CV rate), yet smaller absolute errors (RMSE and mean error in Fig. 13), also thanks to a calibration focused on the downstream virtual stations.

Figure 13 Skills of the run forced by satellite precipitation and evaporation (PE) versus upstream area at the 27 measurement stations. The three markers denote three calibrated parameter sets, obtained with different configurations of PE input and of benchmark discharge (Q). Conventional observational datasets are indicated with “obs”, while “EO” are the satellite-derived datasets. Dashed lines indicate the optimum value of each score.

5 Discussions

A critical evaluation of the results from the experiments performed can help to identify strengths and weaknesses, as well as measures that can be taken to maximize the benefits of satellite observations in Earth system modeling. Overall, hydrological simulations driven by satellite datasets produced encouraging results, with 95 % of KGE of the station–experiment combinations above the no-skill threshold (versus 100 % for the baseline run). The remaining 5 % of combinations with KGE below the no-skill threshold occur in just 3 of 27 stations and only in model configurations including SM2RAIN. Generally, the precipitation dataset is found to have the largest weight on the resulting model performance, with standard deviation (SD) of changes in KGE versus the baseline simulation (SD${}_{\mathrm{\Delta }\mathrm{KGE},\mathrm{SM}\mathrm{2}\mathrm{RAIN}}=$ 0.37) being more than twice that of all the other satellite-driven configurations (SD${}_{\mathrm{\Delta }\mathrm{KGE},\mathrm{RT}\mathrm{1}}=$ 0.16, SD${}_{\mathrm{\Delta }\mathrm{KGE},\mathrm{GLEAM}}=$ 0.09, SD${}_{\mathrm{\Delta }\mathrm{KGE},\mathrm{C}-\mathrm{SNOW}}=$ 0.06). In other words, the simulation performance shows strongest sensitivity to the precipitation forcing, which in fact leads to the largest deteriorations compared to the baseline run. Additionally, some of the largest improvements in KGE, up to ΔKGE${}_{\mathrm{MAX},\mathrm{SM}\mathrm{2}\mathrm{RAIN}}=$ 0.29, are well above all the improvements produced by GLEAM (ΔKGE${}_{\mathrm{MAX},\mathrm{GLEAM}}=$ 0.17) and C-SNOW (ΔKGE_MAX,C-SNOW= 0.12) at any single station. This result is largely in agreement with previous findings (e.g., Jones et al., 2006; Sperna Weiland et al., 2015) and highlights the importance of advances in satellite precipitation estimation for hydrological applications. Qi et al. (2016) showed that model performance can also be impacted by model-precipitation product interactions, though this can partly be mitigated by dedicated model calibrations for each combination of input products. The high-resolution version of SM2RAIN used in this work leads to hydrological performance comparable to that of the best non-gauge-corrected satellite products found in the literature (Camici et al., 2018; Amorim et al., 2020), and local results are better than those obtained with previous coarser resolution versions (see e.g., Beck et al., 2017; Tang et al., 2020). These works also show that satellite-based precipitation datasets that are bias-corrected with ground observations further improve the overall quality, including the performance of hydrological modeling.

With regard to the precipitation forcing, one must also note that the MCM dataset used in the baseline represents a benchmark that is particularly difficult to overcome. The high station density and the merging with the Italian radar composite make MCM a high-quality and detailed product both spatially and temporally. Yet only a few world areas can rely on seamless and nearly unbiased gauge-radar products, while satellite datasets remain prime candidates in ungauged regions, especially for real-time applications, thanks to key features such as extended coverage, high resolution, short latency, and spatial consistency. In addition, satellite datasets are unaffected by country borders, making them suitable for applications in transboundary river basins, especially in countries where data-sharing agreements are not easily implemented.

In contrast, GLEAM and C-SNOW consistently produced moderate improvements, though on a larger number of river sections, with only a minority of stations where skills deteriorated in comparison to the baseline run. Finally, the assimilation of RT1 SM shows contrasting behavior. On the one hand, it deteriorated KGE values throughout most of the stations in the main reach of the Po River basin due to a general increasing negative bias. On the other hand, it shows general benefits in small-size upstream catchments and notably the best improvement in KGE (ΔKGE${}_{\mathrm{MAX},\mathrm{RT}\mathrm{1}}=$ 0.41) among all 216 station–experiment combinations for the Trebbia River at Valsigiara. Our findings confirm the challenges related to implementing a semi-automated assimilation of satellite SM already pointed out in previous research (Laiolo et al., 2016; Wanders et al., 2014), where a range of factors affect and often decrease the assimilation performance, including the presence of complex topography, snow cover, frozen soil, urban areas, as well as differences between modeled and actual vegetation cover and LAI.

A final comment goes to the surprisingly high skills of hydrological simulations at the five validation stations, which on average exceed those at the calibration stations in five of nine experiments (see Fig. 10). The multi-site calibration strategy is designed to find an optimal parameter set for the entire domain, thus reducing the effect of highly variable model performance which is typical of cascading calibrations (e.g., Alfieri et al., 2020). All results are then compared at the calibration and validation stations for the same period 2017–2019. This is twice the duration of the calibration period, implicitly adding a validation component at the calibration stations as well. Higher performance at the validation stations seem to be particularly evident in simulations forced by SM2RAIN, though a connection between these facts is not known and it may simply be related to spatial differences in the skills of the satellite-based precipitation forcing in the sub-catchments where validation stations are located. A noteworthy case is that of the validation station Toce River at Candoglia, in the northwestern part of the Po River basin. It is influenced by a large number of reservoirs upstream and Lake Maggiore, located just downstream, hugely smoothens its runoff characteristics from the rest of the river network. This makes the sub-basin almost disconnected to the rest of the Po River basin. Notwithstanding, simulation performance at Candoglia is higher than those of the calibration stations in all experiments but one (RT1), with the case of SM2RAIN scoring a KGE = 0.74, hence 0.22 points higher than the average calibration KGE among all stations. Given the number of experiments presented, focused on the role of different input data and model parameterization, results are only shown through overall statistics of each model run. Future work will investigate detailed model behavior over specific hydrological processes, regimes, seasonality, and quantiles of the flow duration curve, to better disentangle strengths and weaknesses of the considered satellite products under specific hydrological conditions.

6 Conclusions

This research explored the impact of five high-resolution satellite products in distributed hydrological modeling. In a set of experiments we tested the use of satellite-based precipitation and evaporation as forcing input, data assimilation of satellite SM and snow depth, and satellite river discharge estimates as benchmark for model calibration. We found skillful performances for all simulations including satellite-derived products, with GLEAM evaporation and C-SNOW snow depth respectively yielding an average 2 % and 4 % improvements over a baseline run driven by high-quality ground-based datasets. The skills of model runs including EO data showed considerable variability in space and time. In addition, we found skillful results in a model calibration relying heavily on satellite products, both with regard to forcing input and benchmark discharge. This heralds the use of hydrological models that fully rely on satellite data as an appealing solution for large-scale applications and for regions where ground-based observations are not available, particularly in near real time.