Complementing ERA5 and E-OBS with high-resolution river discharge over Europe

The 0.5 ° resolution of many global observational or quasi-observational datasets is not sufﬁcient for the evaluation of current state-of-the-art regional climate models or the forcing of ocean model simulations over Europe. While higher resolved products are available for meteorological data, e.g. ERA5 reanalysis and the E-OBS vs 22 (EOBS22) datasets, they lack crucial information at the land-ocean boundary. ERA5 is frequently used to force regional climate models (RCMs) or ocean models and both datasets are commonly used as reference datasets for the evaluation of RCMs. Therefore, we extended both datasets with high-resolution river discharge for the period 1979—2018. On the one hand, our discharge data close the water cycle at the land-ocean interface so that the discharges can be used as lateral freshwater input for ocean models applied in the European region. On the other hand, the data can be used to identify trends in discharge that are induced by recent climate change as ERA5 and EOBS22 are rather independent datasets. The experimental setup to generate the discharges was chosen in a way that it could be easily adapted in a climate or Earth system modelling framework. Consequently, the recently developed 5 Min. horizontal resolution version of the hydrological discharge (HD) model was used to simulate discharge. It

climate models (RCMs) or ocean models and both datasets are commonly used as reference datasets for the evaluation of RCMs.Therefore, we extended both datasets with high-resolution river discharge for the period 1979-2018.On the one hand, our discharge data close the water cycle at the land-ocean interface so that the discharges can be used as lateral freshwater input for ocean models applied in the European region.On the other hand, the data can be used to identify trends in discharge that are induced by recent climate change as ERA5 and EOBS22 are rather independent datasets.The experimental setup to generate the discharges was chosen in a way that it could be easily adapted in a climate or Earth system modelling framework.Consequently, the recently developed 5 Min.horizontal resolution version of the hydrological discharge (HD) model was used to simulate discharge.It has already been applied in multiple climate modelling studies and is coupled within several global and regional Earth system models.As the HD model currently does not regard direct human impacts of the river runoff, it is well suited to investigate climate change-related discharge trends.In order to calculate the necessary gridded input fields for the HD model from ERA5 and EOBS22 data, we

Introduction
The 0.5 °resolution of many global observational or quasiobservational datasets is not sufficient for the evaluation of current state-of-the-art regional climate model (RCM) simulations and the forcing of ocean models over Europe.Nowadays, both kinds of models usually operate at higher resolutions.Thus, the ERA5 reanalysis ( Hersbach et al., 2020 ) of the European Centre for Medium-range Weather Forecasts (ECMWF) and E-OBS data ( Cornes et al., 2018 ) are frequently used as reference datasets when RCM results are evaluated on resolutions higher than 0.5 °.In addition, ERA5 data are also commonly used to force regional ocean models ( Bonaduce et al., 2020 ;Muis et al., 2020 ;Wilson et al., 2019 ;Zampieri et al., 2021 ).As ERA5 data do not comprise river discharges, the lateral forcing of freshwater inflow from land is taken from other data sources, such as station data, runoff climatologies, etc.These datasets are not necessarily consistent with the ERA5 forcing over the ocean surface.Moreover, if such data are derived from station data, they are only available for specific rivers and not for all coastal areas.In addition, they might not be representative for the river mouth if the respective station location is too far upstream, which is often the case.In order to mitigate these shortcomings, we decided to extend ERA5 and E-OBS vs. 22 (EOBS22) with simulated high-resolution river discharge.This also allows a consistent assessment of hydrological changes over Europe from these two datasets.
In several recent studies, time series of river discharge were derived from observational data or re-analyses at high-resolution, either globally or specifically over Europe.Tsujino et al. (2018) used river discharge at 0.25 °resolution that was derived from the Japanese 55-year reanalysis ( Kobayashi et al., 2015 ) using the CaMa-Flood model ( Yamazaki et al., 2011( Yamazaki et al., , 2013 ) ).Here, a correction to the runoff was applied ( Suzuki et al., 2017 ) so that the river discharge into the ocean agrees with that reported by Dai et al. (2009) .For the river discharge reanalysis of the Global Flood Awareness System, the LISFLOOD hydrological and channel routing model ( Van Der Knijff et al., 2010 ) was forced with ERA5 data to simulate daily river discharges at 0.1 °global resolution from 1979 until near-real time ( Harrigan et al., 2020 ).Here, LISFLOOD model parameter were calibrated against daily river discharge observations for 1287 catchments globally.The HYdrological Predictions for the Environment (HYPE) distributed hydrological model ( Lindström et al., 2010 ) was forced with daily temperature and precipitation of the ERA-Interim reanalysis ( Dee et al., 2011 ) at 0.75 °resolution to yield daily discharges for 35447 sub-basins with a median size of 214 km ² over Europe ( Donnelly et al., 2016 ).For this E-Hype simulation, precipitation data were adjusted to match the monthly climatological precipitation means from the Global Precipitation Climatology Centre (GPCC) database at 0.5 ° ( Becker et al., 2013 ).E-Hype model parameter were calibrated using observed discharges at 181 river gauges.Recently, the mHM hydrology model ( Samaniego et al., 2010 ) was forced with E-OBS precipitation and temperature data over Europe at 0.625 °resolution from 1950-2019 ( Rakovec and Kumar, 2022 ).Here, model parameters were calibrated with observed discharges focusing on rivers in the catchments of Baltic Sea and North Sea.
These studies commonly have used a full hydrology model (e.g.E-Hype, mHM) or a river routing model (e.g.CaMa-Flood, LISFLOOD) that were calibrated with observed discharges.This approach is useful if the simulated discharges should be as accurate as possible on a large scale.
Here, the calibration may also partially compensate for biases included in the meteorological forcing, e.g. in precipitation.However, if the same setup shall be used in the context of Earth system modelling or climate change studies, it has some disadvantages.With using a full hydrology model, no feedbacks to the atmosphere are considered and the simulated hydrology may be inconsistent with climate model forcing.When using a calibration with observed discharges, the climate change impact on river runoff may be difficult to quantify, as these observations also comprise direct human impacts on the discharge.In addition, the calibration may obscure deficiencies in process representations that may lead to erroneous behavior in ungauged catchments or different climates under global warming conditions ( Hagemann et al., 2020 ).This can be especially the case for rivers where the streamflow is significantly affected by human activities, and where these activities are not sufficiently well represented in the hydrology model.It can also be noted that these high-resolution studies provide hydrology model applications to a single forcing dataset, but no inter-comparison between different forcing datasets.
For our study, we set up two main objectives.The first is to provide consistent high-resolution discharges for the two commonly used datasets ERA5 and E-OBS over Europe.This should be done with a setup that is easily adaptable for climate change studies or in a coupled system modelling approach where the water cycle is closed at the land-ocean interface.Second, we want to analyze the impact of the recent climate change on river runoff.By using two independent forcing datasets based on re-analysis and station data, consistent trend patterns should be a robust measure to indicate this impact.Climate change-related trends are difficult to quantify from observed river discharges at a large scale since most large rivers, especially those for which a long-term streamflow record exists, have been impacted by human influences such as dam construction or land use ( Hartmann et al., 2013 ).Note that Stahl et al. (2010Stahl et al. ( , 2012) ) investigated streamflow trends based on a data set of near-natural streamflow records from more than 400 small catchments in 15 countries across Europe for 1962-2004.However, these studies did not yield any qualitative trends for large European rivers that are based on climate change alone.
In order to extend ERA5 and EOBS22 with high-resolution river discharge, we chose the recently developed version of the Hydrological discharge (HD) model ( Hagemann et al., 2020 ) that is operated at 5 arc minutes horizontal resolution.For the development of this HD model version, no river-specific parameter adjustments were conducted so that the HD model is generally applicable for climate change studies and over ungauged catchments.The HD model has been applied in multiple climate modelling studies and has been coupled to various regional and global Earth system models (ESMs; see Hagemann et al., 2020 ).It is also equipped with an interface for the coupling to an ocean model as the river runoff at the coast can be provided at coastal ocean points on the respective ocean model grid if desired.
The HD model requires gridded fields of surface and subsurface runoff as input with a daily temporal resolution or higher.These variables are not part of EOBS22, as no large-scale observations of these variables exist.ERA5 comprises archived fields of surface and subsurface runoff, but it turned out that the way how the total runoff was distributed between these two runoff fields is not suitable to generate adequate river discharges with the HD model (see Section 2.1 ).Therefore, they need to be calculated by a land surface scheme or hydrology model.Here, we used the HydroPy global hydrology model ( Stacke and Hagemann, 2021 ).Considering discharge simulated from both forcing datasets, we first performed an evaluation with observed discharges at the station locations.Then, we analyzed consistent trend patterns in several discharge statistics at the river mouths.An overview of the used data, models and metrics is provided in Section 2 .In Section 3 , the evaluation of the forcing data and simulated discharges, as well as the trend analysis, are presented and the found trends pattern are discussed.Finally, Section 5 ends with a summary and conclusions.

Data and methods
First, we describe the two atmospheric datasets, ERA5 and E-OBS.This is followed by short introductions to the HD and HydroPy models and the experimental setup used.Then, we refer to the observational data that are used in the evaluation of the model results.Finally, we provide information on the evaluation metrics and trend measures that are analysed.

Atmospheric forcing
2.1.1.ERA5 ERA5 is the fifth generation of atmospheric reanalysis ( Hersbach et al., 2020 ) produced by the European Centre for Medium-Range Weather Forecasts (ECMWF).It provides hourly data on many atmospheric, land-surface, and sea-state parameters at about 31 km resolution.The data are archived in the ECMWF data archive (MARS) and a pertinent subset of the data, interpolated to a regular 0.25 °grid, is available at the C3S Climate Data Store (CDS).ERA5 is produced using 4D-Var data assimilation and model forecasts in the development cycle 41r2 (CY41R2) of ECMWF's Integrated Forecast System (IFS), with 137 hybrid sigma/pressure (model) levels in the vertical and the top level at 0.01 hPa.In the present study, daily surface (or single level) data of precipitation, 2m temperature, downwelling longwave and shortwave radiation, 2m specific humidity, surface pressure and 10m wind are used at the original ERA5 resolution for the period 1979-2018.
Note that ERA5 also comprises archived fields of surface and subsurface runoff.However, ERA5 surface runoff seems to comprise only surface runoff at the top of the soil surface, while all other runoff parts belong to subsurface runoff so that the majority of runoff is allocated as subsurface runoff.This is a hydrological approach valid at the point or field scale, but it is not suitable if larger gridded areas are considered.Here, surface runoff can usually be merged with the fast flow component in the upper layers of soil, such as it is commonly done in climate modelling, where this component is also designated as surface runoff.Therefore, the ERA5 partitioning of total runoff into the two archived fluxes cannot be used to generate adequate river discharges with the HD model.

E-OBS
The E-OBS dataset ( Cornes et al., 2018 ) comprises several daily gridded surface variables at 0.1 °and 0.25 °resolution over Europe covering the area 25 °N-71.5 °N, 25 °W-45 °E.The dataset has been derived from station data collated by the ECA&D (European Climate Assessment & Dataset) initiative ( Klein Tank et al., 2002 ;Klok and Klein Tank, 2009 ).E-OBS is an ensemble dataset that is constructed through a conditional simulation procedure.For each of the members of the ensemble, a spatially correlated random field is produced using a pre-calculated spatial correlation function.The mean across the members is calculated and is provided as the "best-guess" field.For more details on the method, see Cornes et al. (2018) .In the present study, we use the best-guess fields of precipitation and 2m temperature of v. 22 (EOBS22) at 0.1 °resolution for the years 1950-2018.During the evaluation of our results, we noted a problem in the EOBS22 data over Turkey and southern Greece ( Figure 1 ).While the coverage of daily precipitation data is almost 100% over the European land area (except over a small area in northern Turkey that overlaps with the catchments of Sakarya and Kizilirmak), the coverage deteriorates in the year 2005 and gets worse in the year 2014.Similar behaviour is noted for 2m temperature where the coverage is 100% until the year 2005.This deteriorates in the year 2006 and gets worse in 2014, too.Consequently, the related missing data negatively influence the results over the affected catchments in these areas.

The HD model
The Hydrological Discharge (HD) model calculates the lateral transport of water over the land surface to simulate discharge into the oceans.It has been validated and applied in many studies since the publication of its original global 0.5 °version ( Hagemann and Dümenil, 1998 ;Hagemann and Dümenil Gates, 2001 ).Recently, a new version has been developed that can be applied globally at a 5 Min.( ∼8-9 km) resolution.This HD version was applied and validated over Europe by Hagemann et al. (2020) .This study will be referred to as H2020 in the following.The HD model has been coupled to several global and regional ESMs (cf.H2020).It separates the lateral water flow into the three flow processes of overland flow, baseflow, and riverflow ( Figure 2 ).Overland flow and baseflow represent the fast and slow lateral flow processes within a grid box, while riverflow represents the lateral flow between grid boxes.The HD model requires gridded fields of surface and subsurface runoff as input for overland flow and baseflow, respectively, with a temporal resolution of one day or higher.To generate these fields from ERA5 and EOBS22 data, we used the HydroPy global hydrology model (see Section 2.3 ).
Since the time step of the runoff data simulated by the HydroPy model is one day, the HD model time step was also set to one day.However, as the minimum travel time through a grid box is limited by the time step chosen, an internal time step of 0.5 hours is used for river flow.As in H2020, the HD model was set up over the European domain covering the land areas between -11 °W to 69 °E and 27 °N to 72 °N.The domain is shown in Figure 3 together with a few rivers that are specifically noted in this study.These rivers were chosen as they are mentioned in the text for various  reasons.For example, some noticeable behavior was seen over the respective catchments in the analysis of results.
The HD model parameters for overland flow, base flow and river flow are generated as described in Hagemann and Dümenil (1998) and H2020.However, different to the HD model vs. 4 described in H2020, the present vs. 5 ( Hagemann and Ho-Hagemann, 2021 ) utilizes inland water fractions from the ESA CCI Water Bodies Map v4.0 ( Lamarche et al., 2017 ) and wetland fractions from the Global Lakes and Wetlands Database ( Lehner and Döll, 2004 ) instead of the previously used lake and wetlands fractions.A further update in vs. 5 concerns the flow directions and model orography, which were manually corrected for a number of rivers (most of them north of 60 °N).Corrections were based on available GIS data, such as from DIVA ( https://www.diva-gis.org/gdata), CCM River and Catchment Database ( Vogt et al., 2007 ), SMHI (Swedish Meteorological and Hydrological Institute), NVE (Norges vassdragsog energidirektorats), and SYKE (Finnish Environment Institute).
Note that no river-specific parameter adjustments were conducted to generate the HD model parameter so that the HD model is generally applicable for climate change studies in other regions and over ungauged catchments (H2020).As direct human impacts on the discharge are not taken into account, the simulated discharge trends can be considered as caused by the meteorological forcing alone.

HydroPy
HydroPy ( Stacke and Hagemann, 2021 ) is the successor of the Max-Planck Institute for Meteorology's Hydrology Model (MPI-HM; Stacke and Hagemann, 2012 ).It is a state-ofthe-art global hydrology model and its predecessor MPI-HM has contributed to the WATCH Water Model Intercomparison Project (WaterMIP; Haddeland et al., 2011 ) and the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP; Warszawski et al., 2014 ).HydroPy simulates the large-scale water balance in a similar complexity as a modern land surface model, but does not feature a physical energy balance.Instead, processes like snowmelt are represented using empirical relations.The separation of rainfall and snowmelt into surface runoff and infiltration and the generation of drainage (subsurface runoff) is calculated following the improved Arno scheme ( Hagemann and Dümenil Gates, 2003 ).Note that the Arno scheme ( Dümenil and Todini, 1992 ) is widely used in climate and hydrological research (see, e.g., Hagemann and Dümenil Gates, 2003 ).Its infrastructure is improved compared to the MPI-HM, but the hydrological process representation is still very simi-Figure 4 Flowchart showing the main steps of generating simulated discharges using HydroPy and the HD model as well as conducting evaluation and trend analysis.lar.HydroPy is designed to be used with daily data only, as some of its formulations (e.g., the degree-day equation for snowmelt) are not suited for other time steps.We recognize that diurnal variations are important for land-atmosphere interactions.However, using daily data is quite common for large-scale hydrology models ( Telteu et al., 2021 ) and so far, no major simulation deficiencies have been reported for applications on these spatial and temporal scales.Further information about the model performance, evaluation and a comparison between HydroPy and MPI-HM can be found in Stacke and Hagemann (2021) .We choose to include Hy-droPy in our modelling framework as its representation of runoff fluxes is much more suitable to force the HD Model than surface and subsurface runoff archived in ERA5 (see Section 2.1 ).Table 1 provides an overview on the two forcing datasets and their characteristics.Consequently, two experiments were defined that differ in the atmospheric forcing and in how HydroPy was used to generate the HD forcing data of surface and sub-surface runoff:

HD5-ERA5
HydroPy was driven by daily ERA5 forcing data from 1979-2018 to generate daily input fields of surface and subsurface runoff at the ERA5 resolution.It uses precipitation and 2m temperature directly from the ERA5 dataset.Furthermore, potential evapotranspiration (PET) was calculated from ERA5 data in a pre-processing step and used as an additional forcing for HydroPy.Here, we applied the Penman-Monteith equation to calculate a reference evapotranspiration following Allen et al. (1998) .All necessary variables (2m temperature, downwelling shortwave and longwave radiation, 2m specific humidity, surface pressure and 10m wind) were taken from ERA5.This calculation was improved by replacing the constant value for albedo with a distributed field from the LSP2 dataset ( Hagemann, 2002 ).In order to initialize the storages in the HydroPy model and to avoid any drift during the actual simulation period, we conducted a 50-years spin-up simulation by repeatedly using year 1979 of the ERA5 dataset as forcing.Note that we investigated whether the choice of this specific year has a significant impact on the simulation and found no strong effects -neither due to the choice of this year nor due to the choice of our spin-up method in general.

HD5-EOBS
HydroPy was driven by daily EOBS22 data of temperature and precipitation at 0.1 °resolution from 1950-2019.Unfortunately, EOBS22 does not provide all the necessary data fields to calculate PET using one of the more complex conceptual relations like Penman-Monteith ( Allen et al., 1998 ) or Priestley-Taylor ( Priestley and Taylor, 1972 ).Thus, we fall back to a more simple method and compute PET following the approach proposed by ( Thornthwaite, 1948 ).This approach is based on the empirical relation of PET to temperature, represented as a heat index that is calculated using monthly mean temperatures.Furthermore, the average day length at a given location is included.Obviously, this method is less accurate as no land surface characteristics or radiative parameters are taken into account.The period 1950-1978 was used to spin up the model and we only considered the years 1979-2018 for the analysis in this study.

Observed discharge data
Given the resolution of the HD5 model, an adequate simulation of discharge is not expected for small catchments.Therefore, we only considered rivers with catchment areas around 1500 km 2 or larger.Most of the daily discharge data used in this study were provided by the GRDC (Global Runoff Data Centre, 56068 Koblenz, Germany).Further data were obtained from various sources such as listed in Table 2: For each river, the most downstream river gauge was considered for which daily discharge data were available, i.e. from that measurement station, which is closest to the river mouth.

Evaluation metrics
The discharge evaluation of the HD5 simulations was conducted for those grid boxes that correspond to the sta- tion locations within the river network of the respective discharge observations.All evaluation metrics were calculated using simulated and observed time series of daily discharge for the period 1979-2018, thereby covering only those days where observed data are available.We used the mean bias and Kling Gupta efficiency (KGE; Gupta et al., 2009 ;Kling et al., 2012 ) such as described in detail in H2020: The computation of KGE comprises three main components: i) The Pearson correlation coefficient r cor , with an ideal value of one.ii) The variability ratio that is computed by using the coefficients of variation of simulated ( CV sim ) and observed ( CV obs ) data, with an ideal value of one.iii) The ratio between the means of the simulated values μ sim and the observed values μ obs , with an ideal value of one.Note that c 3 corresponds to the mean bias of the simulated values.
KGE ranges from negative infinity to one.Essentially, the closer to one, the more accurate the model is.It also can be noted that a KGE > -0.41 means that the respective simulated discharge represents the observations better than using the observed long-term mean flow ( Knoben et al., 2019 ).The latter is a method that is still used in many ocean model applications.

Trend statistics
For each river, we calculated various trend statistics for the period 1979-2018 at the respective river mouth.The trends were calculated for the variables listed in Table 3 using a linear regression formulation.Note that the number of low flow days is defined as the annual number of days with flows below the environmental water requirement.For every catchment, the environmental water requirement was defined as 30% of its long-term annual mean discharge ( Hagemann et al., 2013 ), thereby utilizing the results of Smakhtin et al. (2004) .

Results
In the following, various metrics and trends were calculated at station locations and river mouths, respectively.However, for the graphical representation, these measures were allocated to the respective catchment areas.

Evaluation of forcing datasets
For the evaluation of precipitation, we used the WFDEI (WATCH forcing data based on ERA-Interim) data ( Weedon et al., 2014 ) as reference dataset.Its monthly means correspond to precipitation data from the GPCC ( Becker et al., 2013 ) that are corrected for gauge undercatch ( Weedon et al., 2011( Weedon et al., , 2014 ) ). Figure 5 compares the annual mean precipitation of ERA5 and EOBS22 to WFDEI data for 1979-2016, the time period for which WFDEI data are available.To allow for an easier comparison, the precipitation was interpolated to the HD5 model grid using conservative remapping ( Figure 5 , left column).In addition, the precipitation bias of ERA5 and EOBS22 is also projected on the respective catchment areas ( Figure 5 , right column).
While the general pattern agrees between the three datasets, EOBS22 exhibits a widespread dry bias over almost the entire European EOBS domain.This dry bias is likely related to the fact that EOBS22 data are based on gauge measurements from the ECA&D station network, but these data have not been corrected for gauge undercatch.The underestimated precipitation also translates in a dry runoff bias for most of the catchments ( Figure 6 , upper right).Exceptions are catchments where large anthropogenic water abstractions occur such as the Ebro ( Merchán et al., 2013 ), Don ( Khublaryan, 2009 ) and many Turkish rivers ( Akbulut et al., 2009 ).As these anthropogenic water withdrawals are not regarded by the HD model, the simulated discharge is considerably larger than the observations.The missing data problem in EOBS22 (cf.Section 2.1 ) also affects the discharge simulations over Turkey.While there is a general overestimation of discharge for Turkish rivers until 2004, this turns into an underestimation of discharge for the second period 2005-2018 ( Figure 6 ).This will lead to erroneous decreasing discharge trends.Here, it must be noted that most observed Turkish discharges are only available until 2011 (a few rivers have data until 2015), so that the discharge bias for the whole period 1979-2018 is strongly determined by the bias for the period before the missing data problem occurs.
ERA5 tends to overestimate precipitation over Scandinavia, the northwestern British Isles, eastern Turkey and larger parts of Eastern Europe compared to WFDEI.This also largely contributes to the wet biases over Eastern Europe and eastern Turkey.However, the simulated discharge of HD5-ERA5 tends to be underestimated for many rivers in Scandinavia and the northwestern British Isles, which points to a general overestimation of evapotranspiration by Hy-droPy using the ERA5 forcing over these areas.
Except for rivers with large abstractions of water, biases in total runoff and discharge agree in their long-term averages.These biases are generated by biases in the forcing data and may be introduced by uncertainties in the land surface representation of the HydroPy model.In Hy-droPy, sources of uncertainty range from spatial resolution and sub-grid representations to process simplifications, e.g.there is no discrimination between different types of surface water bodies.Moreover, for HD5-EOBS, a major limitation is imposed for the calculation of potential evapotranspiration, as only temperature is available.Here, a simpler method than for HD5-ERA had to be used to estimate potential evapotranspiration (see Section 2.4 ).If more variables (e.g.radiation, wind, humidity) would be available, evapotranspiration can be represented much better by HydroPy as shown in Stacke and Hagemann (2021) .
Note that runoff biases impose an upper limit on the accuracy of simulated discharge.Therefore, if the KGE of a specific river with a runoff bias b R is considered, the maximum KGE that can be yielded after simulating the discharge is 1-b R .

Evaluation of simulated discharge
For both simulations, the general behavior of discharge is captured well (KGE > 0.4) for many European rivers, especially in northern Iberia, Western and Central Europe ( Figure 7 , left panel).HD5-ERA5 shows also a good performance over northern Russia.
As expected (cf.H2020), larger deviations of the simulated from the observed discharges occur for rivers that are strongly affected by human impacts such as water abstractions, e.g. for irrigation, and regulation, e.g. by dams.This is the case for many Scandinavian and Turkish rivers as well as the Volga and Don rivers.For most parts of Northern Europe, HD5-EOBS is worse than HD5-ERA5.This is related to the large dry bias of HD5-EOBS ( Figure 6 ), which is likely caused by the underestimated precipitation over this area (see Section 3.1 ).Because evapotranspiration is not moisture limited over Northern Europe ( Teuling et al., 2009 ), absolute biases in precipitation directly translate into runoff biases.In addition, the simpler calculation of PET in HD5-EOBS leads to increased evapotranspiration over Northern and Western Europe compared to HD5-ERA5 (Figure S1).As the HD5-ERA5 evapotranspiration already seems to be

Trends in simulated discharge
Figure 8 shows simulated trends in the annual minimum ( Q min ), mean ( Q mean ) and maximum ( Q max ) discharge.In HD5-ERA5, a general drying is indicated for these statistics over Central, Eastern and Southern Europe.Over the latter area, exceptions can be noted for southeastern Spain, southern Italy and the Balkan (except southern Greece) where Q max and Q mean tend to increase.Here, increases of Q min occur only over southeastern Spain and the Tejo River, the Italian Alpine rivers and a few Balkan rivers.These general trend patterns are also the case in HD5-EOBS, except for the rivers Vistula and Danube that show wetting trends.In addition, the wetting trends over southeastern Spain, southern Italy and the Balkan are more pronounced than in HD5-ERA5 and are also visible in the Q min trends (except for southeastern Spain).
For Northern Europe, trends are more diverse.In both experiments, only Q max tends to get lower, while Q min mainly becomes wetter, especially in HD5-ERA5.In Q mean , no widespread trend signal is seen.Northern UK generally  tends to get wetter in each of the three statistics, which also applies to Ireland for HD5-ERA5.For HD5-EOBS, only Q max increases for Ireland, while Q min and Q mean are decreasing.
Note that the HD5-EOBS trend patterns over Turkey are not robust due to the impact of the missing EOBS22 data after 2004 (cf.Section 2.1 ).The trend patterns for the period 1979-2004 are very similar in HD5-ERA5 and HD5-EOBS (not shown) over Turkey for all statistics.However, the trend statistics for the shorter earlier period are affected by decadal variability, and they partially differ from those for 1979-2018.Therefore, we did not merge trends from the two periods to overcome the missing data problem.
For most areas, the trend pattern in the winter mean discharge ( Figure 9 ) follows the respective pattern of Q min .Only over the Iberian Peninsula, an overall drying is indicated that is more pronounced (HD5-EOBS, lower right) or even opposite (HD5-ERA5, lower left) to the respective trends in Q min .The summer mean trend patterns ( Figure 9 , upper panels) are rather similar to those in Q mean and Q max for both experiments.The most notable difference is a more pronounced drying in the summer mean over Central Europe in HD5-ERA5.In HD5-ERA5, the largest seasonal drying trends occur in summer and autumn, while in HD5-EOBS, these tend to occur in spring (British Isles and Central Europe) and autumn (not shown).
In both experiments, an increase in the number of low flow days ( Figure 10 , 1 st row) and longest low flow period (not shown) is seen for many Southern and Central European rivers, and a general decrease over Northern Europe.Considering the time of the year when 50% of the annual flow volume is exceeded ( Figure 10 , 2 nd row), the event tends to occur earlier for Central and Southern Europe, whereas the trend is more pronounced in HD5-ERA5.Over Northern Europe, no widespread trend pattern is seen.On the one hand, it can be noted that the annual peak flows are reduced over Northern Europe and western Iberia if compared to the respective annual mean discharges ( Figure 10 , 3 rd row).On the other hand, the peak-to-mean relation tends to increase over the British Isles, Western Europe and eastern Iberia.Other parts of Europe do not show a consistent pattern in both experiments.

Comparison of simulated discharges with observed climatology
In order to get an impression on the applicability of the simulated discharges as lateral freshwater forcing for ocean models, we compare them to the observed discharge climatology as this has been commonly used as an ocean model forcing.Note that nowadays such a forcing became less common, especially for models of European marginal seas such as the Baltic Sea (M.Meier, pers. comm., 2022).Therefore, we calculated the KGE for the case that a monthly climatology of observed discharge is used (experiment CLIM; Figure 11 ).For each station, the monthly climatology was calculated from available discharge observations during 1979-2018.Then, this climatology was simply converted into daily data for the same period, i.e. the discharge on each day within a month was set to the respective climatological value of the month.Comparing the CLIM KGE ( Figure 11 ) to the two experiments shows that both experiments yield a better KGE ( Figure 7 ) for most rivers in Western and Central Europe.In Southeastern and Northern Europe, the CLIM KGE is similar to the HD5-ERA5 KGE, while the CLIM KGE is noticeably high for Northeastern European and Russian rivers, so that it is about equal to or better than the KGE of the simulated discharges.
It must be pointed out that this comparison is somewhat unfair as the CLIM bias is zero by default.Consequently, KGE is determined only by the correlation and the variability ratio.The latter is largely underestimated in CLIM due to the missing daily and interannual variability.Consequently, high KGE can be yielded if the discharge curve is rather smooth and has a strong seasonally reoccurring behavior, such as for the Northeastern European rivers that are characterized by large snowmelt-induced peaks in the spring.In addition, using CLIM as freshwater forcing at the river mouths would introduce errors in the freshwater forcing at rivers where the gauge station is not close to the river mouth.For those rivers, the relative quality of the simulated discharge may be comparatively higher as the HD model simulates the discharge also at the river mouth.The relative quality of CLIM compared to the two experiments may also be lower in years where the actual discharge deviates from CLIM, either due to interannual variability or due to missing data in the available discharge observations.

Trends in simulated discharge
We found that for Northern Europe, discharges are increased in winter as well as for Q min .This is consistent as for most areas in this region, the lowest discharge occurs in the winter when frozen soils and solid precipitation prevail.The warming in the recent decades was accompanied by increased winter precipitation ( European Environment Agency, 2017 ) but lower snowfall ratios ( Hartmann et al., 2013 ).Consequently, on the one hand, more discharge is generated during winter, thereby increasing Q min and reducing the number of low flow days.On the other hand, the lower snowfall ratio led to a reduction in snowmelt amounts so that the related peaks in discharge (i.e.Q max ) are reduced.As there are no widespread trend signals in Q mean , these reduced peaks also became lower in compari-son to the mean discharge volumes, which is expressed by the peak-to-mean relation.
Over Central and Southern Europe, we found a general drying in the mean, maximum and minimum discharges that is most pronounced during summer and autumn.The recent warming caused a general lengthening of the growing season ( European Environment Agency, 2017 ) that is also accompanied by an earlier onset of transpiration in the spring.Together with the tendencies of reduced summer and an- nual (only Southern Europe) precipitation ( European Environment Agency, 2017 ), this leads to drier soils and less runoff.Consequently, the number of low flow days has increased.As the drying is most pronounced in summer and autumn, i.e. in the second half of the year, the 50% threshold of the annual flow volume is reached earlier within a year.
However, we noted a few areas in Southern Europe where Q min tends to increase.For the Balkan, this seems to be related to a slight increase in annual and summer precipitation.For southeastern Spain, it can be noted that the recent decrease in annual precipitation over the Iberian Peninsula is rather low and insignificant, i.e. much lower than in western Spain ( European Environment Agency, 2017 ).Moreover, Peña-Angulo et al. (2020) found no consistent long-term precipitation trends in Southwestern Europe.Therefore, they argued that the observed decrease in precipitation during 1961-2000 in the ECA&D data is likely affected by uncertainties introduced in this dataset due to the low density of stations, their varying densities over time, and the lack of temporal homogeneity of some series.Further, the results of Llasat et al. (2021) indicated a predominant increase in convective precipitation over eastern Spain.They also stated that the contribution of convective precipitation to total precipitation in the East of Spain exceeds 10% annually and can reach 100% in late summer.As the summer is the dry season in this area, the runoff created by convective precipitation events is the major source of discharge during this season, so that an increase in those events may explain the increase in Q min .On the other hand, these trends may be caused by a large interannual or decadal variability.To elucidate this, we considered the annual time series of Q min for the rivers Tejo and Guadalquivir ( Figure 12 ).For Tejo, Q min varies largely.HD5-ERA5 misses the wet peak in 1979 and tends to have wetter low flows after 2000, so that it shows a slightly increasing trend while HD5-EOBS shows a negative trend.For Guadalquivir, the increasing trend in HD5-ERA5 seems to be mainly caused by the increasing wetness of the three events around the years 1990, 1997 and 2010, when Q min was noticeable wetter than usual.However, HD5-EOBS includes the same events, but also shows a similar event in 1979 so that no positive trend is shown in Figure 12 .According to discharge observations, this wet event in 1979 actually occurred for both rivers so that the Q min trends of HD5-ERA5 seem erroneous.Note that none of these trends is significant.We will take up the issue of significance in Section 5 .
Note that the found trend patterns are largely consistent with the results of ( Stahl et al., 2010 ) on near-natural European streamflow trends for the period 1962-2004 (cf.Section 1 ).As those data are rather sparse, they cover mostly rather small catchments and no rivers were utilized in Ireland, Poland, Southern Spain, Italy, Hungary and Southeastern Europe.Stahl et al. ( 2010) yielded a regional coherent pattern of annual streamflow trends with negative trends in southern and eastern regions, and generally positive trends elsewhere.Our trends in Q mean and the winter mean discharge are generally similar but increases in Central Europe were not obtained in HD5-ERA5 and HD5-EOBS.However, these might be obscured as we considered only large catchments whose discharge represents an integral over a large area.Our Q min trends match rather well with the trends for mean monthly streamflow for the month of the regime minimum in Stahl et al. (2010) , especially the increasing trends seen in the Alpine region in HD5-ERA5.Note that the comparison to the results of Stahl et al. (2010) is hampered by the use of different periods and different catchment sizes (see above).
As mentioned above, the near-natural rivers considered by Stahl et al. (2010) comprise only rather small catchments and do not cover large parts of Europe.Hence, Stahl et al. (2012Stahl et al. ( ) calculated trends (1963Stahl et al. ( -2000) ) based on an ensemble of eight global hydrology models (GHMs; 0.5 °resolution, only one was calibrated) driven by WATCH forcing data ( Weedon et al., 2011 ).While for Q min , the trend patterns of HD5-ERA5 and HD5-EOBS generally agree with the GHM ensemble mean, there are noticeable differences for Q mean over Central Europe where the ensemble mean shows increasing trends while HD5-ERA5 and HD5-EOBS show no or decreasing trends.For Q max , similar differences also occur over parts of Northern Europe.However, from Figure 4 of Stahl et al. (2012) , it seems that the Q mean and Q max trends from our simulations agree better with the observed locations than with those from the ensemble mean.
We also compared our results with trends obtained from some of the model-based products mentioned in Section 1 .The mHM hydrology model was forced by EOBS vs. 21e data for the same period as used for HD5-EOBS ( Rakovec and Kumar, 2022 ).Consequently, the trends in this mHM-EOBS discharge (Figure S2 -left column) are rather similar to HD5-EOBS.However, the mHM-EOBS Q mean and Q max have more increasing trends over Northern and Western Scandinavia and the decreases over Southwestern Europe are more pronounced.The same applies for the mHM-EOBS Q min , with also the increases over Scandinavia being more pronounced, which is actually more similar to HD5-ERA5 than to HD5-EOBS over this region.Hence, other trend statistics look also rather similar.
For E-Hype ( Donnelly et al., 2016 ), we had access to discharge data from 1981-2010.E-Hype Q mean and Q max show similar patterns (Figure S2 -right column) as in our results only over Western Europe, but more pronounced decreasing trends in southwestern Scandinavia, and increasing trends over many catchments in Eastern and Southeastern Europe.In addition, Q mean shows also larger increases over the Iberian Peninsula.E-Hype Q min shows rather different patterns (Figure S2 -lower right panel), having increases over most eastern European rivers and the Iberian Peninsula, but decreases over Central Sweden.Hence, other trend statistics look also rather different.These differences may partially be related to the rather coarse resolution (0.75 °) of the forcing.A more detailed investigation of the reasons for these differences is beyond the scope of the present study.

Summary and conclusions
In the present study, we generated high-resolution river discharge data over Europe from two rather independent meteorological datasets, ERA5 and EOBS22, for the period 1979-2018.To obtain the results, ERA5 and EOBS22 were used to force the HydroPy hydrology model and the HD model.This experimental setup was chosen, as it can also be used to close the hydrological cycle at the land-ocean interface within an ESM.In this way, we generated river discharges that are consistent with the meteorological forcing, which is especially relevant for ERA5, which is frequently used as an atmospheric forcing for ocean models, see, e.g., ( Bonaduce et al., 2020 ).After a short evaluation, the simulated discharges were then used to identify consistent trend patterns in both experiments, HD5-ERA5 and HD5-EOBS.
Overall, the general behavior of discharge is captured well with a KGE larger than 0.4 in both experiments for many European rivers.These results are consistent with earlier results of Hagemann et al. ( 2020) using different meteorological forcing data.Based on the KGE, the simulated discharge is equal to or better than a freshwater forcing derived from an observed monthly climatology (CLIM) for many rivers, especially in the case of HD5-ERA5.For EOBS22, two problems were identified.On the one hand, the analysis of trends in HD5-EOBS over Southeastern Europe was hampered by missing data in EOBS22 after 2004.On the other hand, the missing undercatch correction in station data used in EOBS22 led to a widespread low bias in simulated discharge of HD5-EOBS.In addition, noticeable wet biases for some catchments in both experiments are related to human water abstractions that are currently not regarded in the HD model.
In both experiments, we found increases in winter and annual minimum discharges in Northern Europe that led to a reduced amount of low flow days.As lower snowfall ratios cause a reduction in snowmelt-related discharge peaks in spring, and, thus in Q max , the peak-to-mean relations also became lower.Over Central and Southern Europe, we noted a general drying in the annual discharge statistics ( Q mean , Q max and Q min ) that is seasonally more pronounced in the summer and autumn discharges.This is accompanied by an increase in low flow days and an earlier exceedance of the 50% annual flow volume threshold.
We emphasize that some of our trend statistics for the considered 40-years period can be strongly impacted by interannual or decadal variability.Large interannual variability is equivalent to a high noise level that makes it difficult to separate weaker signals, i.e. trends, from it.Consequently, if only annual time series at specific river mouths are considered, a significance test will yield that some of these trends are not significant.Therefore, our approach in the evaluation of trends was to consider spatial trend patterns, i.e. spatial clusters of rivers with trends in the same direction in two datasets.In this way, having robust spatial patterns is kind of increasing the sample size beyond the 40 years at each river location.In addition, we looked for plausible causes if specific rivers showed noticeable deviations from these spatial patterns.However, we acknowledge that some of these deviations may still occur by chance, especially if the interannual variability is rather large.In order to allow a more complete picture, we provide patterns with significant trends for Q mean , Q max and Q min in Figure S3 as a supplement.
A large decadal variability may lead to the identification of trends that are not caused by climate change but are merely induced by the sequence of 'high' and 'low' decades in the considered trend statistic.Using 40-years long time series should lead to a reduced probability that such trends occur.In addition, we also checked (cf.Section 4.2 ) whether the trends patterns are consistent across the various trend statistics.However, for some rivers, such misleading trends might still occur, especially if the interannual variability is also high, such as for Q min of the Guadalquivir (see Section 4.2 ).
In the future, we suggest a comparison of trends at gauge station locations to observe discharge trends in order to identify and quantify changes that are not related to climate change.Results from such a study may also help to implement processes into the HD model that represent direct human impacts on the river.This includes the implementation of dams and reservoirs (based on available global databases), their management of river flow regulations as well as modules to simulate water withdrawals, e.g. for irrigation ( Hagemann et al., 2020 ).
It may also be interesting to utilize the WFDE5 dataset ( Cucchi et al., 2020 ) as forcing, where the WATCH Forcing Data (WFD) methodology ( Weedon et al., 2011 ) was applied to surface meteorological variables from ERA5 to yield biascorrected time series.Unfortunately, these data are currently only available at 0.5 °resolution.Here, a utilization at a higher resolution would be desirable, such as e.g. on 0.25 °or the original ERA5 resolution (cf.Section 2.1 ).
As mentioned above, the simulated discharges of both experiments can be used to force an ocean model.Thus, the impacts of discharge trends on marine states and processes may be investigated as well as how uncertainties imposed by differences between the two experiments are propagated into the marine system.However, as the HD model discharge is consistent with the runoff forcing from the driving climate, land surface or hydrology model, it also comprises biases imposed by biases in precipitation and/or uncertainties in the land surface representation of the respective model.For semi-enclosed seas like the Baltic Sea, a bias correction of the simulated discharge may be necessary if the basin-wide biases are too large.For Baltic Sea models, for instance, the mean long-term bias of river discharge needs to be smaller than 7% (M.Meier, pers. comm., 2022).Otherwise, the ocean model will drift into a different climate state compared to the observed state.Using HELCOM data on riverine inflow into the Baltic Sea from 1995-2018 as reference (15752 m 3 /s; Svendsen and Gustafsson, 2020 ), HD5-ERA5 falls into this limit with a small bias of about + 2% during 1979-2018.However, due to the widespread dry bias of HD5-EOBS in this region (-43%; Figure 6 ), a bias correction would be necessary if HD5-EOBS is used to force a Baltic Sea ocean model.
Finally, our setup is also useful to generate highresolution river runoff data consistent with the meteorological forcing for historical periods and future scenarios from any climate model data instead of having to rely on observed time series.First, this this can be done by using only the HD model in a coupled setup, e.g., in the GCOAST AHOI system ( Ho-Hagemann et al., 2020 ).Second, the HD model can be utilized in standalone mode with forcing from hindcasts, analyses or climate simulations, such as generated by the suite of models from the German Weather Service (DWD), i.e.Cosmo-CLM ( Rockel et al., 2008 ) or ICON ( Zängl et al., 2014 ).Third, if only meteorological variables are available in the forcing data, the present setup of combining HydroPy and the HD model may be applied, e.g. using the WFDE5 data as mentioned above.

Figure 1
Figure 1 Daily coverage of EOBS22 data for precipitation (upper left panel) and 2m temperature (upper right panel) from 2005-2018 and 2006-2018, respectively.The lower panel shows the annual time series of this coverage for both variables over the Turkey region encompassing the area between 20 °E-45 °E and 35 °N-45 °N.

Figure 2
Figure 2 Flow processes (in red and yellow) within the HD model.Flow processes labelled in other colors are added for completeness, but are not part of the HD model.

Figure 3
Figure 3 European HD model domain and catchment areas for selected rivers.

Figure 4
Figure 4 summarizes the four main steps of generating simulated discharges and conducting their evaluation and analysis: 1) Preparation of HD model forcing: Choose an atmospheric forcing dataset (ERA5 and EOBS22) and use HydroPy to generate the forcing for the HD model.2) Simulation with the HD model: Interpolate the forcing data of surface and sub-surface runoff to the HD model grid and simulate daily discharges with the HD model.3) Evaluation of results: Compare simulated and observed discharges at station locations and calculate various evaluation metrics.4) Calculation of trends: Calculate various trends measures at the river mouth locations.

Figure 10
Figure 10 Trends in the number of low flow days [d/a] (1 st row), in the time of the year when 50% of the annual discharge volume is reached [d/a] (2 nd row) and in the peak to mean relation [%/a] (3 rd row) for HD5-ERA5 (left column) and HD5-EOBS (right column) from 1979-2018.

Figure 11
Figure 11 Kling Gupta efficiencies for using the observed monthly discharge climatology in CLIM during 1979-2018.

Figure 12
Figure 12 Annual time series of Q min at the river mouths of Tejo (upper panel) and Guadalquivir (lower panel).

Table 1
Model experiments and their forcing data.

Table 2
Data sources for discharge observations.

Table 3
Discharge variables considered for the trend calculation.Annual maximum divided by annual mean = Q max / Q mean [%/a] Time of 50% volume Time of the year when 50% of the annual discharge volume is exceeded [d/a]