High‐resolution simulated water balance and streamflow data set for 1951–2020 for the territory of Poland

Long‐term simulated water balance and streamflow data for the territory of Poland with high spatio‐temporal resolution were reconstructed for the period 1951–2020 with a daily time step, using the Soil and Water Assessment Tool (SWAT) hydrological model. The spatial structure includes 4,381 sub‐basins and 47,725 hydrological response units. The SWAT model was extensively calibrated and validated following a two‐stage, multi‐site calibration method for 88 gauging stations across Poland. The simulations show satisfactory performance with a median Kling–Gupta efficiency of 0.73 and 0.7 in the calibration and validation period respectively. The model performance exceeded that of its predecessor, the SWAT model of the Vistula and Odra basins. The simulated hydrological data set covers daily water balance components at the sub‐basin level, including precipitation total, snowfall, snowmelt, evapotranspiration (actual and potential), surface runoff, soil water, percolation, groundwater recharge, groundwater flow, tile flow and water yield (all components in mm), and daily streamflow (m3/s) at the reach level. The data set is stored in an online, publicly available research data archive 4TU.Research Data (Marcinkowski, P., Kardel, I., Płaczkowska, E., Giełczewski, M., Osuch, P., Okruszko, T., Venegas‐Cordero, N., Ignar, S. & Piniewski, M (2021). PL‐SWAT‐51_20: a high‐resolution simulated water balance and streamflow data set for 1951‐2020 for the territory of Poland, 4TU.ResearchData. dataset, https://doi.org/10.4121/16843183).


| INTRODUCTION
Freshwater resources substantially vary in time and space. They depend on the local climate, topography and geology. Their quantitative assessment is fundamental from the perspective of socio-economic development and functioning of the environment (Xu & Singh, 2004). Throughout decades, the availability of water resources has been shaped by natural factors such as climatic variability, and anthropogenic activities such as water consumption, allocation and storage. Understanding the effects of past or ongoing climate and socio-economic changes on water resources for a specific region or subbasin requires the estimation of a long-term water budget, including different user groups, that is, agricultural, municipal or environmental water users. Depending on the scale of the analysis, different tools and data sets can be used. Although popular, the flow gauge-based approach assessing water budget by means of streamflow records might be inapplicable at a large (national) scale due to low density of gauging stations or missing data records. At the national scale, water resources management requires indepth knowledge of the characteristics of watersheds and river systems, including artificial structures .
Given the limitations of monitoring data, and the need of the process-based understanding of watershed behaviour, the most sensible approach for large-scale water balance estimation should be based on hydrological simulation models (White et al., 2015). Such an approach has been widely reported to be of fundamental importance for the understanding of national and continental impacts of climate and land-use changes, and for sustainable management of water resources (Jung et al., 2012;Werth & Guntner, 2010). Regional and catchment-scale hydrological models of fine (study site-adjusted) resolution are designed for simulating streamflow characteristics using detailed local input data and applying specific calibration procedures. Numerous applications of processbased, (semi-)distributed hydrological models have been reported in the literature at the continental (Abbaspour et al., 2015;Schuol et al., 2008), national (Faramarzi et al., 2009) and large river basin scale (Karabulut et al., 2016).
In recent years, hydrological model Soil and Water Assessment Tool (SWAT; Arnold et al., 1998) was set-up and calibrated for the Vistula and Odra basins (VOB). It covered 89% of Poland and small parts of the neighbouring countries (Piniewski et al., 2017a). It provided estimates of natural streamflow and water balance components for the period 1954-2013. They have been made available for scientific use (Piniewski et al., 2017b). The model has been applied in multiple contexts, for example, for studying the effect of climate change on runoff , high and low flows (Piniewski et al., 2017c), soil moisture drought (Piniewski et al., 2020), flow conditions important for fish (O'Keeffe et al., 2018), wetland habitats (O'Keeffe et al., 2019) and anastomosing rivers (Marcinkowski & Mirosław-Świątek, 2020). Nonetheless, an important limitation of the existing VOB model was the lack of national coverage, and lack of representation of water management operations (i.e. water withdrawals, point sources, agricultural drainage systems and artificial reservoirs). Moreover, recent years have been characterized by increased presence of hydrological and agricultural droughts (e.g. 2015, 2018 and 2019) that were not covered by the previous model. Finally, new, higherresolution source data for various input maps such as soil and land cover maps have recently become available.
The objective of this study is to estimate streamflow and water balance components for the territory of Poland at a high spatial and temporal resolution using the SWAT model. Meeting this objective involved the following steps: (i) development of the model set-up, (ii) design of a spatial calibration approach and (iii) extensive calibration and validation of the model. The SWAT modelling performed in this study addresses the limitations of the previous VOB model by (i) extending the spatial coverage to the union of the entire Polish territory and the Vistula and Odra basins, (ii) extending temporal coverage from the period 1954-2013 to the years 1951-2020, (iii) using climate forcing at 2 km resolution (compared to the previous 5 km resolution), (iv) using new, higher-resolution land cover and soil maps and (v) including water management in the model set-up.

| STUDY AREA
The study area covers the union of the Polish territory and the Vistula and Odra basins (hereafter denoted as 'PL+'; set for 1951-2020 for the territory of Poland, 4TU.ResearchData. dataset, https:// doi.org/10.4121/16843183).

| Hydrological model
The SWAT model is a semi-distributed, process-based hydrological model simulating the movement of water with a daily time step at a catchment scale . In SWAT, river basins are divided into sub-basins (enclosing a single river reach), further partitioned into hydrologic response units (HRUs) constituting a combination of uniform land cover, soil and slope within each sub-basin. Water balance is computed separately for each HRU, and F I G U R E 1 Location of the study area then aggregated at the sub-basin level and routed downstream through the stream network to the outlet. SWAT incorporates a number of equations and sub-models to account for relevant processes determining water balance such as evapotranspiration, snowmelt, surface runoff generation, soil moisture, interception, infiltration, saturated and unsaturated water storage and routing components.

| Model input data sets
The basic geospatial data required to develop the model structure cover land elevation (topography), soil type, land cover and stream network. It should be emphasized that due to presence of transboundary rivers, for parts of the model domain located outside Poland, global data sets characterized by lower resolution had to be used. The specific data used for the PL+ basins are listed in Table 1. Geospatial data sets and maps of land cover and soils were pre-processed before implementing them into the model domain. High-resolution (1-m) DEM (for the territory of Poland) was aggregated to a resolution of 50 by 50 meters, and combined with resampled (also to 50 m resolution) Copernicus EU-DEM and SRTM4 DEM (for the territory outside Poland). The bilinear resampling method was used for aggregating 1-m DEM, and the spline interpolation method for resampling 90-m SRTM4 to 50 m resolution (Her et al., 2015). For the territory of Poland, the standard CORINE Land Cover 2018 (CLC) map was enhanced by including the Imperviousness Density (EU Copernicus database, 2018) layer in the urbanized areas, and refined within heterogeneous agricultural areas using more univocal Sentinel-2 Global Land Cover 2017 (S2GLC) agricultural classes. Because the CLC map does not cover the territory of Ukraine and Belarus, S2GLC map served as the primary land cover map for the part of PL+ model domain located in these countries. Moreover, for the eastern part of Belarus (approx. 4,500 km 2 , 1.3% of total study area), not covered by any of the data sets used in this study, basic land cover classes (urbanized areas, water, arable land, meadows and forests) were digitized from the ESRI orthophoto base map. Then, the actual crop distribution acquired from the agricultural census (2010, Statistics Poland) for 2477 communes in Poland and from ONESOIL database (https://map.oneso il.ai/2018) for the area outside Poland was applied within the arable land polygons. The methodology for the distributionpreserving random scatter of eight crops (spring barley, winter wheat, corn for grain and silage, potato, sugar beet, spring canola and vegetables) within the arable land class of the land cover map was the same as in Piniewski et al. (2015). Compared to the VOB model (Piniewski et al., 2017a), a major advancement in crop number (eight vs five) and level of detail (2477 communes vs coarser 314 districts) was achieved. The soil map within the Polish borders was created by pre-processing and generalizing Polish source maps (1:5,000, 1:25,000 and 1:100,000) to a reference scale of 1:100,000. The generalization process was performed in two steps: first, polygons with an area lower than 5 ha were deleted from the source maps; then, the resulting gaps were filled by adding their geometry to the neighbouring polygon with the longest shared edge. For the territory outside Poland, the generalized soil map was combined with the product characterized by lower resolution (1:1,000,000), that is, the Harmonized World Soil Database (HWSD).
Temporally and spatially variable weather data are the major driving force in SWAT. Daily minimum and maximum temperatures (°C) and precipitation totals (mm) were obtained from a gridded 2 km daily climate data set for the union of the Polish territory and the Vistula and Odra basins (G2DC-PL+, Piniewski et al., 2021). While originally the G2DC-PL+ meteorological data set covers the years 1951-2019, for the purpose of this study it was extended to 2020, following procedures described by Piniewski et al. (2021). Another difference compared to G2DC-PL+ data is that the latter used precipitation data corrected for gauge undercatch following the Richter method, whereas here we used the same gridded precipitation product but with no gauge undercatch correction.
Daily discharge time series with a temporal coverage of 2000-2019 for 88 gauging stations ( Figure 1) were obtained from the Institute of Meteorology and Water Management and used for the calibration and validation of the PL+ hydrological model.

| Hydrological model set-up
In the process of automatic watershed delineation, the model domain was divided into 4,381 sub-basins with an average area of 79.8 km 2 (10th percentile of 8.5 km 2 and 90th percentile of 164 km 2 ), and the total delineated area of 349,766 km 2 . The correctness of DEM-based delineated sub-basins was assessed at 88 flow gauging stations with known drainage areas. The mean absolute percentage error of the drainage areas upstream of the gauges was 1.9%. An overlay of land cover, soil and slope maps within SWAT sub-basins resulted in the designation of 47,725 HRUs with an average area of 7.33 km 2 (10th percentile of 2.8 km 2 and 90th percentile of 13.1 km 2 ). Source maps of land cover, soil and slope were added to the supplementary material (Figures S1, S2 and S3). The SWAT model allows the user to select different methods and equations for simulating particular water balance components, depending on the input data availability. Here, we provide a short description of methods/ equations selected in the SWAT model set-up for PL+. The modified USDA Soil Conservation Service (SCS) curve number method was selected for calculating surface runoff. We selected an option in which the retention parameter of the SCS CN method varies with plant evapotranspiration. For potential evapotranspiration (PET) modelling, the Hargreaves method was selected because it requires exclusively minimum and maximum temperature records. Channel routing was modelled using the Muskingum method.

| Calibration and validation procedures
Sensitivity analysis, calibration and validation were performed with the application of SWAT-Calibration Uncertainty Procedures (SWAT-CUP; Abbaspour, 2008), using the Sequential Uncertainty Fitting Procedure Version 2 (SUFI-2 algorithm; Abbaspour et al., 2004) with a parallel processing module. Kling-Gupta efficiency (KGE) was used as an objective function (Gupta et al., 2009). KGE is a function of the correlation term (linear regression coefficient between the measured and simulated variable), the variability term (the ratio between simulated and measured standard deviation), and the bias term (the ratio between simulated and measured mean).
Where r is the linear regression coefficient (correlation term), ∝ = s ∕ m is the variability term, = s ∕ m is the bias term, with s and m are standard deviations, and s and m are the means of simulated and measured variables respectively.
The percent bias (PBIAS) that measures the average tendency of the simulated data to be smaller or larger than their observed counterparts and coefficient of determination (R 2 ) were also tracked.
In this study, a two-stage multi-site calibration was performed based on daily discharge time series covering the period 2000-2009. The calibration was preceded by sensitivity analysis that provided a list of the most sensitive parameters. At the first stage of actual calibration, we used data from 72 flow gauges representing medium-sized (447 do 3,429 km 2 ), non-nested catchments uniformly distributed across Poland (Figure 1). The calibration catchments were clustered into seven groups, based on similarity in response (i.e. simulated flow) to parameter changes observed upon sensitivity analysis. For each group of calibration catchments, a unique set of best performing parameter values was derived upon calibration. In the next step, the parameter regionalization technique was used in which parameter values were transferred from seven donor catchment groups to target (ungauged) sub-basins ( Figure S4) based on geographical distance (He et al., 2011). Specifically, for each sub-basin that was not part of any of the calibration catchments, the closest group with its corresponding parameter set was assigned.
At the second stage of calibration, we used data from 16 flow gauges located at the largest Polish rivers: four on Vistula, three on Odra and the remaining nine on their main tributaries: Wieprz, Pilica, Bug, Narew, San, Dunajec, Warta and Noteć (Figure 1). In this step, we altered none of the parameter values prescribed across all HRUs and sub-basins in the first stage of calibration. We only adjusted two parameters related to the streamflow routing phase of the hydrological cycle: reservoir parameter NDTARGR (Number of days to reach target storage from current reservoir storage) and channel roughness parameter CH_N2. It was necessary due to the presence of several large reservoirs in the model domain, and specific flow conditions in large rivers that were not observed in the medium-sized calibration catchments. This step of calibration did not affect the objective function values in any of the 72 gauges used in the first step.
Model validation was performed for an independent period of 2010-2019 based on daily flow records for 88 (72 and 16 from the first and second stage respectively) gauging stations.

| Calibration and validation results
The SWAT model performance measured by KGE was satisfactory for almost all catchments (median KGE = 0.73) in the calibration period (Figures 2 and 3). Comparable results were obtained for corresponding catchments in the validation period (median KGE = 0.7, Figures 2 and  3). It is worth emphasizing that the KGE values obtained in 72 calibration catchments in the PL+ SWAT model exceed the values obtained in 80 calibration catchments of the VOB model (Piniewski et al., 2017a). While the difference between the models' performance for the calibration period is inconsiderable (median KGE = 0.73 vs 0.69), it is notably higher for the validation period (KGE = 0.7 vs 0.62). Model performance measured by R 2 was satisfactory for almost all catchments for the calibration (median R 2 = 0.67) and validation (median R 2 = 0.68) periods as well (Moriasi et al., 2015). Regarding PBIAS, according F I G U R E 2 KGE, R 2 and PBIAS values assigned to calibration catchments derived upon calibration (a) and validation (b). Flow gauges marked as black-frame circles refer to stage 2 of the calibration process on the main rivers (see Section 3.5) to Moriasi et al. (2007), the median value for calibration (4.1%) and validation (−5.5%) allows for rating the model performance as good in most of the catchments.
Global data sets, characterized by lower resolution, were used for the part of the study area outside Poland. It is, therefore, important to assess how that affected the model F I G U R E 3 Box plots showing the model performance expressed by KGE, R 2 and PBIAS values in 88 calibration catchments. Groups 1-7 (Stage 1) refer to 72 gauges used at the first stage of calibration, and Stage 2 refers to the remaining 16 gauges located on large rivers performance. It involved the comparison of the goodnessof-fit measures for two transboundary rivers (Bug and Odra) at two gauging stations located immediately downstream of the country border, with median values for all 16 flow gauges located at the largest Polish rivers (KGE = 0.80 and 0.77, R 2 = 0.74 and 0.76, PBIAS = 4.2%, and −2.5%, for calibration and validation periods respectively). For the Bug River (partly draining Ukraine and Belarus), for both the calibration (KGE = 0.76, R 2 = 0.70, PBIAS = 6.3%) and validation (KGE = 0.72, R 2 = 0.66, PBIAS = −13.2%) period, goodness-of-fit measures are only slightly worse, but still satisfactory according to Moriasi et al. (2015). For the Odra River (with part of the upper catchment located in the Czech Republic), model performance measures are also slightly worse (KGE = 0.71 and 0.69, R 2 = 0.72, and 0.64, PBIAS = 25.1%, and 8.3%, for calibration and validation periods respectively). The presented results suggest that the use of global data sets, characterized by lower resolution for part of the study area, only moderately affected the model's performance.

| Spatio-temporal changes in water balance
The calibrated and validated SWAT model of PL+ was used to run a simulation for the entire period of availability of G2DC-PL+ meteorological data . Here, we have extracted selected key water balance components (precipitation, evapotranspiration, water yield and discharge) expressed as decadal (10-years) mean values to analyse their spatio-temporal patterns (Figures 4 and 5). The remaining SWAT output variables are presented in the supplementary material ( Figure S5-S13). Precipitation  (Figure 5a) are the highest in the southern (mountainous) and northern (seaside) parts of Poland. The large central part of the country is generally characterized by lower precipitation and water yield values. Such a spatial pattern seems to be consistent in all decades, but significant differences in magnitude of precipitation and water yield are observed. During the 70 years of simulation, the highest mean precipitation total was observed in the years 2001-2010 (660 mm), and the lowest in the period 1951-1960 (575 mm). Regarding the water yield, the highest values were recorded for the years 1971-1980 (171 mm) and the lowest for the period 1951-1960 (116 mm). Actual evapotranspiration ( Figure 4b) seems to be less variable in time and space, ranging from 450 mm to 487 mm between the decades. The spatial pattern follows land cover and soils to a larger extent than climate. Mean decadal discharge (Figure 5b) was summarized from all basin outlets (Vistula, Odra and other coastal rivers) draining to the Baltic Sea. A significant variation in discharge was observed between the decades. The highest values were recorded for the period 1971-1980 (1,825 m 3 /s), and the lowest for the years [1951][1952][1953][1954][1955][1956][1957][1958][1959][1960] (1,232 m 3 /s).
For better assessment of evapotranspiration, the SWAT results were compared with 0.1-degree resolution ERA5-Land global reanalysis data set (Muñoz-Sabater et al., 2021). Figure S14 suggests good spatial agreement between the data sets across Poland with locally higher ET values in the southern, eastern and northern parts of the country, and lower in the central part. Moreover, Figure S15 presents a similarly increasing temporal trend and inter-annual variability in annual ET totals for both data sets for the time period 1951-2020 (R 2 = 0.41). It is worth emphasizing that R 2 is diverse through decades, and increases over the last 20 years (R 2 = 0.54). There is, however, a general disagreement in the magnitude of ET F I G U R E 5 Mean decadal water yield (a) and discharge (b) at the sub-basin level in the period 1951-2020. Values above the maps denote spatially-averaged values for the entire model domain between the data sets. The ERA5-Land data set presents higher ET values (23% on average) compared to the PL+ SWAT model. Although ERA5-Land is one of the best reanalysis data sets available to date, we think it provides overestimated ET over Poland.

| DATA AVAILABILITY
The PL+ SWAT model hydrology data set was derived for years 1951-2020 at the sub-basin and reach level, with a daily time step. At the sub-basin level, average daily water balance components, namely precipitation total, snowfall, snowmelt, evapotranspiration (actual and potential), surface runoff, percolation, groundwater recharge, groundwater flow, tile flow and water yield (all components in mm/day), and soil water content (mm), were included. At the reach level, average daily streamflow values (m 3 /s) were added ( Table 2). The data set called PL-SWAT-51_20 is stored in an online, publicly available research data archive 4TU.Research Data, https://doi. org/10.4121/16843183 (Marcinkowski et al., 2021). The file format is netCDF files. Additional shapefiles containing river outlets, reaches and sub-basins are also available.

| CONCLUSIONS
This paper presents a detailed description of the development the hydrological SWAT model for the union of the Polish territory and the Vistula and Odra River basins (PL+). A gridded 2 km daily climate (temperature and precipitation) data set (G2DC-PL+) was used to reconstruct the meteorological conditions. The SWAT model was extensively calibrated and validated following a twostage multi-site calibration method for 88 gauging stations across Poland. The simulations show satisfactory performance with a median Kling-Gupta efficiency of 0.73 and 0.7 in the calibration and validation period respectively. The calibrated and validated SWAT model for the PL+ area was used to run a simulation for the period 1951-2020, reconstructing water balance and streamflow data at a daily time step for 4,381 sub-basins. The water balance components are characterized by remarkable spatio-temporal variation during the simulated 70 years. The PL+ SWAT model hydrology data set with a daily time step  includes the following water balance components: precipitation total, snowfall, snowmelt, evapotranspiration (actual and potential), surface runoff, soil water, percolation, groundwater recharge, groundwater flow, tile flow, water yield and streamflow records. The data set is stored in an online, publicly available research data archive 4TU.Research Data (Marcinkowski et al., 2021). The results of the presented paper can be further particularly used as the basis for water resources management and planning, but also for other sectors such as agriculture, energy production and freshwater ecology. The model can be used as a tool for assessment of climate change impacts on hydrology, as well as other scenarios such as land use or water management change.