HERA: a high-resolution pan-European hydrological reanalysis (1951–2020)

Since 1950, anthropogenic activities have altered the climate, land cover, soil properties, channel morphologies, and water management in the river basins of Europe. This has resulted in significant changes in hydrological conditions. The availability of consistent estimates of river flow at the global and continental levels is a necessity for assessing changes in the hydrological cycle. To overcome limitations posed by observations (incomplete records, inhomogeneous spatial coverage), we simulate river discharge for Europe for the period 1951–2020 using a state-of-the-art hydrological modelling approach. We use the new European set-up of the OS LISFLOOD model, running at 1 arcmin (≈1.8 km) with 6-hourly time steps. The hydrological model is forced by climate reanalysis data (ERA5-Land) that are bias-corrected and downscaled to the model resolution with gridded weather observations. The model also incorporates 72 surface field maps representing catchment morphology, vegetation, soil properties, land use, water demand, lakes, and reservoirs. Inputs related to human activities are evolving through time to emulate societal changes. The resulting Hydrological European ReAnalysis (HERA) provides 6-hourly river discharge for 282 521 river pixels with an upstream area >100 km². We assess its skill using 2448 river gauging stations distributed across Europe. Overall, HERA delivers satisfying results (${\text{median KGE}}^{\prime }=\mathrm{0.55}$), despite a general underestimation of observed mean discharges ($\text{mean bias}=-\mathrm{13.1}$ %), and demonstrates a capacity to reproduce statistics of observed extreme flows. The performance of HERA increases through time and with catchment size, and it varies in space depending on reservoir influence and model calibration. The fine spatial and temporal resolution results in an enhanced performance compared to previous hydrological reanalysis based on OS LISFLOOD for small- to medium-scale catchments (100–10 000 km²). HERA is the first publicly available long-term, high-resolution hydrological reanalysis for Europe. Despite its limitations, HERA enables the analysis of hydrological dynamics related to extremes, human influences, and climate change at a continental scale while maintaining local relevance. It also creates the opportunity to study these dynamics in ungauged catchments across Europe. The HERA hydrological reanalysis and its climate and dynamic socio-economic inputs are available via the JRC data catalogue: https://doi.org/10.2905/a605a675-9444-4017-8b34-d66be5b18c95 (Tilloy et al., 2024).

1 Introduction

In the last century, Europe has experienced a growth in its population, economy, and urbanized area (Li et al., 2021; Paprotny and Mengel, 2023). Recent decades have also witnessed a rapid rise in global air temperature that is attributable to anthropogenic activities (IPCC, 2023). These evolving conditions have significantly changed flows in European streams and rivers (Barker et al., 2019; Gudmundsson et al., 2021; Vicente-Serrano et al., 2019; Wang et al., 2024), leading to multiple challenges for hydrological sciences related, for example, to long-term variability, climate change, extremes, or human alterations of the water cycle (Blöschl et al., 2019a). In order to assess the impacts of these changes, hydrologists need consistent, reliable, and long hydrological series. Observations, despite continuous improvements (Blöschl et al., 2019b; Ekolu et al., 2022), can hamper the analysis of pan-European long-term trends due to sparse spatial distributions in some regions and temporal discontinuities. One option for overcoming these limitations is to rely on a suite of models (climate, hydrological, and land use) to simulate past hydrological conditions and interpret changing dynamics in the hydrological cycle in connection with rapidly changing human systems (e.g. Richards and Gutierrez-Arellano, 2022). This article introduces the Hydrological European ReAnalysis (HERA) for the period 1951–2020, providing consistent estimates of river flow for European rivers at high spatial and temporal resolution.

Hydrological models are essential tools for understanding and characterizing processes related to the water cycle (e.g. flood and drought forecasting). In the past 3 decades, there have been efforts to develop models that are able to simulate hydrological processes at a large scale (continental to global scales). A myriad of these global hydrological models (GHMs), differing in their conceptualization, now exist (Beck et al., 2017; Sood and Smakhtin, 2015; Kauffeldt et al., 2016; Prudhomme et al., 2011). The nature of GHMs implies that they are usually run at a coarse spatial resolution (e.g. 0.5°), limiting their relevance for local and regional water resource problems (Sood and Smakhtin, 2015). Nonetheless, the development of GHMs has been fuelled by continuous improvements in remote sensing technologies and processing power (Yang et al., 2021). Remote sensing technologies provide high-resolution input for hydrological models, such as land use and vegetation properties. The advancements in computational capabilities have allowed us to refine the spatial and temporal scale of hydrological models, enabling a more accurate representation of surface and sub-surface processes and reducing modelling uncertainties (Wood et al., 2011). In this context, HERA falls within a global effort towards the development of hyper-resolution (1 km and below) land surface and hydrological models at the continental (Hoch et al., 2023; O'Neill et al., 2021) and global (Hanasaki et al., 2022) scales.

A key hindrance to simulating past river flows has been the availability of meteorological inputs for hydrological models. Of the potential inputs, climate reanalysis offers several advantages: temporal coverage (typically spanning several decades) and a large number of variables (e.g. precipitation, wind speed, or temperature) that are physically consistent with a homogeneous spatio-temporal resolution. The reanalysis data are outputs of climate models calibrated on observed data worldwide (Brönnimann et al., 2018). Here we use ERA5-Land, the land component of ERA5 (Muñoz-Sabater et al., 2021). One main advantage of ERA5-Land compared to ERA5 is its horizontal resolution, which is 9 km globally, compared to 31 km in ERA5. This enhanced resolution is obtained by downscaling meteorological variables from ERA5. The temporal resolution is hourly as in ERA5. Nonetheless, reanalysis data are obtained from short-term model forecasts and can be affected by forecast errors (Pfahl and Wernli, 2012). Variables produced in ERA5 are averages over grid cells. This averaging combined with the relatively coarse resolution of ERA5/ERA5-Land often smooths local extremes (Donat et al., 2014; Tilloy et al., 2022). To tackle this issue, we downscale and bias-correct ERA5-Land with a gridded observational dataset, EMO-1 (Thiemig et al., 2022) (Sect. 2.2).

In the context of the European Flood Awareness System (EFAS), an operational system for European flood monitoring and forecasting (https://www.efas.eu, last access: 2 December 2024), there have been recent efforts to develop more detailed surface fields (e.g. land use or vegetation) (Choulga et al., 2024) and observational climate inputs (Thiemig et al., 2022) at a spatial resolution of 1 arcmin (1^′, 0.0167°, typically 1.5–3 km² over Europe). These developments have come with improvements in the OS LISFLOOD hydrological model underpinning EFAS. OS LISFLOOD is a spatially distributed grid-based hydrological and channel routing model which was initially developed for flood forecasting and flood risk assessment (Burek et al., 2013). However, it is also able to model effects of land use change, climate change, and river regulation measures and has been used in a wide range of hydrological applications, such as mapping populations under water stress in relation to how much water is reserved for the environment (Vanham et al., 2021) and projecting droughts in view of climate change (Cammalleri et al., 2020b). It is also used in the generation of the GLOFAS-ERA5 hydrological reanalysis (Harrigan et al., 2020).

Therefore, this article brings together improvements from diverse fields (i.e. remote sensing, climate modelling, machine learning, and hydrology) to generate a state-of-the-art hydrological reanalysis for a European domain that covers EU27 countries, the UK, Switzerland, Iceland, Norway, and the Balkan countries (Serbia, Montenegro, Bosnia-Herzegovina, Kosovo, North Macedonia, and Albania) over the past 70 years. HERA aims to reproduce as accurately as possible the evolution of the hydrological landscape of Europe by using the latest development of OS LISFLOOD (improvements in processing speed, spatial and temporal resolution, and calibration), which is also used in the generation of the latest EFAS v5.0 reanalysis (1991–2022) (Decremer et al., 2023) (Sect. 2.1). Climate inputs are derived from ERA5-Land and are bias-corrected and downscaled to 1 arcmin to improve the representation of extremes (Sect. 2.2). We generated dynamic socio-economic inputs (water demand, land use, and reservoir maps) to capture the effect of human activities on the water cycle (Sect. 2.3). These developments make this dataset the first publicly available long-term pan-European hydrological reanalysis that takes into account the evolving socio-economic conditions that have altered the hydrological cycle since 1951. In Sect. 3, we assess the performance of HERA against observations from 2448 river gauges in Europe.

2 Method

The modelling framework developed to generate the HERA dataset is presented in a flowchart in Fig. 1. The framework is organized around the OS LISFLOOD hydrological model that is used to simulate river discharge. For this run, we use calibrated parameters for the European setting of OS LISFLOOD developed by ECMWF in the context of the EFAS v5.0 calibration (CEMS-Flood online documentation, 2023). We first introduce OS LISFLOOD and its calibration procedure (Sect. 2.1). Figure 1 also displays the main input of OS LISFLOOD: high-resolution climate inputs (Sect. 2.2) and state-of-the-art static (Sect. 2.3.1) and dynamic socio-economic maps (Sect. 2.3).

Figure 1 Flowchart of the framework employed in the generation of HERA. The numbers relate to the section in which each component of the framework is presented.

2.1 Hydrological modelling

2.1.1 The OS LISFLOOD model

Here, we simulate sub-daily continuous streamflow time series over Europe by means of the OS LISFLOOD model (Burek et al., 2013; Van Der Knijff et al., 2008). This is a spatially distributed, semi-physical rainfall-runoff model combined with a routing module for river channels (Dottori et al., 2022). The model has been developed by the Joint Research Centre (JRC) since the late 1990s and is used operationally for large-scale flood forecasting in EFAS and the Global Flood Awareness System (GLOFAS). OS LISFLOOD has also been used in drought monitoring (Cammalleri et al., 2020a, 2017) to assess the effect of flood adaptation measures, environmental flow protection, or climate change (Mentaschi et al., 2020; Vanham et al., 2022). Since 2019, the model has been open-source and available on GitHub along with a set of auxiliary tools (https://github.com/ec-jrc/lisflood-code, last access: 2 December 2024). OS LISFLOOD is composed of the following main components:

• three soil layers (superficial, upper, lower) for water balance modelling;

• sub-models for the simulation of groundwater and subsurface flow (using two parallel interconnected reservoirs);

• a sub-model for the routing of surface runoff to the nearest river channel; and

• a sub-model for the routing of channel flow.

Other processes, e.g. snowmelt, infiltration, rainfall interception, leaf drainage, evaporation and water uptake by vegetation, surface runoff, and exchange of soil moisture between soil layers, are also simulated by the model (European Commission, Joint Research Centre (JRC), 2025a). OS LISFLOOD is also able to model lakes and reservoirs.

In this work, we use the latest version of OS LISFLOOD (v4.1.2, January 2023), which includes upgrades compared to previous versions in the hydrological routines and improvements in the modelling of water abstraction for anthropogenic use. Moreover, OS LISFLOOD v4.1.2 benefits from improvements in the management of large inputs and computational performance. Figure 2 displays the domain for which data were retained in HERA. This comprises 42 European countries and excludes non-EU countries of the former Soviet Union, countries in northern Africa and the Middle East, and Türkiye, which are included in the EFAS domain. Moreover, HERA uses the same domain as the Historical Analysis of Natural Hazards in Europe (HANZE) database (Paprotny and Mengel, 2023; Paprotny et al., 2023). We run the model using the 1^′ grid in EFAS v5.0 (Decremer et al., 2023). The temporal resolution of the simulation is 6-hourly, which has been the standard for EFAS since 2020. Due to the size and spatial resolution of our domain, combined with the 6-hourly time steps, we divide the simulations into 71-year chunks based on calendar year starting on 3 January 1950. To estimate the initial model state, we performed a 71-year pre-run. In particular, we used the pre-run to initialize the soil and upper groundwater zone storages and to derive the average inflow into the lower zone and discharge, which represent the theoretical steady-state storage. Due to the rapidly evolving socio-economic conditions in catchments of Europe, we change the input socio-economic maps at the start of every new calendar year of the simulation (Sect. 2.4). This differs from the standard EFAS settings, which assume static land use and a reservoir network, and only varies the water demand values. At the start of every calendar year, the model is initialized with the state variables from the last time step of the previous year (warm start). As water volumes at the first time step in the channels are not known, the model sets a conventional initial volume (OS LISFLOOD uses half-bankfull), leading to unrealistic initial discharge in some catchments. We therefore removed the first simulation year (1950) from the final dataset. Further, we only retained simulations for river pixels with an upstream area greater than 100 km², resulting in simulations in the 282 521 river pixels displayed in Fig. 2.

Figure 2 River network (rivers with an upstream area >100 km²) on which discharge data have been generated. The HERA domain (in which data are provided) is confined by the red-bordered area.

2.1.2 Model calibration

In this work, we also take advantage of the new EFAS v5.0 calibration that was completed in December 2022 by ECMWF. The calibration was performed using the EMO-1 meteorological dataset (Thiemig et al., 2022) over the period 1990–2021, with a focus on high flows. The modified Kling–Gupta efficiency (KGE^′; Gupta et al., 2009; Kling et al., 2012) was used as a skill metric. Discharge data at 1903 stations, identified through a selection process based on several criteria (CEMS-Flood online documentation, 2023), were used to calibrate the OS LISFLOOD model over Europe. Sub-daily data are always preferred when available (994 of the 1903 stations). For stations where only daily observations were available, the 6-hourly discharge simulations were first aggregated to daily steps (daily mean) before evaluating the objective function. The calibration was performed at the catchment level, with the 1903 selected stations covering 69.6 % of the HERA domain. A map showing the calibrated catchments is provided in Fig. S1 in the Supplement. The calibration was performed on 14 parameters that influence the modelling of snowmelt, water infiltration into the soil, surface water flow, groundwater flow, and lake and reservoir dynamics. A list of the calibration parameters is provided in Table S1 in the Supplement. Parameter values were identified using the Distributed Evolutionary Algorithm for Python (DEAP, Fortin et al., 2012) within a physically realistic range. The calibration protocol went from head catchments to downstream catchments in a top–down manner, prescribing physical dependencies between upstream and downstream catchments within the same basin.

Coastal and endorheic catchments with drainage areas smaller than 150 km², representing 6.5 % of the HERA domain, are modelled with default parameter values. Parameter values for the other ungauged catchments were estimated by parameter regionalization. These catchments are mostly located near the coastlines, with a high concentration in southern Italy and Greece, and represent 23.9 % of the HERA domain. The parameter regionalization here consists of transferring parameter values (except for the ones linked to reservoirs and lakes) from a calibrated catchment to an ungauged catchment. Catchments are matched according to climatic and geographical similarities (Beck et al., 2016). We discuss the impact of calibration on the skill of HERA in Sect. 3.1.1. For more information on the calibration of EFAS v5.0, we refer the reader to the online documentation of the Copernicus Emergency Management Service for floods (CEMS-Flood online documentation, 2023).

2.2 Climate inputs: bias-adjusted climate reanalysis data

To force the OS LISFLOOD hydrological model, we used a bias-adjusted and downscaled climate dataset based on the ERA5-Land climate reanalysis (Muñoz-Sabater et al., 2021). The main steps involved in the preparation of the climate inputs are summarized in Fig. 3. The following variables are retrieved from ERA5-Land at an hourly temporal resolution for 1950–2020:

• total precipitation (tp);

• mean temperature (ta);

• mean zonal and meridional wind speeds (u,v);

• mean dew point temperature (td); and

• total surface solar radiation downwards (ssrd).

Figure 3 Climate input pre-processing scheme, including temporal aggregation, bias adjustment, statistical downscaling, and processing of evapotranspiration.

Precipitation and temperature data were aggregated to 6-hourly resolution and the other variables to daily resolution (Fig. 3). All the variables were averaged, except for precipitation, which was summed to reach the target temporal resolution. Minimum and maximum daily temperatures were also calculated, while dew point temperature was converted into relative humidity and actual vapour pressure.

Our setting of OS LISFLOOD requires meteorological data with a 1^′ resolution. To downscale ERA5-Land data from $\mathrm{0.1}\mathit{°}={\mathrm{6}}^{\prime }$ to 1^′, we performed statistical downscaling and bias adjustment using ISIMIP3BASD v3.0.0 (Lange, 2019; Lange et al., 2024; Frieler et al., 2024). The ISIMIP3BASD method was initially developed for phase 3 of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) and aims to provide robust bias adjustment of extreme values, preservation of trends across quantiles, and a clearer separation of bias adjustment and statistical downscaling compared to its predecessors (Lange, 2019). We used the new EMO-1 gridded observational dataset (1^′ version of EMO-5; Thiemig et al., 2022) developed for the operational EFAS v5.0 as the high-resolution reference dataset. EMO-1 covers the period 1990–2020 and has also been used directly as a climate input in the calibration (Sect. 2.1.2). We used 1990–2020 as the training period for the algorithm since both datasets overlap for this period. The trained algorithm is then applied to ERA5-Land to produce high-resolution data for both the training period and 1950–1989, when high-resolution data comparable to EMO-1 are not available. The resulting climate data consistently cover 1950–2020. The ISIMIP3BASD method is applied to the following variables:

• daily mean near-surface relative humidity (hurs), obtained from actual vapour pressure (vp);

• daily and 6-hourly total precipitation (pr);

• daily total surface downwelling short-wave radiation (rsds);

• daily mean near-surface wind speed (ws);

• daily and 6-hourly mean near-surface air temperature (tas);

• diurnal near-surface air temperature range ($\text{tasrange}=\text{tasmax}-\text{tasmin}$); and

• diurnal near-surface air temperature skewness ($\text{tasskew}=(\text{tas}-\text{tasmin})/\text{tasrange}$).

Here, tasmin and tasmax are the daily near-surface air temperature minimum and maximum, respectively.

Version 3.0.0 of ISIMIP3BASD differs technically from version 2.5.0, which was used to produce the climate forcing data for phase 3b of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP3b; Frieler et al., 2024), yet both versions produce the same results, and we apply version 3.0.0 using the same climate-variable-specific parameter settings as for the ISIMIP3b data production (Lange et al., 2024; Frieler et al., 2024). ISIMIP3BASD was designed for daily data, but it is applied here to bias-adjust and statistically downscale sub-daily (6-hourly pr and tas) data as if these are daily values. For the bias adjustment, a parametric trend-preserving quantile mapping method was applied to pr, sfcwind, tas, and tasrange, while non-parametric quantile mapping was applied to hurs, rsds, and tasskew. The bias adjustment was done at the spatial resolution of ERA5-Land, 6^′, using spatially aggregated EMO-1 data (spatial averaging). Data resulting from the bias-adjustment were then statistically downscaled to 1^′ spatial resolution by using an algorithm based on the multivariate bias correction (MBCn) bias-adjustment method (Cannon, 2018) (Fig. 3). The downscaling method is conservative in the sense that the 1^′ output data would be identical to the 6^′ input data in case the former are spatially aggregated back to 6^′ resolution.

Finally, potential evapotranspiration (et₀), potential open-water evapotranspiration (e₀), and potential bare soil evapotranspiration (es₀) are computed with bias-adjusted and downscaled data at the pixel level using an approach based on the Penman–Monteith equation with the LISVAP model (LISVAP online documentation, 2023).

2.3 Surface field maps

OS LISFLOOD requires a set of surface field maps. Depending on the model set-up, it can ingest up to 108 surface fields divided in six categories:

• (i)

  catchment morphology and river networks;

• (ii)

  vegetation cover types and properties;

• (iii)

  soil properties;

• (iv)

  land use;

• (v)

  water demand; and

• (vi)

  lake and reservoir information.

The first three categories, hereafter referred to as static maps, were taken directly from the CEMS_SurfaceFields_2022 open-source dataset of the Copernicus Emergency Management Service developed for the European domain at 1 arcmin resolution, which can be found in the JRC Data Catalogue (Choulga et al., 2024). The last three categories were derived from CEMS_SurfaceFields_2022 and modified to take into account socio-economic changes (hereafter referred to as dynamic socio-economic maps). This section briefly presents each of the map categories, with an emphasis on dynamic socio-economic maps, which are original to this work.

2.3.1 Static maps

Static maps include surface fields of morphology and channel shapes (14 maps), vegetation properties (18 maps), and soil properties (29 maps).

Morphology and river network information were used directly for the computation of snowmelt, temperature scaling, river routing, and open-water evapotranspiration. Morphological information was derived from elevation and includes the elevation gradient, the within-grid standard deviation of elevation, and Manning's roughness coefficient. Maps representing channel shapes and river networks provide information on the grid cell area (which varies with latitude as the grid projection is WGS84), local drainage direction, upstream area, and channel dimensions. All the morphologies and river network maps were derived from Multi-Error-Removed Improved-Terrain Digital Elevation Model v.1.0.3 (MERIT DEM) (Yamazaki et al., 2019) and Catchment-based Macro-scale Floodplain (CaMa-Flood) Global River Hydrodynamics Model v4.0 maps (Yamazaki, 2023).

Vegetation cover types and property maps are involved in the computation of precipitation interception, evaporation, transpiration, surface runoff, and root water uptake. These properties are described by four variables: crop coefficients (transpiration), crop groups (water uptake), Manning roughness (surface runoff), and leaf area index (interception and evaporation). Each of these variables was mapped for three different land cover types: forest, irrigated, and other. Further, maps of planting and harvesting days for rice, which has specific water demands, are also available. Vegetation properties were derived from several data sources, including the Copernicus Global Land Service (CGLS) leaf area index (LAI) at 1 km (Copernicus, 2021), the Spatial Production Allocation Model (SPAM) Global Spatially-Disaggregated Crop Production Statistics Data for 2010 (Yu et al., 2020; International Food Policy Research Institute, 2019), and the Food and Agriculture Organization (FAO) of the United Nations Irrigation and Drainage Paper No. 56 (Allen et al., 1998).

Soil properties refer to physical characteristics of the soil and aim to describe the water dynamics through a vertical soil profile. In OS LISFLOOD, the soil profile is composed of three layers: superficial (0–5 cm), upper (5–varying (30–50) cm), and lower soil layers. For each layer, variables representing soil hydraulic properties (e.g. soil moisture content or pore size index) are provided. Similarly to vegetation property maps, variables were mapped for two categories of land cover, “forest” and “other”. Soil properties were derived from the International Soil Reference and Information Centre (ISRIC) global gridded SoilGrids dataset (release 2017) available at 250 m (Hengl et al., 2014), which is based on more than 150 000 observation sites and covariate data.

A table summarizing all the static and dynamic surface field maps used to produce HERA is provided in Table S2 in the Supplement. For more details on these surface field maps, their production, and the input datasets used, we refer the reader to Choulga et al. (2024).

2.3.2 Dynamic land use

OS LISFLOOD includes six land use classes as inputs: rice, other irrigated land, forest, sealed surfaces, open water, and other (non-irrigated agriculture, non-forest natural, pervious artificial); these land use classes are mostly based on the CLC-Refined 2006 dataset of Batista e Silva et al. (2013) in the default setting. Of the hydrological processes, interception, evapotranspiration, infiltration, and surface runoff respond differently to each land use type. With the aim of better representing complex rainfall-runoff processes, OS LISFLOOD accounts for the sub-grid variability in land use. Therefore, the spatial distribution of each land use class is defined as a percentage of the whole represented area of a given pixel (European Commission, Joint Research Centre (JRC), 2025a). The magnitude of the variation of hydrological responses is tied to the magnitudes of the changes in land cover. De Roo et al. (2001), for instance, investigated the effects of land use changes on floods in two European catchments and identified different results depending on the magnitude of the land cover change. While such changes tend to have a limited impact on river discharge, they can locally increase flood magnitude (Merz et al., 2021; Sajikumar and Remya, 2015; Van Lanen et al., 2013; Van Loon, 2015). Here we modified the grid cell fractions of each land use class using HANZE-Exposure land use maps at 100 m resolution (Paprotny and Mengel, 2023) for 42 countries in the study area. In the remaining part of the domain, we used coarser, 5^′ resolution maps from HYDE 3.2 (Klein Goldewijk et al., 2017) to modify the 2006 values. The temporal evolution of the land area of each class is displayed in Fig. 5a. There has been a strong increase in sealed surfaces (+40 %), while for the other relevant land use classes the changes are less than 10 %, with more land occupied by irrigated agriculture (except rice), water surfaces (due to reservoir construction), and forests.

Figure 4 Variation in socio-economic inputs in the hydrological model, averaged over the entire EFAS domain: (a) land area by use category, 1950=100; (b) number of existing reservoirs; (c) water demand by sector in millimetres per grid cell per year; and (d) shares of land use between the different classes in 2020.

Figure 5 HERA hydrological skill for the 2448 selected stations in terms of (a) KGE^′ and its three components: (b) Pearson correlation, (c) bias ratio, and (d) variability ratio. In panel (a), the green dashed vertical line represents the benchmark KGE^′ value (−0.41). The red vertical line represents the ideal values and the blue dot represents the median for all the stations.

2.3.3 Dynamic water abstraction

Human water use, representing water withdrawal from the natural environment (e.g. rivers, reservoirs, or groundwater) for human needs, is grouped into four main sectors: livestock, domestic, manufacturing industry, and energy production. In OS LISFLOOD, water use is supplied by surface water bodies and groundwater, depending on the sector (Choulga et al., 2024). A considerable increase in water abstraction in a region can diminish surface water resources within the same area. The model also accounts for groundwater abstraction for human use, except for flooded irrigation and cooling processes. Increased groundwater abstraction can locally reduce (or halt) baseflow. To derive monthly historic sectoral water withdrawal maps, we followed the methodology of Huang et al. (2018) and used the FAO AQUASTAT sectoral water withdrawal data (Food and Agriculture Organisation, 2023) as a starting point. These data were subsequently spatially and temporally disaggregated using a variety of datasets. These include the Global Human Settlement Layer (Schiavina et al., 2019; Florczyk et al., 2019) for population estimates, the Global Change Analysis Model (GCAM; Calvin et al., 2019) for regional water withdrawal and electricity consumption, and the Gridded Livestock of the World (GLW; Gilbert et al., 2018) for livestock distribution. Additional datasets included Multi-Source Weather (MSWX; Beck et al., 2022) for air temperature data, United States Geological Survey (USGS) water withdrawal estimates, and Vassolo and Döll (2005) industrial and thermoelectric withdrawal maps. More information on the water demand and input datasets used is provided in Choulga et al. (2024).

We extrapolated the water withdrawal maps to the period 1950–1978 using annual gridded 0.5° data from ISIMIP 3a (Frieler et al., 2024; Wada et al., 2016) that were downscaled to 1^′ resolution using historical population data from HANZE (Paprotny and Mengel, 2023) and HYDE 3.2 (Klein Goldewijk et al., 2017) for the other parts of the domain. More precisely, the ratio between EFAS high-resolution water demand maps and the ISIMIP 3a dataset for 1979 was used to adjust the water withdrawal data in each grid cell. Intra-annual (monthly) cycling of water use in the energy and domestic sectors was estimated for 1950–1978 using the same approach as for 1979–2020, informed by temperature data from our input meteorological dataset (Sect. 2.3.1). Livestock water use was assumed to be constant before 1979. Water demand and use for irrigation were computed directly by the hydrological model based on land use data and available water. The evolution of water use by sectors between 1950 and 2020 is displayed in Fig. 5c and Table S4 in the Supplement. The total water use peaked in 1990 after more than doubling since the 1950s and before declining due to a drop in demand from the manufacturing and energy sectors. Nonetheless, there are usually much stronger trends at the country or catchment levels.

2.3.4 Dynamic reservoir maps

Reservoir maps contain locations and identifiers of reservoirs and are linked to tables containing metadata on storage capacity, construction year, and a set of values associated with reservoir operation rules. Normal reservoir outflow rates were further adjusted through the model calibration (Sect. 2.1.2). The year of construction for each reservoir was taken from the EFAS reservoir database HANZE (Paprotny and Mengel, 2023), the Global Reservoir and Dam Database (GRanD) v1.3 (Lehner et al., 2011), or additional manual research for reservoirs not covered by the three datasets. The reservoir maps are updated every simulation year (1 January) by adding newly built reservoirs. When a reservoir is added, it is considered empty and fills up according to its associated metadata. Figure 5b shows the evolution of the number of reservoirs in Europe during the period 1950–2020. The number of reservoirs in the model increased 6-fold from 244 in 1950 to 1419 in 2020, though few were built after the late 1980s.

3 Results

3.1 Technical validation

We evaluated our hydrological reanalysis by comparison against a dataset of daily river discharge observations from 3442 stations across Europe. Of the data obtained, 60 % were from the Global Runoff Data Centre (GRDC) and 40 % from national public datasets of France, Norway, Poland, Spain, Sweden, and the UK. Furthermore, this dataset was compiled independently of the one used in the EFAS calibration (Sect. 2.1.2). The stations' record durations vary between 1 and 71 years. The selection of the stations used for validation is based on several criteria:

• Spatial matching. To link stations to their corresponding river pixel, we scanned the nine modelled pixels around the river gauge location. When information on the upstream area was available (for 60 % of the stations), we retained the pixel with the upstream area closest to the reported one. For pixels without information on the upstream area, we retained the one with the simulated mean discharge (Q_mean) closest to the observed one. For more accurate spatial matching, we used the available OS LISFLOOD coordinates from the EFAS calibration (1026 stations). A total of 546 stations did not match OS LISFLOOD river pixels, mostly due to their upstream area being lower than 100 km².

• Upstream area verification. The spatial matching selected the closest upstream area for stations where we have information on catchment area. It is however possible that the reported catchment differs greatly from its matched pixel upstream area. We removed stations where the difference between the pixel and observed upstream areas was larger than 50 % (51 stations).

• Mean discharge comparison. For some stations, the ratio between the observed and simulated Q_mean was suspicious. This could be due to an erroneous spatial match (i.e. matching of a river with a station on a tributary). As uncertainty grows with smaller streams, we decided to remove those with a suspicious Q_mean ratio (r_Qmean>6 or r_Qmean>3 if Q_mean,obs>10 m³ s⁻¹) (49 stations).

• Manual check. A manual verification was performed on 66 stations with ${\text{KGE}}^{\prime }<-\mathrm{0.41}$. Each station and its matching pixel were checked individually, resulting in the removal of 13 more stations due to wrong spatial matching, erroneous station location, and doubtful observations. The corresponding river pixel was set manually for eight stations. Manually checked stations and the reason for their exclusion or inclusion are provided in Table S5 in the Supplement.

• Finally, we removed stations with a record length shorter than 30 years (334 stations). This enabled a meaningful comparison between different locations in the validation process.

This procedure resulted in the selection of 2448 river stations across Europe, with an upstream area ranging from 100 to 785 421 km². Of these stations, more than half (1507) have an upstream area of less than 1000 km², and one-fifth (498) have an upstream area of less than 200 km².

HERA comes at a sub-daily resolution (6-hourly), but the performance could only be evaluated at the daily time step of the observational dataset. Discharge data from HERA were therefore aggregated (daily mean) for the technical validation. We expect the performance to be slightly higher at the daily scale, as the temporal aggregation tends to increase the correlation between observed and modelled discharge. The performance was assessed using the KGE^′ on discharge data (Gupta et al., 2009; Kling et al., 2012). KGE^′ was used as the standard performance metric in EFAS and GLOFAS (Harrigan et al., 2020; Cammalleri et al., 2020a) as well as in other hydrological model assessments (Lin et al., 2019; Harrigan et al., 2020; Beck et al., 2017) and has three components: correlation, bias errors, and variability errors:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\text{KGE}}^{\prime }=\mathrm{1}-\sqrt{(r-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\gamma }-\mathrm{1}{)}^{\mathrm{2}}},\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\beta }={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{s}}}{{\mathit{\mu }}_{\mathrm{o}}}},\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\gamma }={\displaystyle \frac{\frac{{\mathit{\sigma }}_{\mathrm{s}}}{{\mathit{\mu }}_{\mathrm{s}}}}{\frac{{\mathit{\sigma }}_{\mathrm{o}}}{{\mathit{\mu }}_{\mathrm{o}}}}},\end{array}$$

where r is the Pearson correlation coefficient between the simulated (s) and observed (o) flow, β is the bias ratio, γ is the variability ratio, μ is the mean discharge, and σ is the discharge standard deviation. KGE^′ and its three components are dimensionless, with an optimal value of 1. It is important to note here that KGE^′ values should not be interpreted like the more traditional Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970). Indeed, for KGE^′ the mean flow benchmark has a value of ${\text{KGE}}^{\prime }=\mathrm{1}-\sqrt{\mathrm{2}}=-\mathrm{0.41}$. Any value above −0.41 therefore exceeds the benchmark (Knoben et al., 2019), meaning that the model performs better than simply taking the mean.

In Sect. 3.1.1, we assessed the model performance across space, time (1951–2020), and catchment size in order to identify the strengths and weaknesses of HERA. Despite covering many aspects of the performance of hydrological models, KGE^′ mainly focuses on mean values and gives a higher weight to high extremes compared to low ones. As this dataset also aims for use in long-term analysis of hydrological extremes, we evaluated how well high and low extremes are reproduced, including their timing and seasonality.

3.1.1 Hydrological performance

Here we quantified the overall performance of HERA in terms of KGE^′ as well as the decomposition of this indicator into its three components: correlation, bias, and variability. Figure 5 displays the distribution of KGE^′ and its three components across the 2448 validation stations. We obtained ${\text{KGE}}^{\prime }>-\mathrm{0.41}$ for 2411 (98.5 %) of them, meaning that the reanalysis is skilful for these stations (Fig. 5a). The median KGE^′ across all the catchments is 0.55, while the mean is 0.46, although this value varies widely across the catchments (Figs. 5a and 6a). The mean correlation value is relatively high ($\overline{r}=\mathrm{0.69}$), with 90 % of the stations having r>0.5 (Fig. 5b). From Fig. 5c and d, we can observe that there is a tendency to slightly underestimate flows ($\overline{\mathit{\beta }}=-\mathrm{13.1}$ %) and flow variability ($\overline{\mathit{\gamma }}=-\mathrm{14.2}$ %). The bias ranges between 0.8 and 1.2 (0.5–1.5) in 50 % (91 %) of the river gauges, which is considered very good for hydrological reanalysis (Harrigan et al., 2020; Alfieri et al., 2020; Lin et al., 2019; Yang et al., 2021).

Figure 6 KGE^′ and its three components: (b) Pearson correlation, (c), bias ratio, and (d) variability ratio at the 2448 river gauges considered in the validation of HERA. Point sizes are proportional to catchment sizes.

Figure 6 shows the spatial performance of the model in terms of KGE^′ and its components. The highest skill can be observed in central and north-western Europe. The vast majority of stations in the UK, Germany, France, Austria, and Switzerland (which together account for 51 % of all 2448 stations) exhibit a good (>0.5) to very good (>0.75) KGE^′. On the other hand, performance is relatively poor in Spain, Cyprus, Scandinavia, and northern Poland. Factors that can explain the poor performances in southern Europe include the combination of arid climates and the strong influence of lakes and reservoirs (Fig. 7c). Dry catchments where precipitation events are separated by long dry spells are in general very difficult to model (Cantoni et al., 2022). In Scandinavia, the negative bias (Fig. 6c) could be linked to an underestimation of precipitation and snowmelt in Scandinavian mountains (Beck et al., 2017, 2020). Figure 6d presents the variability ratio of simulated to observed flow. Overall, our reanalysis exhibits lower variability than observations, with 83 % of the catchments having a variability ratio below 1. The underestimation of variability was also found in the EFAS v5.0 run, although it is more pronounced in HERA. This could be explained by the different meteorological forcings used in the two runs.

Figure 7 Boxplot of the HERA KGE^′ according to different classifications of the 2448 river stations used in the validation: (a) time, (b) catchment area, and (c) reservoir impacts. The numbers inside the boxplot represent the number of river gauges for each category, while the colours of the boxplot represent the median performance of the group from low (light blue) to high (dark blue).

We validate HERA on stations with a wide range of catchment areas (mean upstream area of 7615 km²), which has an impact on OS LISFLOOD performance (Harrigan et al., 2020). The set of 2448 validation stations includes stations that were used in the calibration process (596) as well as stations in uncalibrated catchments (36) (Fig. S1). In Fig. 7, we break down the performance of the reanalysis according to the different attributes of each catchment: time (Fig. 7a), catchment area (Fig. 7b), reservoir impact (Fig. 7c), and calibration status (Fig. 7d).

Overall, the skill of HERA shows a slight increase through time of 21 % in ${\text{KGE}}_{\text{med}}^{\prime }$ between the 1950s and the 2010s. The skill increases between 1951 and 1980 and then stabilizes from 1981 to 2020, though the results are influenced by changes in the gauge data availability over time. This could also be driven by improved climate inputs. Figure 7b shows that the model skill increases with catchment size, from ${\text{KGE}}_{\text{med}}^{\prime }$ values of 0.44 (inter-quartile range IQR 0.25–0.59) for the 498 smallest catchments (<200 km²) to 0.77 (IQR 0.68–0.84) for the 28 largest catchments (>100 000 km²). Such patterns have already been observed at global scales (Harrigan et al., 2020). It is important to note here that the majority of stations used in this validation (62 %) have an upstream area below 1000 km², and the median upstream area of the 2448 stations is 583 km². This is half of the median upstream area of the 1903 stations used in the calibration of EFAS v5.0 (CEMS-Flood online documentation, 2023).

We also divided stations according to reservoir influence. From the 1420 reservoirs active in 2020 (which represent the maximum amount over the considered time window), we estimated the impact of reservoirs on streamflow at the grid cell level. This was done by computing the ratio (c; –) of reservoir volume to mean discharge (Nilsson et al., 2005) at every grid cell. The ratio has been computed with the accuflux function from PCRaster and compares the upstream cumulative reservoir capacity (m³) and the cell-specific annual volume of annual streamflow (m³) (Zajac et al., 2017). This ratio varies between 0 and 1608 downstream of Embalse de Finisterre in central Spain. Most of the river grid cells highly impacted by reservoirs are found in southern Europe, particularly in Spain and Bulgaria. Figure 7c highlights the influence of reservoirs on the skill of the reanalysis. The river cells affected (medium and high, c>0.5) only represent 6 % of the stations and grid cells in the domain (Fig. 2). The median skill is lowest for highly impacted (c>1) stations with ${\text{KGE}}_{\text{med}}^{\prime }=\mathrm{0.24}$, whereas minimally impacted stations have a ${\text{KGE}}_{\text{med}}^{\prime }$ of 0.55. This highlights the difficulty of large-scale hydrological models such as OS LISFLOOD in accurately simulating reservoir outflows (Zajac et al., 2017).

Finally, we investigated the influence of calibration on model skill. In Fig. 7d, river gauges are divided into four groups according to their calibration status. As displayed in Fig. S1, 83 % of the stations considered in the validation fall into the domain calibrated for EFAS v5.0. We find better performance for calibrated stations (${\text{KGE}}_{\text{med}}^{\prime }=\mathrm{0.64}$) and comparable skill for stations within the calibrated domain (${\text{KGE}}_{\text{med}}^{\prime }=\mathrm{0.52}$) and stations benefitting from the parameter regionalization (${\text{KGE}}_{\text{med}}^{\prime }=\mathrm{0.47}$). The performance is much worse for catchments with default parameters, which here are limited to small (<150 km²) coastal and endorheic catchments.

3.1.2 Reproduction of extremes

Large-scale hydrological models forced by climate reanalysis often fail to reproduce extreme hydrological event characteristics, in part due to the coarse spatial and temporal resolution (Brunner et al., 2021; McClean et al., 2023). Here, we analyse how well HERA reproduces different flow quantiles (q05, median, and q95) through the Pearson correlation coefficient and the coefficient of determination (R²) (Fig. 8) for the 2448 considered catchments. The ability to capture annual maxima and minima as well as their seasonality is also assessed (Fig. 9).

Figure 8 Scatterplot of observed and simulated river flow quantiles (m³ s⁻¹): (a) 5th quantile, (b) median (q50), and (c) 95th quantile (q95) for the 2448 river gauges.

Figure 9 Assessment of the ability of HERA to reproduce the timing of annual maximum and minimum flows. (a) Violin plot of errors in the mean day of occurrence of annual maxima (daily discharge) and minima (30 d averaged discharge) computed with circular statistics. Inside each violin plot, boxplots display the median and 1st and 3rd quantiles. (b) Difference between the modelled and observed mean annual maximum dates (a positive value means a later occurrence in HERA). (c) Difference between the modelled and observed mean annual minimum dates (a positive value means a later occurrence in HERA).

Figure 8 displays scatterplots of observed and simulated quantiles. Each point represents 1 of the 2448 stations. We observe that low (5th quantile: Q₀₅) and median (Q₅₀) flows are generally well represented with R²>0.99 (Fig. 8a and b), especially for higher discharge values. However, despite this generally good agreement, there is a more pronounced deviation of simulated values from observations for lower flow values, expressed by a higher dispersion for Q₀₅. These deviations can be attributed to bias in climate inputs (McClean et al., 2023), the hydrological model (Feyen and Dankers, 2009), errors in flow measurements (especially for Q₀₅) (Despax, 2016; Tomkins, 2014), and anthropogenic impacts on low- and median-flow regimes (Brunner, 2021) that are not accurately represented in the model (Fig. 8c). The number of stations with large deviations in the reproduction of high-flow statistics (Q₉₅) is minor compared to Q₀₅ and Q₅₀. Nonetheless, despite a relatively high R² (0.99), there is a general underestimation in the simulated values (Fig. 8c), which is common for large-scale hydrological models. Similarly to low and median flows, errors in high-flow statistics can be due to biases and smoothing of extremes in climate inputs and errors in the hydrological modelling. Uncertainty associated with flow measurements also plays a major role in high flows, as river discharges are usually not measured directly during floods (Despax, 2016). Finally, the spatial and temporal resolution of the model can affect its ability to reproduce high flows, particularly for flash floods in small catchments.

We also assessed the ability of the reanalysis to reproduce the timing of annual maxima and minima of discharge as well as their overall seasonality. As the daily temporal scale is not the most relevant factor when it comes to drought analysis with discharge data (Hannaford and Marsh, 2006; Kohn et al., 2019), annual minima were computed from 30 d moving average flows. Figure 9a displays the mismatch in the mean day of occurrence computed with circular statistics following Berghuijs et al. (2019). We observe that the median error in the mean day is very close to 0 for both maxima ($\text{median}=\mathrm{0.1},\text{IQR}=-\mathrm{12}$–18) and minima ($\text{median}=-\mathrm{1},\text{IQR}=-\mathrm{28}$–41) but with a much higher dispersion for annual minima compared to annual maxima. The higher dispersion for low flows is due to the slow-onset nature of these events (Brunner, 2021). Figure 9b shows the difference in timing between simulated and observed annual maxima across the 2448 considered stations. Differences in timing are smaller over the Atlantic coast, though a particularly high lag (simulated maxima delayed by 30 d or more in HERA) is observed over Poland and central Spain. For low flows (Fig. 9c), delays in central Europe are more than 30 d, while in Scandinavia the timing can be up to several months too early. This can be explained by the high hybridity of river regimes (several high- and low-flow seasons) in these regions, which may be captured with varying accuracy in HERA.

In addition to the validation protocol presented in this section, we compared the reported performances of HERA with other recent hydrological datasets and carried out a comparison between HERA and another recent hydrological simulation done with the grid-based conceptual mesoscale Hydrological Model (mHM) (Kumar et al., 2013; Samaniego et al., 2010, 2019; Thober et al., 2019) for Europe for the period 1960–2010. More details on the comparison are provided in Figs. S3–S6 in the Supplement.

3.2 Usage notes

HERA brings together several improvements (climate, scale, and socio-economic dynamics) to better simulate river discharge in catchments of Europe over the past 70 years. Despite still covering a relatively short period of time compared to human history on Earth, these 70 years capture a very intense period of climate and socio-economic change, often called the Anthropocene, and offers multiple research opportunities:

• assessment of long-term trends in European river regimes;

• provision of benchmark data for “data-poor” areas;

• generation of catalogues of flood and drought events;

• identification of spatial and temporal correlations between European catchments;

• identification of changes in hydrological extreme characteristics (frequency, magnitude, and timing);

• combination with other data products for compound hazard analysis; and

• provision of scenarios for flood inundation simulations.

In this section, we briefly present possible usage of the data, addressing changes in regime for diverse rivers across Europe (Fig. 10).

Figure 10 Changes in flow regime between 1951–1981 (blue) and 1990–2020 (red) for six diverse European rivers: (a) Ardèche, (b) Schelde, (c) Ebro, (d) Rhône, (e) Danube, and (f) Vistula. The regime is computed here as the 30 d moving average. Shaded coloured areas represent the IQR of discharge for every day of the year. The grey-shaded area represents the absolute difference between the two regimes corresponding to different periods.

Figure 10 displays hydrological regimes, represented here as the mean of a 30 d average moving window over a given period, for six European rivers. These rivers differ in terms of hydrological regimes, with three main regimes represented:

• a Mediterranean pluvial regime for the Ardèche (a), with its recognizable high flows in autumn;

• a pluvial or oceanic regime for the Schelde in Ghent (b) and the Ebro in Zaragoza (c); and

• a nival regime for the upper Rhône in Lyon (d), the Danube in Vienna (e), and the Vistula in Warsaw (f).

These six rivers also vary in terms of catchment area, geographical location (France, Austria, Poland, Belgium, or Spain), climate (Mediterranean, continental, oceanic, or alpine), and geomorphological conditions. For each river, the flow regimes for 1951–1981 (in blue, first 30 years of HERA) and 1990–2020 (in red, last 30 years of HERA) are shown. By comparing the two regimes, one can observe diverging patterns of changes between these rivers. For the two pluvial rivers, the Schelde and the Ebro (Fig. 10b and c), we observe opposite patterns of change: the Schelde saw an increase in its average discharge throughout the year, while the Ebro experienced a downward shift in its regime. For the upper Rhône and Danube (Fig. 10d and e), which are influenced by snowmelt in their upper catchments, we see lower and earlier flow peaks in spring and summer. The Vistula (Fig. 10f) saw an overall increase in flow throughout the year. Finally, the Ardèche (Fig. 10a) has seen reduced flow throughout the year, with a notable decrease in late winter which can be associated with the reduction in snowfall in the Massif Central, where the Ardèche has its upper waters (François et al., 2023). The timing of the autumn peak seems to have shifted slightly towards earlier dates, in agreement with a recent study on trends in Mediterranean floods (Tramblay et al., 2023).

4 Discussion

Recent developments in diverse fields, including climate, hydrology, remote sensing, and computational sciences, have made the generation of high-resolution reanalysis products possible (Aerts et al., 2022; Hanasaki et al., 2022; Hoch et al., 2023). In this context, HERA brings discharge data for all European rivers with upstream areas larger than 100 km² for the period 1951–2020. With its refined spatial and temporal resolution, HERA represents hydrological processes in Europe with more detail than previous publicly available hydrological reanalysis products (Harrigan et al., 2020; Schellekens et al., 2017). Calibrating hydrological models can significantly improve river flow simulation (Beck et al., 2017; Kauffeldt et al., 2016). Parameters in 93.5 % of the HERA domain were adjusted with a calibration process (Sect. 2.1.2) or by parameter regionalization (Beck et al., 2016). This is a very high calibration coverage for a GHM (Beck et al., 2017), which can be explained by the relatively high coverage in river gauging stations in Europe.

Table 1 Description of the HERA dataset.

Table 2 List of inputs and outputs of OS LISFLOOD provided in the HERA database (link here).

It is difficult to compare HERA with other recent hydrological reanalyses such as GLOFAS-ERA5 (Harrigan et al., 2020) and GRFR (Yang et al., 2021), for several reasons: (i) spatial coverage (global vs. continental), (ii) spatial resolution (0.25°, 0.05°, 0.0167°), (iii) temporal coverage, and (iv) dynamic vs. static socio-economic conditions. However, we provide a short summary of the reported performances of HERA, GLOFAS-ERA5, GRFR, and a European-scale hydrological simulation with the mHM model (EUmHM) in Table S6 in the Supplement. While the reported performances of HERA are higher than its global counterparts, they are very close to the performance of EUmHM. In a more detailed comparison with EUmHM over 515 European river gauges (Figs. S3–S6), we show that HERA generally outperforms the EUmHM run in terms of KGE^′ (Fig. S4), but both models exhibit strengths and weaknesses, spatially (Fig. S5) and in terms of the components of KGE^′ (Fig. S6). Differences in performances between the HERA and EUmHM runs can be attributed to the many different features in the two runs, such as meteorological forcing, resolution, calibration, and flow routing within the hydrological model. Conversely, HERA shares a great number of features with the EFAS v5.0 reanalysis (Decremer et al., 2023), with a slightly lower performance (not shown here). Nonetheless, EFAS v5.0 only covers the period 1990–2022 and assumes static socio-economic conditions (land use, water abstraction, and reservoirs).

Similarly to the other aforementioned hydrological reanalyses, HERA exhibits reduced performance in cold and semi-arid catchments. This could be related to deficiencies in the representation of snow processes within OS LISFLOOD or the underestimation of precipitation at northern latitudes (Beck et al., 2017, 2020). Semi-arid environments are notoriously challenging areas for hydrological models due to the highly non-linear rainfall-runoff response and lower precipitation data quality (Cantoni et al., 2022). GHMs tend to represent runoff poorly in small- to medium-sized catchments (10–10 000 km²) (Harrigan et al., 2020; Sood and Smakhtin, 2015), and nearly 90 % of the catchments used in the validation of HERA (Sect. 3.1) are small- to medium-sized catchments. The drop in performance with smaller catchment areas in HERA remains, however, moderate compared to that of the GLOFAS-ERA5 global hydrological reanalysis (Harrigan et al., 2020). The presence of reservoirs also influences the performances of the reanalyses. While including reservoirs in the hydrological modelling has a positive impact on model performance (Zajac et al., 2017), there is still a high level of uncertainty regarding the operating rules of each reservoir. Moreover, the 1422 reservoirs used to generate HERA most likely represent just fractions, mainly the largest ones, of all the operational reservoirs in the modelled domain (Speckhann et al., 2021). In summary, the main strength of HERA lies in its relatively low bias in comparison to the other hydrological datasets considered here (Table S6 and Fig. S6), while its performances are hampered by its underestimation of variability.

HERA is generated through hydrological modelling, which brings a range of uncertainties that can be divided into four categories: (i) model inputs, (ii) model structure, (iii) parameter values, and (iv) observations. It remains challenging to quantify these uncertainties. However, the quality of the inputs and in particular the climate inputs is often an important factor in uncertainty (Beck et al., 2017; Sood and Smakhtin, 2015). Despite efforts in bias correction and downscaling of the climate input, it seems that, on average, HERA slightly underestimates river discharges, with a more pronounced bias for high flows. As reported in other studies, negative biases can be related to an underestimation of precipitation in the climate inputs, in particular for extreme events (McClean et al., 2023; Mahto and Mishra, 2019), at high latitudes, and in (semi-)arid catchments (Beck et al., 2016; Sood and Smakhtin, 2015; Hirpa et al., 2018). Model structure can also play an important role, as shown in Fig. S6, where EUmHM is the best model in terms of correlation, while HERA exhibits a smaller bias ratio. This can be the result of different choices made in the main equation behind the two models, resulting in different responses to forcings and calibration. The large impact of model selection on streamflow and trend estimates is now increasingly acknowledged (Karlsson et al., 2016; Clark et al., 2016). Calibration generally improves streamflow simulations (Hirpa et al., 2018), and HERA also shows better performance for the stations used in the calibration process (Fig. 7d). The negative biases and variability ratios can be related to the different meteorological forcing (EMO-1) used in the calibration, although an underestimation of the variability was also found in the EFAS v5.0 run (which is forced by EMO-1). The method, parameters, and skill metrics used for calibration further affect the uncertainties. Despite its qualities, the skill metric used for the calibration presented in Sect. 2.1.2 (KGE^′) is known to result in an underestimation of variability (Brunner et al., 2021) and to place more weight on high values (Garcia et al., 2017).

This could partly explain the reduced performances in reproducing the extreme low flows observed in Figs. 8 and 9. Other uncertainties can arise from surface field maps (Sect. 2.3) and measurement of river discharges (instruments and rating curves). With sparser gauging and more complex hydraulic conditions for high and low flows, uncertainty rises (Despax, 2016).

5 Code and data availability

The HERA hydrological reanalysis and its climate and dynamic socio-economic inputs are available in the JRC data catalogue: https://doi.org/10.2905/a605a675-9444-4017-8b34-d66be5b18c95 (Tilloy et al., 2024). Table 1 provides a brief description of the dataset, and Table 2 gives a general overview of the content of the dataset.

The dataset consists of three distinct folders that are described here and in Table 2:

• Climate inputs. This is a folder containing the climate forcing for the OS LISFLOOD hydrological model. Out of the five variables provided, three are at daily temporal resolution, i.e. potential evapotranspiration, potential evaporation, and potential evaporation from bare soil (obtained with LISVAP; LISVAP online documentation, 2023), while two have a 6-hourly time step, precipitation, and temperature. The spatial resolution of the climate inputs is 1^′. The files are in NetCDF format, with one file per year per variable for a total of 355 files (2.3 TB of data).

• Socio-economic inputs. This is a folder containing the dynamic surface field maps (Sect. 2.3), divided into three categories: land use, reservoirs, and water demand. The land use sub-folder contains 426 yearly files (4.6 GB) of land use fraction maps for all six land use classes. The reservoir sub-folder contains 71 yearly files (3.6 GB) of reservoir locations and identifiers. Reservoirs are added or discarded from the simulation every year according to their construction or destruction data. Finally, the water demand sub-folders contain four files (3.9 GB) representing water demand for the considered sectors (Sect. 2.3.3). Each file contains monthly maps of water abstraction for a given sector. All socio-economic inputs are provided in NetCDF format.

• River discharge. This folder contains river discharge NetCDF files for each year at a 6-hourly time step for all European rivers with an upstream area greater than 100 km² (2.3 GB per file, 166 GB total).

All the data share the same projection (WGS84) of grid and spatial resolution (1^′). Static surface field maps were retrieved directly from the OS LISFLOOD static and parameter maps for Europe (2024) dataset, which were developed in the context of the new EFAS deployment (Decremer et al., 2023). It is important to note that HERA simulates discharge in a slightly smaller domain than the original EFAS domain, and the mask used for HERA is also provided in the dataset.

The LISFLOOD OS hydrological model used in this research is released as open-source software and is available at https://ec-jrc.github.io/lisflood/ (European Commission, Joint Research Centre (JRC), 2025a). Version 4.1.2 of the code was used in this study (https://github.com/ec-jrc/lisflood-code/tree/v4.1.2; European Commission, Joint Research Centre (JRC), 2025b). The R and Python scripts used to assess the performances of HERA against observed discharge and to generate the figures of this article are available in this GitHub repository (https://github.com/Alowis/HERA, last access: 20 January 2025) or from the following Zenodo repository: https://doi.org/10.5281/zenodo.14718275 (Alowis, 2025).

6 Conclusions

Despite the limitations discussed above, HERA represents a state-of-the-art, high-resolution, and long-term hydrological reanalysis for Europe in the form of homogeneous river flow data generated with the OS LISFLOOD model. To our knowledge, no other publicly available hydrological reanalysis currently provides discharge data at similar scales and spatio-temporal coverage for Europe. The inclusion of dynamic socio-economic conditions provides a more realistic reanalysis of river flows in heavily managed European catchments. The increased spatial resolution improves the performance due to a better representation of hydrological processes and inputs required to simulate them, including the river network (Hoch et al., 2023; Thober et al., 2019). HERA advances the reanalysis of extreme hydrological events, notably by the sub-daily temporal resolution and high-resolution bias-corrected climate input. The magnitude and seasonality of extremes are fairly reproduced, even if biases exist in some regions (e.g. central Poland or southern Spain). The dataset covers 70 years and is therefore suited for the analysis of long-term trends of several hydrological signatures. The modelling framework developed here further forms a basis for creating alternative (counter-factual) time series of river discharges where climatic or socio-economic conditions can be kept static, enabling the attribution of changes in hydrological regimes across Europe (Kreibich et al., 2019; Sauer et al., 2021; Scussolini et al., 2024).