Operational Implementation of Satellite-Rain Gauge Data Merging for Hydrological Modeling

Systems exposed to hydroclimatic variability, such as the integrated electric system in Uruguay, increasingly require real-time multiscale information to optimize management. Monitoring of the precipitation field is key to inform the future hydroelectric energy availability. We present an operational implementation of an algorithm that merges satellite precipitation estimates with rain gauge data, based on a 3-step technique: (i) Regression of station data on the satellite estimate using a Generalized Linear Model; (ii) Interpolation of the regression residuals at station locations to the entire grid using Ordinary Kriging and (iii) Application of a rain/no rain mask. The operational implementation follows five steps: (i) Data download and daily accumulation; (ii) Data quality control; (iii) Merging technique; (iv) Hydrological modeling and (v) Electricity-system simulation. The hydrological modeling is carried with the GR4J rainfall-runoff model applied to 17 sub-catchments of the G. Terra basin with routing up to the reservoir. The implementation became operational at the Electricity Market Administration (ADME) on June 2020. The performance of the merged precipitation estimate was evaluated through comparison with an independent, dense and uniformly distributed rain gauge network using several relevant statistics. Further validation is presented comparing the simulated inflow to the estimate derived from a reservoir mass budget. Results confirm that the estimation that incorporates the satellite information in addition to the surface observations has a higher performance than the one that only uses rain gauge data, both in the rainfall statistical evaluation and hydrological simulation.


Introduction
The renewable contribution of the electric energy matrix in Uruguay has been increasing steadily during the last decades, with hydroelectric, wind and solar components that have different inherent variability and predictability. This poses both a challenge and an opportunity to optimize planning at different embedded timescales and, ultimately, dispatch. The interconnected Electric System Simulator (SimSEE [1]) is used for these purposes [2], from the management of the spot market to long-term analysis of the evolution of the generation capacity, with intermediate seasonal and nested weekly planning. Particularly, in the case of the hydroelectric generation, we used a coupled hydrological and electric system modeling approach in order to generate and process a hydrological ensemble forecast for the largest reservoir of the system [3]. The ability to forecast the hydrological inflow contributes to the optimal use of each energy source, with the corresponding economic and environmental benefits.
In most applications that use operational hydrological models, the spatial and temporal variability of precipitation constitutes one of the dominant factors and with greater associated uncertainty. In this context, remote sensing products are ideal instruments to through comparison with an independent historic rain gauge dataset. Furthermore, a hydrological application was implemented using the GR4J model [24] at daily time step and compared to the estimated "theoretical" inflow to the hydroelectric reservoir. We expect that this work will contribute to the understanding of the reliability of the latest NRT satellite-based precipitation products and provide a reference for their applications in operational hydrological simulation and water resource management.

Study Area
The upper Rio Negro basin, in northeastern Uruguay, has a surface area of about 40,000 km 2 ( Table 1), taking Gabriel Terra hydroelectric plant (G. Terra) as its closure point. Downstream G. Terra in the Rio Negro, we find Baygorria and Constitución hydroelectric plants. The binational (Argentina-Uruguay) hydroelectric plant Salto Grande, in the Uruguay River, completes the total hydroelectric capacity that collectively represents a third of the current installed power in the country's electric system and contribute with more than the 50% of the mean total generated electricity [25] with large interannual variability. In the Uruguayan system, the main storable resource is the water in the hydropower reservoirs, particularly in G. Terra, since it has the highest storage capacity ( Table 2). Considering the growing pressure for water demand from both agricultural and forestry expansion, together with the continuous increase of electric energy demand, this highlights the need for adequate tools for the management of water resources in the Rio Negro basin [26].  Figure 1 shows the location of existing hydroelectric plants (black triangle) and the delimitation of G. Terra basin. We also included the location of the rain gauges available at NRT (red square), used for the operational implementation of the merging approach and the historic rain gauge data (blue dot) used for validation purposes. Both datasets are presented in Section 3.

Rain Gauge Data Available in Near-Real-Time
Precipitation data available in NRT, used for the operational implementation of the merging technique, comes from the public electric utility (UTE) network. After a data quality control, we selected 19 automatic stations within the Rio Negro basin. The data quality control included the identification of missing data and outlier values, the implementation of plausibility checks based on Scherrer et al. (2011) [27], as well as the evaluation of the accumulated and mean annual rainfall, the average number of wet days (having nonzero rainfall) and the length of the longest dry spell. The period analyzed is 31 January 2010 to 31 May 2020. Daily rainfall totals are taken at 1000 UTC. Figure 1 shows the location of the selected pluviometric stations (symbolized with red squares). It highlights that they are not uniformly distributed, which surely influences the performance of the merging technique to be implemented.

Historic Rain Gauge Data for Validation Purposes
The historic reference data used to evaluate the proposed methodology comes from a relatively dense and uniformly distributed network of 95 stations provided by UTE, the National Institute of Meteorology (INUMET), the National Institute of Agricultural Research (INIA) and the National Water Agency (ANA) of Brazil. Figure 1 presents the spatial distribution of these stations (blue dots) in comparison with the location of the automatic stations selected for the operational implementation of the merging approach (red squares).
The validation period is 1 February 2017 to 31 May 2020, during which the satellite rainfall estimates selected (presented in Section 3.3) are also available. Daily rainfall totals are taken at 1000 UTC.

Satellite Rainfall Estimates
In view of the new generation of global precipitation satellite products, which integrate multiple platforms and previously existing algorithms, with high spatial and temporal resolution and better performance than the predecessor products [17,28], the following products were selected: GSMaP: Global Satellite Mapping of Precipitation [10,11] of the Japan Aerospace Exploration Agency (JAXA), version GSMaP-Gauge-NRT v7 (https://sharaku.eorc.jaxa.jp/ GSMaP/). IMERG: Integrated Multi-satellitE Retrievals for GPM, Global Precipitation Measurement [12] of the National Aeronautics and Space Administration (NASA), version Level 3 V06, NRT Late Run (https://pmm.nasa.gov/data-access/downloads/gpm).
We also took into account in the selection the data latency and accuracy of the different available products. Table 3 presents the main characteristics of the selected datasets (spatial and temporal resolution, data latency and period of availability). Although both products are available on an hourly frequency, the present study is limited to daily precipitation totals since this is the information needed for the hydrological modeling (presented in Section 4.2). As a first exploration to evaluate the satellite-rainfall estimates at daily time step, the root mean squared error (RMSE) and the probability of detection (POD) for a precipitation threshold of 5 mm were calculated for both, GSMaP and IMERG, against the observed records. To this end, we only considered those grid boxes of the satellite grid that contained at least one gauge observation for the specific day (collocated gauge-satellite data pairs). Table 4 presents the results obtained for the period 1 February 2017 to 31 May 2020. Additionally, the coverage of each satellite product for the entire period is included, as expressed as the percentage of pixels × day of available data. It shows that, in both cases, the satellite estimates present a satisfactory performance. Both products have a very good coverage in the analyzed period (close to 100%).
However, as mentioned earlier, previous experience with this type of products in the study area [22,23] confirms the need to implement a bias removal scheme based on available surface observations prior to any application.
As an example, Figure 2 shows the comparison of the rain gauge observations and satellite-rainfall estimates for a particular day (15 December 2019).

Other Data
The following datasets are used for the hydrological modeling (presented in Section 4.2): Precipitation forecast: a 14-day ensemble precipitation forecast is obtained from the Global Ensemble Forecast System (GEFS v11.0) produced by the National Centers for Environmental Prediction (NCEP-NOAA). The ensemble is composed of the control run and 20 perturbed members and has a spatial resolution of 1 • × 1 • [29].
Potential evapotranspiration (PET): the mean annual cycle of PET was calculated from the records of 9 meteorological stations belonging to INUMET and INIA for the period 1991-2015, using the Penman-Monteith method.
Amount of water storage capacity (SC) in the soils present in the G. Terra basin: the SC for each soil type was obtained from the CONEAT soil map at scale 1:40.000 [30] of the Office of Natural Resources of the Ministry of Livestock, Agriculture and Fisheries of Uruguay (DGRN-MGAP). Then, it is weighted by area to obtain a representative value for each sub-basin. A digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM-NASA) with a resolution of 90 m was used to perform watersheds delineation and characterization.
Additionally, for the evaluation of the hydrological model, we used the daily series of estimated inflow to G. Terra reservoir provided by UTE (represented with grey dots in Figure 3). This series is called "theoretical" since it consists of an estimation based on a water balance in the reservoir and it is not a direct observation. Specifically, the estimated "theoretical" inflow is obtained (indirectly) from the water surface elevation at the dam and the turbinated and discharged flows. Therefore, this estimation is sensitive to the representation of the reservoir and the effect of the wind on its surface. Indeed, Figure 3 shows negative inflow values, which may be due to the compensation effect of excessively high values particular of the methodology (possibly associated with the action of the wind in the reservoir). Figure 3 also includes the time series of the 7-days filtered estimated "theoretical" inflow (blue line), considering that the model is used as a tool to support the decision-making of the weekly dispatch.

Merging Approach
The satellite-rain gauge data merging technique considered is based on the universal model of spatial variation [31,32]. As one of the hybrid geostatistical models, Regression Kriging (RK) is a spatial interpolation technique that combines a deterministic model (regression) with a statistical model (Ordinary Kriging of the regression residuals). It uses a deterministic model to estimate a value of the variable (precipitation) by using actual ground measurements to calibrate a model for the satellite estimates and then refines the estimate analyzing the residuals for spatial correlation; finally, it combines the statistical fitting and deterministic modeling [33].
The 3-step proposed model is summarized as follows ( Figure 4). Regression of the station data on the satellite data using a Generalized Linear Model (GLM). A GLM model is implemented in order to fit the satellite estimates to the rainfall observations at station locations. For the GLM, we use a spatially correlated residual structure that is fit to the available data. For each day, we calculate both exponential and spherical spatial correlations and choose the one with the highest Akaike information criterion (AIC). Several regressors alternatives were tested, including both satellite products and each one separately. Based on the performance statistics obtained (not shown), we decided to use, for each day the individual satellite product, either IMERG or GSMaP, with the highest Pearson correlation between the rain gauge observations and the collocated satellite values.
Interpolation of the regression residuals at station locations to the entire grid using Ordinary Kriging. Once the regressed satellite estimation is obtained, we calculate the error (residual) between it and the observations at the station locations. Then, the interpolation of the regression residuals to the entire grid is done through Ordinary Kriging [34], which exploits the spatial correlation in the residuals and this is added to the regressed satellite estimation in order to obtain the "unmasked" merged product.
Application of a rain/no rain mask (RNR mask). We apply an RNR mask to the merged product to prevent overestimation of the occurrence of rainfall in the interpolated field. The mask is obtained using the same merged precipitation estimate (RK) technique but switching the target observations to binary rain/no rain observations. Satellite estimates are used as regressors to forecast this binary field with the same RK technique described in 1 and 2. We use a threshold of 0.3 in the output of RK, a continuous field, to delimitate rainy region for the mask. Finally, the unmasked product is multiplied by the RNR mask to obtain the final masked merged product. Figure 5 shows an example of the application of the RNR mask for a given day (6 July 2017). The middle column corresponds to the unmasked OK and RK products while the rightmost column shows the masked versions. As can be seen, there is a large purple region of slightly positive values (zeros are transparent) in the unmasked product, while in the masked products this region is forced to zero. The merging algorithm in this study was written in R and is available on GitHub [35,36].

Rainfall-Runoff Model
The G. Terra basin (39,500 km 2 ) was discretized into 17 sub-basins with areas smaller than 7000 km 2 . Figure 6 presents the delimited catchments and their characteristics including: basin area (km 2 ), slope (%), concentration-time (Tc) (h) and SC (mm). To simulate the hydrological inflows to G. Terra reservoir we use a daily hydrological model (GR4J) coupled with a hydrological transit model (Muskingum). The GR4J model is a daily lumped four-parameter rainfall-runoff model developed by Perrin et al. (2003) [24]. The Muskingum model [37] is a two-parameter hydrologic flood routing method, based on the storage continuity equation.
In a previous study, Narbondo et al. (2020) [38] present a successful application of the GR4J daily rainfall-runoff model at 13 watersheds of Uruguay. They proposed an improved regionalization approach to predict runoff time series in ungauged catchments at country scale. Particularly, they found the optimal set of parameters of the GR4J model and, in addition, they found the relationships between them and watershed-physiographic factors. Table 5 shows the description of the "GR4J-Muskingum" model parameters and the values adopted in each case following these recommendations.

Goodness-of-Fit Indicators
The performance of the merged precipitation estimate (RK) was statistically evaluated through comparison with an independent rain gauge network, relatively dense and uniformly distributed (referred to as "historic data" in Figure 1). We also included in the evaluation the estimation based on the Ordinary Kriging interpolation from NRT rain gauge observations (OK) at the same grid as the satellite data (0.1 • × 0.1 • ), which serves as the baseline for comparison with the merged product. Both estimates (RK and OK) were compared with the rain gauge observations belonging to the historic reference dataset. The performance statistics used for the comparison are the mean error (ME), the RMSE, the frequency bias (FBS), the POD and the false alarm ratio (FAR) for a precipitation threshold of 5 mm [39,40].
Furthermore, several verification indices were used to quantitatively assess the hydrological utility of the precipitation estimates based on the estimated "theoretical" inflow to G. Terra reservoir (Figure 3), including the difference of total accumulated inflow (∆V), the Nash-Sutcliffe efficiency (NSE), the Kling-Gupta efficiency (KGE), the coefficient of determination (R 2 ) and RMSE [41]. Additionally, we also conducted a first-level catchment water balance using the runoff ratio (RR), defined as the ratio of the precipitation that contributes to runoff [20]. The RR values calculated using the different outputs from both estimates (RK and OK) were compared to known values from the literature [42].
In all cases, the period analyzed is from 1 February 2017 to 31 May 2020.

Operational Implementation
The 5-step operational implementation of the coupled hydrological and electric system modeling approach is presented next (Figure 7). Data download and daily accumulation. The required precipitation input data are adequately collected: records of NRT stations, GSMaP-NRT, IMERG-NRT Late Run and GEFS ensemble forecast. Daily rainfall totals are accumulated at 1000 UTC.
Data quality control. Prior to the merging algorithm, a data quality control from both NRT rain gauges and satellite estimates is performed based on the Climate Data Tools (CDT-IRI) [43]. Data quality control focuses on outlier detection for the purpose of elimination of data contamination, including the implementation of spatial-plausibility checks based on Scherrer et al. (2011) [27]. The threshold values used in the controls were adjusted manually, looking to eliminate the most obvious suspicious values in the available historical data set.
Merging technique. The satellite-rain gauge data merging technique is implemented in order to obtain the RK precipitation estimate over the Rio Negro basin.
Hydrological modeling. Based on the RK estimate and the GEFS precipitation ensemble forecast, the GR4J rainfall-runoff model is implemented at the 17 sub-catchments of the G. Terra basin. The runoff output is then routed along the river network using the Muskingum model to simulate the daily inflow ensemble to G. Terra reservoir.
Electricity-system simulation. The simulated inflow ensemble is integrated to the existing model of the interconnected electric system (SimSEE), particularly into the synthesizer model (CEGH), through biases and noise attenuators per time step adjusted through maximum likelihood [44].
The implemented model was integrated into SimSEE's on June 2020 and has since run under the responsibility of the Electricity Market Administration (ADME) of Uruguay. The application (called VATES) is continuously updating and executing a SimSEE Room with the representation of the Uruguayan generation system, in order to obtain the dispatch of the following seven days with hourly detail. The results and information relevant to the operation are published automatically on ADME's website [45]. They also provide the required statistical information for the design of exchange offers with neighboring countries and the energy spot market. Table 6 presents the comparison of the performance metrics for the OK (stations only) and the RK (merged product) precipitation estimates. The results obtained are global values, integrated both spatially (among the 95 stations in the Rio Negro basin) and temporally (averaged over the analyzed period). Overall, both estimates have a good performance but RK performs slightly better. This indicates an improvement in the accuracy of the precipitation estimation by the incorporation of satellite data.  Figure 8 shows the spatial distribution of RMSE at the reference data locations obtained with both estimates, OK and RK, averaged over the analyzed period. As can be seen, both maps exhibit similar distribution patterns but there are some differences on the border with Brazil, where there are practically no NRT stations (see Figure 1). The RK estimate in that region has a better performance than OK with differences in RMSE between 10% and 20%.

Rainfall Model Performance
As an example, Figure 9 compares the interpolated precipitation fields obtained with OK (stations only) and RK (stations and satellite) for a given day (15 December 2019). Black dots represent the NRT rainfall observations used for the interpolation. As expected, both maps present a similar spatial distribution of daily precipitation but, particularly in the region towards Rio Grande Do Sul (Brazil), the two estimates show significantly different values, with OK showing amounts of the order of 40 mm and RK of 25-30 mm. In this region, the OK estimate does not appear to be natural, with uniform high values in a smooth zone that gradually fades, as a result of the weighted sum estimation, rather than irregular zones with intensely high peaks as observed in the RK estimate. This highlights the advantage of satellite data in the representation of spatial rainfall variability, particularly in data-sparse regions, as is the case here.

Rainfall-Runoff Model Performance
In this section, the precipitation estimates are incorporated into the hydrological model that runs operationally in ADME (see Section 4.2) and the performance is assessed as compared to the estimated "theoretical" inflow to the hydroelectric reservoir ( Figure 3). The "GR4J-Muskingum" model was forced by the two precipitation estimates (OK and RK) using the model parameters presented in Table 5 to simulate the daily inflows to G. Terra reservoir for the period 1 February 2017 to 31 May 2020.
The simulated and "theoretical" hydrographs are shown in Figure 10 and the statistical comparisons are summarized in Tables 7 and 8. Considering that the model is used as a tool to support the decision-making of the weekly dispatch, we also included the comparison of the 7-days moving average inflows (Figure 11).   Figure 11. Comparison between theoretical and simulated (Ordinary Kriging and Regression Kriging) 7-days filtered inflows to G. Terra reservoir.
As shown in Figures 10 and 11, both simulations have generally good agreement with the "theoretical" streamflow both at daily and weekly time step, although overestimation and underestimation of the peaks are evident in some cases.
Both estimates present a slight underestimation of total accumulated inflow, with a difference of −5.3% and −2.5% for OK (station only) and RK (merged product) output respectively (Table 7); which correspond to a very good performance according to Moriasi et al. (2015) [41].
DINAGUA [42] reports an annual runoff ratio (RR) between 0.37 and 0.43. This reference value is very close to those obtained with the OK and RK estimates from both the estimated "theoretical" inflow (RR Est) and the simulated inflow (RR Sim) ( Table 7). On the one hand, when using the "theoretical" inflows to calculate the RR, we verified that the performance of the precipitation estimates is satisfactory. On the other hand, when considering the simulated inflows series, we confirmed that the hydrological model achieves a good representation of the rainfall-runoff transformation process, regardless of the precipitation estimate considered.
As context, we present the general performance ratings for the adopted statistics recommended by Moriasi et al. (2015) [41], simulated inflow using both estimates (RK and OK) have a satisfactory performance at daily time step (0.60 < R 2 ≤ 0.75 and 0.50 < NSE ≤ 0.70) and a good performance at weekly time step (0.75 < R 2 ≤ 0.85 and 0.70 < NSE ≤ 0.80). Furthermore, for all statistics considered, the RK simulated inflow has a better performance than the OK, for both daily and weekly time step.

Discussion
These results confirm that, as expected, the estimation that incorporates the satellite information in addition to the surface observations (RK) has a higher performance than the one that only incorporates the rain gauge data (OK), both in the rainfall statistical evaluation and hydrological simulation of the basin.
However, the magnitude of the improvement in the rainfall estimation is relatively small as expressed by the global indicators shown in Table 6, averaged both in time and space. Figures 8 and 9 already suggest that the magnitude of the original error with OK and improvement with RK, might be larger in the upper part of the basin, where the density of rain gauges is notably lower and, given the lack of stations on the other side of the water part, extrapolations are required to cover the basin. This is verified in Table 9 where we limit the RMSE indicator to the higher sub-basins (see Figure 6). We limited the analysis to RMSE because it is the most robust statistic, as well as the most relevant for the application and does not require the definition of a precipitation threshold like FBS, POD and FAR. Table 9 also includes the percentage of improvement achieved with RK, which increases from 3% for the global indicator to approximately 20% in the more poorly monitored border sub-basins. While synoptic frontal systems are prevalent through the year in the region and are responsible for most of the rainfall, convective scale storms become a relevant contributor to precipitation totals during the warm season. Of course, many times the latter are embedded in the former generating the multiscale structure of precipitations fields. However, it is well known that precipitation daily totals decorrelate with distance faster in the warm season as compared to the cold one [22]. This motivated an analysis of the seasonality of the improvement in skill when the satellite estimates are incorporated (RK). Table 9 shows the basin averaged RMSE limited to the warmest semester (October through March) and the peak of the warm season: December-January-February (DJF). Even with relatively high density of surface observations, as is the case on average in the region of study, the impact of incorporating satellite information increases as the precipitation field acquires larger amplitude in smaller scales during the warm season, from 3% up to 11%.
These analyses give an insight of the potential for improvement in skill that can be obtained with the merging methodology proposed as a function of rain gauge density and characteristic of the precipitation systems.

Summary and Conclusions
In this study, we developed and implemented a methodology that combines rain gauge observations and satellite-rainfall estimates at daily time step to improve the rainfall monitoring in NRT. The proposed methodology involves 3 steps: (1) regression of station data on the satellite estimate using a Generalized Linear Model, (2) interpolation of the regression residuals at station locations to the entire grid using Ordinary Kriging and (3) an application of a rain/no rain mask. The merged precipitation field thus obtained is then used in a hydrological modeling of the Rio Negro basin whose output is, in turn, coupled with an electric system modeling that guides planning and dispatch decisions for the following seven days.
The performance of the proposed merged precipitation estimate was statistically evaluated through comparison with an independent historic rain gauge dataset. The incorporation of satellite information enhances the representation of spatial variability, particularly in data-sparse regions with reductions in RMSE of up to 20%, although the overall improvement is statistically marginal.
As far as the operation of the energy system is concerned, it is the input to the reservoirs that most directly affect the electric system simulations and, in turn, management optimization. The GR4J hydrological model, with a daily time step, was implemented at 17 sub-catchments of the G. Terra basin with routing up to the reservoir. Model performance was assessed comparing model output to the estimated "theoretical" inflow to G. Terra computed from a mass budget to the reservoir and rendered satisfactory statistics: 0.60 < R 2 ≤ 0.75 and 0.50 < NSE ≤ 0.70. The estimation that incorporates the satellite information in addition to the surface observations has a higher performance, for all statistics considered, compared to the one that only incorporates the rain gauge data.
In an operational setting, simplicity and robustness of the implementation are as important as accuracy. All steps are currently implemented and run on a daily basis at the Electricity Market Administration (ADME): data download and quality control, merging algorithm, hydrological modeling and electric system simulation. The presented implementation improves the estimation of the precipitation field and carries that information all the way to the decision-making stage, with its corresponding socio-economic and environmental benefits.