*Article* **Warming Air Temperature Impacts Snowfall Patterns and Increases Cold-Season Baseflow in the Liwiec River Basin (Poland) of the Central European Lowland**

**Urszula Somorowska**

Department of Hydrology, Faculty of Geography and Regional Studies, University of Warsaw, Krakowskie Przedmie´scie 30, 00-927 Warsaw, Poland; usomorow@uw.edu.pl

**Abstract:** The rapidly changing climate affects vulnerable water resources, which makes it important to evaluate multi-year trends in hydroclimatic characteristics. In this study, the changes in cold-season temperature (November–April) were analyzed in the period of 1951–2021 to reveal their impacts on precipitation and streamflow components in the Liwiec River basin (Poland). The temperature threshold approach was applied to reconstruct the snowfall/rainfall patterns. The Wittenberg filter method was applied to the hydrograph separation. The Mann–Kendall test and Sen's slope were applied to estimate the significance and magnitude of the trends. An assessment of the similarity between trends in temperature and hydroclimatic variables was conducted using the Spearman rank-order correlation. The shift-type changes in river regime were assessed via the Kruskal–Wallis test. The results revealed that temporal changes in both snowfall, rainfall, and baseflow metrics were significantly associated with increasing temperature. Over 71 years, the temperature rose by ~2.70 ◦C, the snowfall-to-precipitation ratio decreased by ~16%, the baseflow increased with a depth of ~17 mm, and the baseflow index rose by ~18%. The river regime shifted from the snow-dominated to the snow-affected type. Overall, this study provides evidence of a gradual temperature increase over the last seven decades that is affecting the precipitation phase and streamflow component partitioning in the middle-latitude region.

**Keywords:** cold season; 1951–2021; trends; snowfall-to-precipitation ratio; baseflow index; river regime shift; lowland river basin; middle-latitude region

#### **1. Introduction**

The warming climate, which is considered one of the most important factors that affects stream-flow regimes in many regions of the world, has environmental and socioeconomic implications, particularly with respect to the vulnerability of water resources [1–3]. With the increasing air temperature, altered precipitation patterns change the water quantities, thereby contributing to runoff components as quick flow and baseflow [4,5]. In addition to the amount of precipitation, the precipitation phase (snowfall or rainfall) plays a critical role in runoff generation [6,7]. Snowfall, if persistent, stores cold-season precipitation into the spring months and keeps the hydrological system dormant, while rainfall infiltrates soils, recharges aquifers, and feeds streams and rivers. The warming air temperature might reduce snowfall and amplify snow melt, thereby resulting in a decline in water storage in snowpacks, earlier snow, and soil thawing, and, as a consequence, causing a shift in the hydrological regime [8–10]. Since the warming climate modifies the ratio of snow to precipitation (S/P ratio) in many parts of the world [11–13], the detailed characterization of temperature and precipitation changes is a high priority in ongoing research [14–16]. However, significant uncertainties remain regarding the current and future trends in the hydrologic implications of warming due to the high sensitivity of hydrological processes to climate variability and change [17].

**Citation:** Somorowska, U. Warming Air Temperature Impacts Snowfall Patterns and Increases Cold-Season Baseflow in the Liwiec River Basin (Poland) of the Central European Lowland. *Resources* **2023**, *12*, 18. https://doi.org/10.3390/ resources12020018

Academic Editor: Diego Copetti

Received: 29 October 2022 Revised: 12 January 2023 Accepted: 13 January 2023 Published: 17 January 2023

**Copyright:** © 2023 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Divergent trends in snow characteristics occur in different regions of the world. When considering snowfall, climate warming might increase humidity, which enhances extreme snowfall, whereas the rising temperature reduces the likelihood of snowfall [18]. For example, the authors of [19] found that the frequency of daily snowfall events tended to decrease across much of the Northern Hemisphere except at the highest latitudes such as in northern Canada, northern Siberia, and Greenland. Divergent snowfall patterns were also uncovered in [18]; it was found that while higher-latitude regions experience increasing extreme snowfall percentiles, decreasing extreme snowfall percentiles are characteristic of lower-latitude regions in Western Europe. Observational evidence from the pan-European in situ data was provided in [20]; it was revealed that over the period of 1951–2017, the mean snow depth decreased more than the extreme snow depth, and widespread decreases in the maximum and mean snow depth were found over Europe except in the coldest climates. It is worth noting that interannual variability in the extent of snow is high, and new extremes in maximum snow metrics over Eurasia have occurred in recent years [16,21]. Generally, the increases in winter temperatures have resulted in a decrease in the S/P ratio and an increase in winter snowmelt for most of the Northern Hemisphere; this was particularly significant in the middle-latitude regions [11,22].

Divergent assessments also concern the hydrologic implications of rising temperature reported in global and regional studies. In the coldest river basins, the response to warming is manifested by an increase in the spring streamflow peak, whereas for the transitional basins, the spring runoff decreases [1]. This is because transitional river basins face large increases in winter streamflow. While many studies have generally focused on mountainous and arctic regions [23–26], relatively fewer studies have reported on snow hydrology changes across low-relief topography regions [27,28]. These are regions in which the river runoff is also sensitive to the effects of the changing patterns of snow accumulation and melt. In [1], regions with snowmelt-dominated runoff were selected using the ratio of accumulated annual snowfall divided by annual runoff (S/QT ratio or snowfall-torunoff ratio), and the criterion of S/QT > 0.5 over the global land regions was applied. According to this criterion, vast areas of the Central European Lowland located within the borders of Poland are not classified as snowmelt-dominated areas. However, this may be due to the short observation period, which covered the years 1980–1999, and the coarse spatial resolution of the gridded data (0.5◦ × 0.5◦ latitude/longitude) used in the analysis. Therefore, it seems reasonable to evaluate the S/QT ratio using ground-based data with a higher spatial resolution for a longer period of time. Such an analysis might unmask multi-year changes in the S/QT ratio over time. In the east and northeast parts of Poland, the hydrological regime of lowland rivers is influenced by relatively persistent seasonal snowpack in winter, with the highest streamflow occurring in spring when the snow mass melts and feeds the rivers. The strong influence of the snowmelt-dominated river regime is manifested by the relatively high Pardé coefficient in the spring months, which reaches or exceeds 180%. Therefore, it is unclear whether this region also did not have a snowmelt-dominated character before the 1980s. While some previous studies investigated selected aspects of decreases in snow-cover depth driven by the rising winter temperature in Poland in the years 1966/1967–2019/20 [29], no studies have considered the trends in wintertime runoff components and their dependence on rising temperatures and changing snowfall patterns. To fill this gap, this paper focused on quantifying the changes in snowfall, rainfall, baseflow, and quick flow, as well as their association with temperature during the cold season. The multi-year trends were examined in an exemplary typical mesoscale river basin of the Central European Lowland. This basin is situated in the middle eastern region of Poland, where the most intensive temperature increase was recorded at selected weather stations followed by a decrease in snow cover depth [29]. Decreasing tendencies were found in the snow metrics; however, the variability in the snow characteristics within both the winter season and the multi-year period was high [30].

The key scientific questions that needed to be answered in this study were as follows:


It was hypothesized that the snowfall has remarkably decreased in the last seven decades, thereby leading to a decrease in S/P ratio. It was also hypothesized that the increasing trends in cold-season temperature are responsible for the changes in streamflow, with a higher component of baseflow being recharged by the increased fraction of rainfall to total precipitation.

#### **2. Materials and Methods**

#### *2.1. Study Area*

This study concerned the Liwiec River basin, a left tributary of the Bug River, which is situated in the Mazovian Lowland in central eastern Poland (Figure 1). The study area is located between 52◦00 N–52◦30 N and 21◦30 E–22◦0 E and belongs to the Central Poland Lowland, which is a part of the Central European Lowland [31]. The stream gauge is situated at the Łochów cross-section and closes the river basin at 2471 km2. The elevation ranges from 83 m a.s.l. near the outlet to 227 m a.s.l. in the southern part of the basin. The Liwiec River is about 142 km in length. The basin is influenced by a snowy, humid climate (Dfb) with a warm summer (see the updated Köppen–Geiger classification [32]). Monthly precipitation and air temperature are evenly distributed across the catchment and differentiate seasonally in four meteorological seasons: spring (March–May), summer (June–August), autumn (September–November), and winter (December–February). The annual precipitation is about 550 mm, of which 250–300 mm falls in the winter half-year (November–April). The average annual air temperature is 7.5 ◦C. The lowest temperatures are recorded in January; the average monthly temperature for this month is −4 ◦C. The highest values are recorded in July (an average monthly value of 19 ◦C). The studied area is characterized by a greater amplitude of air temperature compared to central and western Poland due to moderate continental impacts. The river regime is considered to be nival, with the highest streamflow rates usually occurring in March; the lowest streamflow rates usually occur in July and August and sometimes continue throughout September–October. Low flows might also occur in winter when snow cover blocks the groundwater recharge and the groundwater resource gradually depletes. When snow melts, a pronounced peak flow occurs. The aquifers are found in the Quaternary formations. They are directly recharged by precipitation, and the Liwiec River and its tributaries constitute a zone of natural drainage. The depth of the first groundwater table varies depending on the lithology of the surface formations, but in most of the area it does not exceed 5 m. In the highlands, aquifer sands are often found under less permeable loams. In the river valleys and in the flat areas in the north and south of the basin area, the unconfined aquifers form the subsurface. The land use is dominated by agricultural land (63%), which includes arable land (48%) and meadows (15%). The forest type is the second-most dominant (33%), while the rest of the territory is covered by artificial surfaces (3%) and other categories (1%).

**Figure 1.** (**a**) Geographic location of the study area. (**b**) Elevation map of the Liwiec River basin according to EU-DEM v1.1 (acquired from https://land.copernicus.eu; accessed on 13 October 2022), and distribution of the E-OBS gridded temperature and precipitation points (acquired from https://www.ecad.eu; accessed on 19 June 2022).

#### *2.2. Datasets*

For the analysis of climate conditions, the air temperature and precipitation dataset was employed for 1950–2021 from version 25.0e of the station-based E-OBS gridded dataset (https://www.ecad.eu; accessed on 19 June 2022) available from the European Climate Assessment and Dataset Project [33]. It comprised daily precipitation (P) and air temperature (T) values acquired from a regular latitude/longitude grid of 0.1◦ × 0.1◦ (Figure 1b). For the Liwiec River basin, the data subset was extracted as a 3D-gridded dataset of 6 × 12 × 25,933 dimensions, which corresponded to 6 latitude grid cells, 12 longitude grid cells, and 25,933 daily solutions. Using the daily gridded values, the basin scale, monthly air temperature, and precipitation estimates were calculated, and then the average values for the months November–April were calculated to represent the cold-season conditions. This approach followed the convention of the water year, which was designated in Poland as a 12-month period (November–October) that consists of the "winter half-year" (cold season) and "summer half-year" (warm season) [34,35]. It is worth noting that the cold season is considered to be the 6-month period starting on November 1 of a particular year and lasting until the end of April in the next calendar year. For example, the values representing the cold season of the year 1951 were calculated using the monthly values starting in November 1950 and ending in April 1951. Using these monthly values, the cold-season precipitation and air-temperature time series were calculated and tested for the presence of a trend as explained in Section 2.3.2.

For the analysis of snowfall patterns, the air-temperature threshold approach was applied to partition the precipitation into rain and snow fractions [36]. Since the E-OBS data did not differentiate between rain and snow, it was assumed that all precipitation fell as snow below a threshold of 1.20 ◦C and as rain above this threshold. This threshold is near the upper bound of the transition temperature across the Liwiec River basin as revealed in [36]. The threshold temperature across the basin is evenly distributed in a narrow range of 1.19 to 1.22 ◦C. The threshold data covering the territory of Poland were prepared at a resolution of 0.5◦ × 0.625◦ (latitude × longitude) based on 89 meteorological stations distributed across the country. The spatial resolution of the threshold temperature data was coarser than the E-OBS dataset; thus, resampling to a common resolution of 0.1◦ × 0.1◦ was required to keep the information from the finest layers. Using the 3Dgridded air-temperature dataset of 6 × 12 × 25,933 dimensions, a binary mask was prepared; grids with a temperature above the threshold were set to zero while the remaining grids received a value of 1. This binary representation of snowfall occurrence was used to differentiate between snow and rain in daily precipitation data stored in the 3D matrix of 6 × 12 × 25,933 dimensions. Moreover, the daily snowfall and rainfall values were accumulated into monthly values and averaged for the entire river basin. For the 6-month cold season, the air temperature (T), snowfall (S), and rainfall (R) time series consisted of 71 values covering the period of 1951–2021.

To analyze the possible impacts of the warming climate on the river regime, the daily streamflow data for gauging station no. 152210120 (Liwiec River, Łochów cross-section) were acquired for the water years of 1951–2021 from https://danepubliczne.imgw.pl/data (accessed on 30 August 2022). The data were prepared and verified by the Institute of Meteorology and Water Management—National Research Institute in Warsaw, Poland (IMGW-PIB). The procedure used to separate the hydrograph into baseflow and quick flow is described in Section 2.3.1. For the cold season, the streamflow (QT), baseflow (QB), and quick flow (QQ) time series were tested for the presence of trends as explained in Section 2.3.2.

#### *2.3. Methods*

#### 2.3.1. Hydrograph Separation into Quick Flow and Baseflow Components

The daily streamflow (QT) time series were analyzed for the water years of 1951–2021, and the two components of baseflow (QB) and quick flow (QQ) were separated from the total flow (QT). The HYDRORECESSION toolbox was used for the streamflow recession analysis [37]. Aksoy and Wittenberg's method [38] was applied to the extraction of recession segments, in which negative dQT/dt values were considered to represent the recession of the hydrograph composed of the baseflow. Here, the minimum duration of the recession segment was set to 10 days, and the filter criterion (removed days) was assumed to be 5 days. The Wittenberg filter method was applied to the daily streamflow for baseflow separation, which assumed that the baseflow recession segment satisfied the nonlinear relationship of S = aQb. The optimal model parameters were determined via linear regression (least squares). The cold-season and monthly values of QB and QQ were extracted and examined. Then, the baseflow index was calculated as BFI = QB/QT, and the quick flow index was determined as QFI = QQ/QT.

#### 2.3.2. Trend Detection

A non-parametric Mann–Kendall (MK) test [39,40] was applied to detect trends in the precipitation, air temperature, and streamflow data time series. Moreover, the components of precipitation (snowfall and rainfall) and streamflow (baseflow and quick flow) were also tested for the presence of temporal trends. The MK test was used to test the null hypothesis of no trend (H0) against the alternative hypothesis (H1) that there was an increasing or decreasing monotonic trend. The trends were tested at a significance level of alpha = 0.05. The magnitude (slope) of an existing trend (as change per year) was calculated using the directional coefficient expressed by the Theil–Sen estimator [41,42]. A positive slope value suggested an increasing trend, and a negative slope value indicated a decreasing trend. If change was not statistically significant but showed an inclination, it was called a tendency. The Climate Data Toolbox (CDT) for MATLAB [43] was used to calculate the MK standardized test statistic (Z) and the *p*-value.

#### 2.3.3. Spearman Rank-Order Correlation for Similarity Assessment between Trends

As a measure of the strength of the link between trends in the monthly time series of temperature and other variables (including snowfall, rainfall, baseflow, and quick flow), a non-parametric Spearman's rank-order correlation coefficient (RS) was calculated and evaluated at a significance level of alpha = 0.05. This assessed the relationship between

two variables without making any assumptions about the frequency distribution of the variables. The Spearman rank-order correlation was equal to the Pearson correlation between the rank values of the two variables and ranged between −1 and 1. For each variable, the monthly trend rates were ranked and the strength of the link between trends in air temperature and other variables was assessed. The tested hypothesis was that with the rising temperature, snowfall would decrease, rainfall would increase, and significantly more baseflow would occur. Statistically significant results were given a *p*-value < 0.05. The pairs with positive-correlation coefficients tended to increase or decrease together, while negative-correlation coefficients indicated a relationship between two variables in which an increase in one variable was associated with a decrease in the other. The MATLAB Statistics and Machine Learning Toolbox (Release 2021) was used for all statistical analyses.

#### 2.3.4. Shift-Type Changes in the River Regime

In this study, the ratio of the accumulated annual snowfall to the annual total runoff (S/QT ratio) was used to examine the role of snowmelt in the seasonal streamflow patterns. This metric was introduced in [1] to determine whether or not runoff was snowmeltdominated using the criterion of S/QT > 0.5. The shift-type changes were analyzed, and the time series of S/QT was partitioned into two subseries. The separation into two subperiods was achieved by minimizing the sum of the residual (squared) error of each subset from its local mean and finally returning the index, which in this case was the year in which the change occurred. The change point was identified using the MATLAB function "findchangepts". A complete 71-element S/QT time series was divided into subperiods of differing lengths that consisted of 20 and 51 records covering the years of 1951–1970 and 1971–2021, respectively. Finally, a Kruskal–Wallis test was applied to test for statistically significant differences between the subseries in the two selected subperiods [44,45]. This was a non-parametric test that compared the mean ranks (i.e., medians). For this test, the null hypothesis was that the subseries medians would be equal (versus the alternative of a difference between them).

#### **3. Results**

#### *3.1. Cold-Season Trends in Hydroclimatic Variables over the Period of 1951–2021*

Figure 2 shows the course of the cold season's hydrometeorological variables over the years of 1951-2021, including the precipitation (P11-04), air temperature (T11-04), snowfall (S11-04), snowfall-to-precipitation ratio (RSP11-04 = S11-04/P11-04), rainfall (R11-04), rainfall-toprecipitation ratio (RRP11-04 = R11-04/P11-04), streamflow (QT11-04), baseflow (QB11-04), quick flow (QQ11-04), and baseflow index (BFI11-04). Statistically significant changes occurred in the time series of air temperature (increasing trend; Figure 2b), snowfall-to-precipitation ratio (decreasing trend; Figure 2d), rainfall-to-precipitation ratio (increasing trend; Figure 2f), baseflow (increasing trend; Figure 2h), quick flow (decreasing trend; Figure 2i), and baseflow index (increasing trend; Figure 2j). The T11-04 showed a trend rate of ~0.38 ◦C/decade, while a trend rate of ~−2.27%/decade was detected in RSP11-04. Consequently, RRP11-04 showed an increase of ~2.27%/decade. Thus, the gradual temperature increase was accompanied by a decrease in the snow-to-precipitation ratio and an increase in the proportion of the liquid phase in the precipitation (RRP). The warmest cold season occurred in 2020 with a T11-04 of 4.4 ◦C, while the coldest occurred in 1963 with a T11-04 of −3.2 ◦C (Figure 2b). Generally, the lowest values of T11-04 were recorded in the first half of the analyzed period, while the highest were recorded in the last two decades. In response to the slightly increasing tendency in P11-04 and the significant increase in T11-04, QT11-04 did show a slight decreasing tendency of ~0.048 m3s−1/decade (Figure 2g), which was equivalent to 0.06 mm/decade, while QB11-04 increased with a trend rate of ~0.39 m3s−1/decade (Figure 2h), which was equivalent to 2.46 mm/decade. Thus, the changes in QB11-04 reached ~17.44 mm over 71 years.

**Figure 2.** Changes in cold-season (**a**) precipitation (P11-04), (**b**) air temperature (T11-04), (**c**) snowfall (S11-04), (**d**) snowfall-to-precipitation ratio (RSP11-04 = S11-04/P11-04), (**e**) rainfall (R11-04), (**f**) rainfall-toprecipitation ratio (RRP11-04 = R11-04/P11-04), (**g**) streamflow (QT11-04), (**h**) baseflow (QB11-04), (**i**) quick flow (QQ11-04), and (**j**) baseflow index (BFI11-04). The MK test statistic is denoted as Z. The presence of a trend was determined at a significance level of alpha = 0.05; in cases with a *p*-value > 0.05, the changes were not statistically significant. Sen's slope is expressed in the variable unit per year.

#### *3.2. Trends in Monthly Hydroclimatic Variables over the Period of 1951–2021*

Temperature is a prime factor that determines the occurrence of precipitation as rain or snow [46]. Therefore, this study primarily examined the air-temperature trends over the multi-year period of 1951−2021 as an area-weighted average across the entire river basin. Figure 3 shows the course of monthly air temperature (T) in the cold season that lasted from November to April. The multi-year rate of change in the monthly T was equal

to 0.24, 0.27, 0.32, 0.49, 0.53, and 0.31 ◦C/decade for consecutive months of the cold-season period, respectively.

**Figure 3.** Changes in monthly air temperature in (**a**) November (T11), (**b**) December (T12), (**c**) January (T01), (**d**) February (T02), (**e**) March (T03), and (**f**) April (T04). The MK test statistic is denoted as Z. The straight red line in each figure represents the linear trend (Sen's slope) of the monthly average temperature (expressed in ◦C/y).

For the months of November, February, March, and April, these changes were statistically significant; the highest rates occurred in February and March. In December and January, the change also showed an increase.

Examining changes in the monthly snowfall and snowfall-to-precipitation ratio yielded substantially consistent results (Figure 4). With an increasing temperature, in all six months of the cold season, a decreasing tendency was seen regarding the snowfall amount, which was statistically significant in February. In January, February, and March, the snowfallto-precipitation ratio decreased ~4%/decade in each month (Figure 4f,h,j); therefore, in 71 years, it decreased by ~28%. When considering the change in the six-month cold season, S11-04 decreased at a rate of ~2.5mm/decade (Figure 2c) for a total decrease of ~18 mm in 71 years, and RSP11-04 decreased by 2.3%/decade (Figure 2d), which accounted for 16% in 71 years. Consequently, the monthly rainfall increased in five of the six months of the cold season (from December to April, with the highest increase in March) by ~2.56 mm/decade (Figure 5i), which represented a rate of ~18 mm in 71 years. The highest increase in RRP occurred in the months of December, January, February, and March within a range of ~3.2–4.4%/decade (Figure 5d,f,h,j), which represented a range of 23–31% in 71 years.

Changes in the monthly baseflow and BFI are shown in Figure 6, while the changes in the quick flow and QFI are presented in Figure 7. A statistically significant increase in QB occurred in January, February, and March; while in November, December, and April, growing tendencies were registered. In all six months, BFI showed an increase that was

statistically significant in March and April (Figure 6j,l). Inverse changes were noted for QFI, which showed statistically significant decreasing changes in March and April (Figure 7j,l).

**Figure 4.** Changes in monthly snowfall in (**a**) November (S11), (**c**) December (S12), (**e**) January (S01), (**g**) February (S02), (**i**) March (S03), and (**k**) April (S04); and changes in the monthly snowfallto-precipitation ratio in (**b**) November (S11/P11), (**d**) December (S12/P12), (**f**) January (S01/P01), (**h**) February (S02/P02), (**j**) March (S03/P03), and (**l**) April (S04/P04)). The MK test statistic is denoted as Z. The straight red line in each figure represents the linear trend (Sen's slope) of monthly snowfall (expressed in mm/y) or the monthly snowfall-to-precipitation ratio (expressed in %/y).

**Figure 5.** Changes in monthly rainfall in (**a**) November (R11), (**c**) December (R12), (**e**) January (R01), (**g**) February (R02), (**i**) March (R03), and (**k**) April (R04), and changes in monthly rainfall-toprecipitation ratio in (**b**) November (R11/P11), (**d**) December (R12/P12), (**f**) January (R01/P01), (**h**) February (R02/P02), (**j**) March (R03/P03), and (**l**) April (R04/P04). The MK test statistic is denoted as Z. The straight red line in each figure represents the linear trend (Sen's slope) of monthly rainfall (expressed in mm/y) or monthly rainfall-to-precipitation ratio (expressed in %/y).

**Figure 6.** Changes in monthly baseflow in (**a**) November (QB11), (**c**) December (QB12), (**e**) January (QB01), (**g**) February (QB02), (**i**) March (QB03), and (**k**) April (QB04), and changes in monthly baseflow index in (**b**) November (BFI11), (**d**) December (BFI12), (**f**) January (BFI01), (**h**) February (BFI02), (**j**) March (BFI03), and (**l**) April (BFI04). The MK test statistic is denoted as Z. The straight red line in each figure represents the linear trend (Sen's slope) of the monthly baseflow (expressed in m3s−1/y) or monthly baseflow index (expressed in %/y).

**Figure 7.** Changes in monthly quick flow in (**a**) November (QQ11), (**c**) December (QQ12), (**e**) January (QQ01), (**g**) February (QQ02), (**i**) March (QQ03), and (**k**) April (QQ04); and changes in quick flow index in (**b**) November (QFI11), (**d**) December (QFI12), (**f**) January (QFI01), (**h**) February (QFI02), (**j**) March (QFI03), and (**l**) April (QFI04). The MK test statistic is denoted as Z. The straight red line in each figure represents the linear trend (Sen's slope) of the monthly quick flow (expressed in m3s−1/y) or the monthly quick flow index (expressed in %/y).

#### *3.3. Similarity between Trends over the Period of 1951–2021*

The Spearman rank-order correlation coefficients (RS) expressed the strength of a link between trends in the monthly temperature time series and other variables, including snowfall, rainfall, baseflow, and quick flow (Table 1). With the rising temperature, the tendency toward increases in the rainfall, rainfall-to-precipitation ratio, baseflow, and baseflow index was confirmed by the positive values of RS. Negative values of RS manifested that rising temperature where the snowfall, snowfall-to-precipitation ratio, quick flow and quick flow index had a tendency to decrease. Strong statistically significant results were revealed for the snowfall-to-precipitation ratio, rainfall-to-precipitation ratio, and baseflow.

**Table 1.** Spearman rank-order correlations (Rs) between multi-year trends (Sen's slopes) in the monthly temperature (T) and hydroclimatic variables of snow (S), snow-to-precipitation ratio (S/P), rainfall (R), rainfall-to-precipitation ratio (R/P), baseflow (QB), baseflow index (BFI), quick flow (QQ), and quick flow index (QFI).


<sup>1</sup> Correlation was statistically significant at a significance level of alpha = 0.05.

#### *3.4. Shift in River Regime*

A shift-type change in the snow-to-runoff ratio (S/QT) occurred in 1971 as detected by the change point analysis (Figure 8). The Kruskal–Wallis test confirmed the significance of the difference between the median values of S/QT in the two subperiods. The mean values decreased from 0.61 in the subperiod of 1951–1970 to 0.44 in the subperiod of 1971–2021. In the first snowy period (1951–1970), the highest values of S/QT reached 0.95, while the lowest reached 0.35. In 1971–2021, the highest S/QT did not exceed 0.82 and very often dropped below 0.35; the lowest value of 0.09 occurred in 2020. When considering the mean S/QT in the two subperiods and the threshold criterion of S/QT = 0.5 applied in [1], it was concluded that the river regime shifted from the snow-dominated to the snow-affected type with a mixed recharge by both the snow and rainfall precipitation phases.

**Figure 8.** Snow-to-runoff ratio (S/QT) in two subperiods (1951–1970 and 1971–2021) separated by a change point in 1971.

It is worth noting that the term "snow-dominated" refers to the S/QT ratio. As shown in Figure 2d, the snowfall (S11-04) rarely exceeded rainfall (R11-04). In the subperiod of 1951–1970, the S11-04/P11-04 ratio exceeded 50% in 1958 (63%), 1968 (53%), 1969 (56%), and 1970 (53%); the multi-year mean value was 43%. In the subperiod of 1971–2021, the mean of S11-04/P11-04 dropped to 33%. The S11-04/P11-04 had maximum values in 1979 (62%), 1996 (76%), 2005 (53%), and 2013 (66%) and exhibited strong interannual variations. On average, rainfall (R11-04) exceeds snowfall (S11-04) (Figure 2f), and the S11-04/P11-04 ratio gradually decreased (Figure 2d).

#### **4. Discussion**

This study demonstrated that the air temperature remarkably increased in the last seven decades at both seasonal and monthly time scales over the cold season of 1951–2021. The results of the analyses supported the study's hypotheses. Warming winters have gradually reduced the snowfall amount and snowfall fraction of total precipitation. However, this did not exclude the occurrence of extremely snowy winters, an example of which was the snowfall in 2012/2013, which was the highest in the entire 70 years in this river basin (Figure 2c). However, the 2019/2020 season was marked by extremely low snowfall and the lowest share of snow in the precipitation, which was caused by the exceptionally high temperatures that occurred throughout the cold season. This was the mildest winter on record across Europe, particularly in the north and east [47]. Overall, in the Liwiec River basin, the cold-season temperature rose by ~2.70 ◦C over 71 years, and the snowfall-to-precipitation ratio (S/P ratio) decreased by ~16% over 71 years. The warming winter temperatures across Poland were previously confirmed; it was found that at the majority of weather stations in Poland, the snow-cover depth significantly decreased in recent decades [29,30]. It is worth mentioning that atmospheric thaws alternately occur with cool and frosty periods and are characteristic features of the climate of Poland. These are caused by the variable weather conditions in winter seasons [48]. In the analyzed river basin, in the years of 1960/1961–2009/2010, the mean number of days with atmospheric thaw in December–February was in the range of 10–12 days [48]. However, the extreme thaw-start and thaw-end dates could differ by more than three months due to the high interannual variability. Thus, with increasing temperatures, it might be expected that the frequency of atmospheric thaws would increase, thereby accelerating snowmelt and making the snow cover less persistent. An increased activity of the hydrological system is expected to manifest with the amplified infiltration process occurring on large flat surfaces, which recharges the groundwater from which an increased baseflow is generated.

The follow-up hypothesis was also confirmed; it was found that the relationship between the monthly trends in air temperature and baseflow was strong and statistically significant. Hence, the increasing cold-season temperature trend contributed to the changes in streamflow marked by an increase in baseflow and baseflow index (BFI). The baseflow component of streamflow increased to a depth of ~17 mm over 71 years, and the baseflow index rose by ~18% in 71 years. The BFI increased from 0.51 in the subperiod of 1951–1970 to 0.66 in the subperiod of 1971–2021. The opposite results were found regarding the quick flow metrics (QFI); the time series of the quick flow and quick flow index showed decreasing trends. This may have been due to more frequent snow- and soil-thawing periods and the lack of a sudden amount of melting snow quickly entering the river in the form of overland flow. The results of this study seemed to be consistent with findings in [49], in which it was proved that runoff along the Vistula River in the winter season has become more uniform and shows decreasing maxima and increasing minima of daily flows and a stable mean runoff volume. The study in [4] also made significant contributions to the knowledge of global trends by showing that changes in both the baseflow and BFI were significantly region-dependent. The rivers in the eastern part of Poland were not examined, but the western part (covering the Oder River basin) showed a decreasing baseflow and BFI for the winter season (December–January–February) when evaluated for the period of 1970–2016. Such contrasting results were, as proved in [4], region-dependent; the eastern

part of Poland where the Liwiec River basin is situated has a much cooler and humid climate than the western part.

#### **5. Conclusions**

The main conclusions of this work can be summarized as follows:


Overall, this study provided evidence of the link between trends in temperature and other hydroclimatic characteristics. However, the obtained results only applied to the analyzed river basin and require further verification for a broader range of geographic and climatic conditions in the middle-latitude region. Thus, future studies should consider other river basins. In Poland, the existing historical climate and hydrometric data enable further analysis of this area. Comparative studies that examine other river basins could explain the regional similarities or differences related to snow-hydrology and river-regime transitions. Such research could provide an answer regarding which river basins are still characterized by a snow regime and which have already undergone this transformation.

The changes in the hydrological regime analyzed in this paper only concerned the amount of snowfall related to the total annual runoff as expressed by the snowfall/runoff ratio. When considering the river regime changes, shifts in the time of the maximum winter/spring flows and their impact on the occurrence of low summer flows are also important. Such an analysis would be particularly valuable from the perspective of seasonal water resource availability.

The obtained results suggested that diminishing snowfall is expected to alter the groundwater recharge and streamflow dynamics if warming trends continue. Thus, the analysis presented here could be followed by a consideration of climate change scenarios and their impacts on snowfall and streamflow seasonality. Although this analysis was limited to showing similarities in the studied hydroclimatic variable trends, it seemed to provide valuable insights into the drivers and causes of changes in the river regime in the middle-latitude lowland region.

**Funding:** This research was funded by the Faculty of Geography and Regional Studies, University of Warsaw, Poland, under grant nos. SWIB66/2022 and SWIB85/2022.

**Data Availability Statement:** All of the data used during this study are available at the locations cited in the Acknowledgements section.

**Acknowledgments:** The author thanks the editors and reviewers for their insightful comments and suggestions. The author acknowledges the E-OBS dataset (version 23.1e) from the EU-FP6 project ENSEMBLES (http://ensembles-eu.metoffice.com) and the data providers in the ECA&D project (http://www.ecad.eu); version 25.0e of the dataset was downloaded from the Copernicus Climate Change Service (https://surfobs.climate.copernicus.eu/dataaccess/access\_eobs.php) (accessed on 19 June 2022). The author also acknowledges the gridded Northern Hemisphere 50% rain–snow Ts threshold product (a formatted version of the observational dataset) that was downloaded from the DRYAD: spatial variations in the rain–snow temperature threshold across the Northern Hemisphere (https://doi.org/10.5061/dryad.c9h35) (accessed on 31 January 2022). The author also acknowledges the streamflow dataset that was prepared and verified by the Institute of Meteorology and Water Management—National Research Institute in Warsaw, Poland (IMGW-PIB) and downloaded from https://danepubliczne.imgw.pl/data (accessed on 30 August 2022).

**Conflicts of Interest:** The author declares no conflict of interest. The funder had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Forecasting Monthly River Flows in Ukraine under Different Climatic Conditions**

**Renata Graf 1,\* and Viktor Vyshnevskyi <sup>2</sup>**


**Abstract:** River-flow forecasts are important for the management and planning of water resources and their rational use. The present study, based on direct multistep-ahead forecasting with multiple time series specific to the XGBoost algorithm, estimates the long-term changes and forecast monthly flows of selected rivers in Ukraine. In a new, applied approach, a single multioutput model was proposed that forecasts over both short- and long-term horizons using grouped or hierarchical data series. Three forecast stages were considered: using train and test subsets, using a model with train-test data, and training with all data. The historical period included the measurements of the monthly flows, precipitation, and air temperature in the period 1961–2020. The forecast horizons of 12, 60, and 120 months into the future were selected for this dataset, i.e., December 2021, December 2025, and December 2030. The research was conducted for diverse hydrological systems: the Prut, a mountain river; the Styr, an upland river; and the Sula, a lowland river in relation to the variability and forecasts of precipitation and air temperature. The results of the analyses showed a varying degree of sensitivity among rivers to changes in precipitation and air temperature and different projections for future time horizons of 12, 60, and 120 months. For all studied rivers, variable dynamics of flow was observed in the years 1961–2020, yet with a clearly marked decrease in monthly flows during in the final, 2010–2020 decade. The last decade of low flows on the Prut and Styr rivers was preceded by their noticeable increase in the earlier decade (2000–2010). In the case of the Sula River, a continuous decrease in monthly flows has been observed since the end of the 1990s, with a global minimum in the decade 2010–2020. Two patterns were obtained in the forecasts: a decrease in flow for the rivers Prut (6%) and the Styr (12–14%), accompanied by a decrease in precipitation and an increase in air temperature until 2030, and for the Sula River, an increase in flow (16–23%), with a slight increase in precipitation and an increase in air temperature. The predicted changes in the flows of the Prut, the Styr, and the Sula rivers correspond to forecasts in other regions of Ukraine and Europe. The performance of the models over a variety of available datasets over time was assessed and hyperparameters, which minimize the forecast error over the relevant forecast horizons, were selected. The obtained RMSE parameter values indicate high variability in hydrological and meteorological data in the catchment areas and not very good fit of retrospective data regardless of the selected horizon length. The advantages of this model, which was used in the work for forecasting monthly river flows in Ukraine, include modelling multiple time series simultaneously with a single model, the simplicity of the modelling, potentially more-robust results because of pooling data across time series, and solving the "cold start" problem when few data points were available for a given time series. The model, because of its universality, can be used in forecasting hydrological and meteorological parameters in other catchments, irrespective of their geographic location.

**Keywords:** river flow; XGBoost algorithm; trends; multistep-ahead forecasting; climate variability; Ukraine

**Citation:** Graf, R.; Vyshnevskyi, V. Forecasting Monthly River Flows in Ukraine under Different Climatic Conditions. *Resources* **2022**, *11*, 111. https://doi.org/10.3390/ resources11120111

Academic Editors: Diego Copetti and Demetrio Antonio Zema

Received: 11 October 2022 Accepted: 26 November 2022 Published: 30 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Some of the most popular scientific issues of modern hydrology are climate variability and change and its impact on the water regime, which is of primary importance for states with generally scarce water resources. Modelling and forecasting the time series of river flows are essential elements in the assessment of the hydrological regime and the management of water resources [1–3]. The results of analyses of multiyear flow measurement series enable an assessment of the reaction rivers to supply factors or their limitation. The decomposition of the time series of flows allows us to capture the change trends, seasonality, and a random factor, which is most often a factor disturbing the forecasting of time series [4]. The results of the forecasts, among others, are important in preventing floods and droughts and reducing the effects of their occurrence [5].

Typically, a statistical and deterministic approach is used in hydrological forecasting. Forecasting models fall into two main groups: the physics-based numerical models and the data-driven prediction models [6]. Physics-based models perform mathematical modelling to simulate dynamic processes, such as floods. Data-driven models are also widely used to model and forecast the flow or rainfall–runoff relationship [7]. It has been shown that data-driven approaches can also achieve a comparable performance to that of physicsbased models in predicting flows and extreme hydrological phenomena. According to Tu et al. [8] data-driven models, especially machine-learning models, are becoming alternative approaches to hydrological and hydraulic models. Predictive machine-learning algorithms are used in many studies in the field of applied and data-driven hydrology [9–11]. Shen et al. [12] emphasize that machine-learning models provide improvements to nonlinear modelling over other data-driven techniques. The effectiveness of artificial neural network (ANN) models in hydrological forecasting often exceeds the effectiveness of traditional conceptual or mathematical models based on the modelling of complex hydrological processes [13].

The modelling of temporal flow series is usually performed for both short-term forecasting [14] and long-term forecasting and often for high or low flows [15]. The concept of a short-term real-time daily unsteady flow forecast refers to the relationship between the instantaneous states (water gauge compounds), also without taking into account rainfall data [14,16,17]. Among the long-term forecasting methods, both qualitative and quantitative methods are used [18]. Quantitative methods (referred to as traditional methods) often include statistical methods, correlation, and regression analysis [19–23]. There are many methods available for predicting river flow, including process-based models, which, however, require a lot of data [24]. Models based on computational intelligence and machine-learning algorithms are gaining popularity [6]. Most often, good flow forecast results are obtained by comparing several models [16,25]. In the development of flow forecasts and the assessment of their quality, the conventional linear autoregressive relationship (AR), ANN models (e.g., three-layer feedforward neural network), recursive neural networks (RNN)), and a number of hybrid models are used. Accurate river-flow forecasts are obtained with the RNN model, which often also has the greatest ability to generalize results, showing similar forecast quality in independent tests.

Among the neural networks, the model for the short term works well for river-flow prediction [26–29]. The ANN models have also been used to model rainfall-runoff, riverflow, and flood forecasting, among others, by Imrie et al. [26] for selected catchments in the UK; by Kim and Barros [30] for selected rivers in Pennsylvania (the US); and by Toth et al. [25] for rivers in the Apennines (Italy). Comparisons of prediction and prognostic models of river flows are often made on the basis of autoregressive techniques, neural networks and adaptive neuro-fuzzy inference systems (ANFIS) [31,32]. In the short-term to long-term forecasting of river flows, heuristic optimization algorithms hybridized with ANFIS work well. All developed hybrid algorithms significantly outperformed the traditional ANFIS model performance for all prediction horizons, which is a major advantage compared with the classical black box machine-learning models [33,34]. The performance of neural network models in river-flow forecasting has been compared with linear regression

models [35] and stochastic models [36–38]. Forecasts have been made for one river by using different models, or assessments have been made by using different methods (alternative methods) for contrasting catchments. The forecasts are based on flows observed in the cross sections of the river system with a daily delay and without taking into account any rainfall data. Rainfall-runoff simulation modelling is also used, taking into account RNNs [39] as well as deep learning with a long short-term memory (LSTM) network approach [13,40]. The LSTMs are a type of recursive neural networks capable of learning sequence dependencies in sequence prediction problems. Toth et al. [25] have conducted a comparison of short-term rainfall prediction models for real-time Floyd forecasting.

Multimodel data fusion is a tool commonly used in river-flow forecasting. Abrahart and See [16] have conducted data fusion for two different catchments using arithmetic averaging, the probabilistic method (in which the best method is used to generate the current forecast), the last time step model, two different neural network operations, and two different soft computing methodologies. Each site demonstrated several options and potential benefits of using data fusion tools to produce better estimates of hydrological forecasts.

With the improvement of data-driven modelling in hydrological applications, teamlearning methods (e.g., gradient boosted decision tree, or GBDT) are widely used [41]. Ensemble machine-learning methods are employed to solve problems related to simulation and prediction in hydrology, including resampling methods (bagging, boosting, and dagging), model averaging, and stacking [42]. Boosting methods (e.g., boosting, adaboost, and extreme gradient boosting) are becoming more and more effective at modelling and forecasting runoff, flooding, and drought and at predicting ice phenomena in rivers [43]. They are used in hydrological modelling to achieve better performance by combining many weaker models into a stronger model [44]. After developing a new streamflow forecasting model based on the modular model, Ni et al. [45] confirmed that models based on extreme gradient gain (XGBoost) give satisfactory forecasts. The results of these studies for the Yangtze River monthly flows indicate that XGBoost is applicable to river-flow forecasting and generally performs better than the support vector machine (SVM). Extreme gradient boosting (XGBoost) has been used to predict flood stages in a data-driven flood alert system (FAS) for a flood-prone watershed in Houston (Texas) [41]. It has been shown that an XGBoost-based FAS can operate continuously to automatically detect flooding, without needing activate an external start-up trigger, as is usually required in conventional eventbased warning systems. Flood warning systems built on predictive models constitute a proactive method of flood risk assessment and management [46].

Forecasting changes in the water-balance and hydrological regime is particularly important and valuable for countries struggling with a shortage of water resources and the effects of climate variability and change on such resources. The impact of the human factor is largely uncontrolled, which is noted in Ukraine, where hostilities have been taking place since February 2022. Water resources in Ukraine are becoming increasingly scarce, with the water infrastructure being increasingly exposed to war. In Ukraine, the subject of river-flow forecasting in relation to multiyear measurement series or high (flood) or low flows is undertaken on a different scale. In many studies conducted on Ukrainian rivers, forecasting refers to maximum flows and spring floods, in which the authors rely on the traditional approach, i.e., the use of quantitative methods [47]. In contrast, Khrystiuk and Gorbachova [15] have conducted long-term forecasting on extraordinary spring floods by means of the commensurability method on the Dnipro River near Kyiv, Ukraine. The method used was the Weng Wen-Bo information method, which is a qualitative forecasting method that makes it possible to identify the periods and specific years in which subsequent extraordinary spring river floods may occur. The impact of climate on water resources has also been studied for selected river catchments in Ukraine, e.g., for the Western Bug River in Western Ukraine [48], this by using the water-balance Turk model and the regional model (REMO) [49], climate change scenarios [50] and ecohydrological modelling [51]. Didovets et al. [52] and Vyshnevskyi and Donich [53] have studied one climate change

impact on regional water runoff in the Carpathian region, and Loboda and Kozlov [54] rated the water resources of Ukrainian rivers by using average statistical models for the period 2021–2050.

The aim of the research was to determine long-term changes and forecast the flows of selected rivers located in different regions of Ukraine during the period 1961–2020, taking into account changes in precipitation and air temperature. The research covered the following rivers: the Prut (a mountain river), the Styr (an upland river), and the Sula (a lowland river), taking into account the forecast horizons of 12, 60, and 120 months. Monthly river-flow forecasting was carried out in three stages, using input data sets specific to the XGBoost algorithm: (1) forecasting by train and test subsets, (2) forecasting using a model with train-test data, and (3) forecasting by training with all data. The algorithm used allowed direct multistep-ahead forecasting with multiple (grouped or hierarchical) time series, which is important in the case of taking limited measurements of daily river flows. The research results can be used for the management and rational planning of water resources in the studied regions of Ukraine.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Ukraine is located in the south-eastern part of Europe. The larger part of the country (about 95%) is covered by plains, while the rest is mountainous. For researching long-term changes and predicting river flow in Ukraine, three medium-size river catchments with an area of 500–15,000 km<sup>2</sup> were selected: the Prut River in the mountains (in the Carpathians) and the Styr and Sula rivers on the plains (Figure 1). Three main requirements governed the selection of these rivers: the location within river basins of at least two meteorological stations, a long-term observation period, and a slight impact of human activity on river runoff [55]. We may add that the last demand is especially important in the case of Ukraine, because an important influencing factor, namely irrigation, is present in the south of the country.

In the upper course the Prut River up to the Yaremche hydrological station (with a catchment area of 597 km2) is a typical mountain river. Moreover, the source of the river is located near the Hoverla mountain, which is the highest in Ukraine (with an elevation of 2061 m) (Figure 1). A major part of the catchment area is covered by forests (primarily spruce). The climate of the catchment area of this river essentially depends on elevation. At the Pozhezhevska meteorological station, at a height of 1451 m a.s.l., the mean annual air temperature is equal to 3.5 ◦C, the mean air temperature in January is minus 5.8 ◦C, and the mean air temperature in August (the warmest month) is 13.3 ◦C. On the other hand, the air temperature at the Yaremche meteorological station (at an elevation of 531 m, near the Yaremche hydrological station) is significantly higher: in the period 1991–2020, the mean annual air temperature equalled 7.9 ◦C. The coldest month is January (minus 2.2 ◦C), and the warmest one is July (17.8 ◦C). Under cold conditions, the snow cover in the upper part of the mountains may persist until early June. The mean annual precipitation in the catchment area is above 1000 mm; at the Pozhezhevska meteorological station, it is 1548 mm; and at the Yaremche station, it is 1013 mm. There are waterfalls on the river, including the Probiy waterfall, several kilometres upstream from Yaremche hydrological station. At the same time, there is no pond or any other water reservoir upstream the station. The water regime of the Prut River depends mainly on the snow cover depth and heavy rains. As a rule, the largest river flow is observed in April, during the snow-melting season. Meanwhile, the largest water discharges are caused by heavy rains in the summer period. The main flow characteristics of the Prut River (Yaremche water gauge) are shown in Table 1.

**Figure 1.** The location of the catchment areas of the studied rivers (**a**), the hydrological (marked in blue) and meteorological (marked in red) stations (**b**): 1, the Prut–Yaremche, 2, the Styr–Lutsk, and 3, the Sula–Lubny. Note that the name of the catchment area includes the name of the river and the water gauge stations (catchment closing cross section).

**Table 1.** The main characteristics of the discharge conditions of the studied rives, data modified, according to [55].


The Styr River is located in the north-western part of Ukraine (Figure 1). It is a tributary of the Pripyat River, and the Pripyat River is the largest tributary of the Dnipro River. The upper, southern part of the catchment area upstream features the town of Lutsk, on Volyn-Podilsk Upland, while the lower part is on the Polesian Lowland. The climate of the catchment area is moderate with cool winters and warm summers. During the period 1991–2020, the mean annual air temperature at the Lutsk meteorological station was equal to 8.5 ◦C; the mean air temperature in January was minus 2.9 ◦C, whereas in July, it was 19.7 ◦C. The air temperature at other stations, located upstream, was a little higher (by 0.1–0.2 ◦C). The precipitation in the catchment area upstream from the town of

Lutsk amounts to an average of about 600–700 mm. The water regime of the Styr River at the Lutsk hydrological station (its catchment area is 7200 km2) is rather stable. The largest discharges observed in April and May are about twice as large as those observed at low water in August and September. In the upper course of the river, there are some reservoirs and many ponds. The reservoirs are rather small, and they do not significantly impact the water regime. The same can be said about ponds, many of which are older than 100 years. The main flow characteristics of the Styr River (Lutsk water gauge) are shown in Table 1.

The Sula River catchment is on the Dnipro Lowland with plain relief (north-east Ukraine), in the forest steppe zone, which actually has little forest (Figure 1). The climate of the catchment area is moderate, with cool winters and warm, sometimes hot, summers. During 1991–2020, the mean annual air temperature at the Lubny meteorological station was equal to 8.6 ◦C, the mean air temperature in January was minus 4.1 ◦C, and the mean air temperature in July was 21.3 ◦C. The air temperature at another station, Romny, located upstream, was lower by approximately 0.5 ◦C. The precipitation in the catchment area upstream from the town of Lubny amounts to an average of about 600 mm. The water regime of the river, comparable to the Styr River, is much more unstable. The mean discharges at the Lubny hydrological station (with a catchment area of 14,200 km2), observed in April and May, are approximately 10 times larger than those observed at low water in September. The flow regulation of the Sula River is small. The main flow characteristics of the Sula River (Lubny) are shown in Table 1.

#### *2.2. Dataset Characteristics*

The dataset consisted of the average monthly river flows, precipitation, and air temperatures from the period 1961–2020, acquired for three studied catchment areas. The monthly river-flow data were analysed for the Prut–Yaremche, the Styr–Lutsk, and the Sula–Lubny hydrological stations, and for the meteorological stations (Figure 1b), data on the average monthly air temperature and precipitation were analysed. For the Prut River basin, data were used from the Pozhezhevska meteorological station, in the upper part of the river basin, and from the Yaremche meteorological station, near the water gauge. For the Styr River catchment, data from three meteorological stations were used: Brodu (in the south), Dubno (in the east), and Lutsk( near the hydrological station). For the Sula River basin, data from the Romny (in the north-east part of the river catchment) and Lubny (in the south) meteorological stations were used. The time series consisted of the average monthly data for all points of observation.

#### *2.3. Descriptive Statistics*

Descriptive statistics of distribution of river-flow variables showed that the threshold values of skewness and kurtosis were exceeded (Table 2); therefore, for the purposes of further analysis, nonparametric measures of central tendency, i.e., Median (Quartiles 1–3), were used. As for the distribution of precipitation and air temperature, the skewness values did not exceed 2.0 and the kurtosis −7.0 (Table 3); hence, the nonparametric measures were used for further analysis (mean and standard deviation).

**Table 2.** Descriptive statistics of river flow (m3/s) grouped for the period 1961–2020 (N = 720).


Explanations: N—sample size; M—mean; SD—standard deviation; Mdn—median; IQR—interquartile range (IQR = Q3 − Q1); Min—minimum value; Max—maximum value; Skew.—skewness; Kurt.—kurtosis.


**Table 3.** Global descriptive statistics of precipitation and air temperature grouped by river catchment area for the period 1961–2020 (N = 720).

Explanations as in Table 1. The time series consisted of the average monthly data for all points of observation.

Compared with the Prut and Styr rivers, the Sula River was characterized by a high standard deviation of monthly flows from the average value, which means that the monthly flows showed high variability over the analysed period (Table 2). The large variance in the flow series was also confirmed by the fact that the highest IQR (22.36) to median ratio (second quartile Q2) was also determined for the Sula–Lubny and the lowest for the Styr–Lutsk.

Table 3 shows that the catchment area of the Prut River had the highest mean precipitation (about 100 mm per month); for the other two rivers' catchments, the precipitation was approximately 50 mm. The degree of variability was similar for all rivers and represented approximately 60% of the average value. The mean air temperature at the Yaremche station during 1961–2020 was about 7.4 ◦C, and for the other two rivers, the temperature was in the range 8.4–8.5 ◦C (Table 3). The degree of variability was similar for all rivers and was about 90–100% of the average value.

#### *2.4. Forecast Algorithm*

River-flow and meteorological variables forecasting was carried out in several stages, as shown in flowchart (Figure 2). The algorithm used allowed direct multistep-ahead forecasting with multiple (grouped or hierarchical) time series by using the built-in methods of the "R-forecast ML" package [56] inspired by Bergmeir et al. [57]. The applied direct forecasting approach included the following steps: (1) building a single multioutput model that simultaneously forecasts over short- and long-term forecast horizons, (2) assessing the generalization performance of the model over a variety of available data sets over time, and (3) selecting the hyperparameters that minimize forecast error over the relevant forecast horizons and re-train.

The description of the algorithm is broken down into individual steps in Figure 2.

The first step was to create model training and forecasting datasets with lagged, grouped, dynamic, and static features. For this purpose, the create\_lagged\_df() function of the "R-forecast ML" package was used. The basic idea behind creating custom-feature lags was to improve model accuracy by removing noisy or redundant features in highdimensional training data. For this purpose, three horizons were defined with forecasts for 1:12, 1:60, and 120 months into the future. A numerical vector indicating the lags in the dataset rows for the creation of the lagged features (lookback) was defined at levels of decreasing range (features from 1 to 36 months, from 54 to 78 months, and from 120 to 126 months in the past) with monthly frequency. The variables month and year were defined as dynamic variables with no static variables in the dataset. Based on multiple time series, the variable "river\_id" identified the groups (hierarchies) that were used as model features but were not lagged. The same procedure was performed for the test data set in the model and forecast stage.

**Figure 2.** General flowchart of a multistep-ahead forecasting algorithm.

The next step involved creating time-contiguous validation datasets for model evaluation. Given that measurements are taken monthly, each validation window consisted

of 24 months, skipping the next 12 months. As a result, the dataset contained 10 validation windows. Custom validation subsets with skipped river-flow data are shown in Figure S1 (Supplementary Materials). A similar approach was used for precipitation and air temperature.

Training the user-defined model across forecast horizons and validation datasets was carried out in the next step. The user-defined model consisted of all data transformations, hyperparameter tuning, and inner loop cross-validation procedures with XGBoostspecific [58] input datasets created within this wrapper function.

The next step involved computing the forecast error across forecast horizons and validation datasets (Figure S1 in Supplementary Materials). Estimating forecast error metrics on the test datasets was based on the calculation of the following: mean absolute error (MAE), mean absolute percentage error (MAPE), median absolute percentage error (MDAPE), symmetrical mean absolute percentage error (SMAPE), root mean squared error (RMSE), and normalized root mean squared error by standard deviation (NRMSE), which have been discussed in detail in the works of, i.a., Bahrami-Pichaghchi and Aghelpour [59], and Aghelpour and Norooz-Valashedi [60].

The final step involved combining the multiple direct-horizon forecast models to produce a single h-step-ahead forecast. Prediction intervals were estimated by calculating the α/2 and 1-α/2 percentiles of the simulated data for each prediction horizon by using bootstrap residuals and recursive multivariate prediction. The one-step-ahead prediction error was defined as:

$$e\_t = y\_t - \hat{y}\_{t|t-1} + \varepsilon \tag{1}$$

where *yt* represents the actual prediction, *y*ˆ*t*|*t*−<sup>1</sup> represents the prediction conditioned by previous prediction,  represents the estimation error predictions, which were simulated by including the residuals (sampling from the collection of errors from past prediction iterations). The expected variance was represented by the difference in predictions.

The approach was to combine forecasts across models in such a way that short-term models would produce short-term forecasts and long-term models would produce longterm forecasts. This implies that for 120 months, the ahead forecast included (1) the 12-step-ahead model forecasts from the next month through 12 months, (2) the 60-stepahead model forecasts from months 13 through 60, and (3) the 120-step-ahead model forecasts from months 61 through 120.

Analyses were performed with the statistical language R (version 4.1.1) [61].

#### XGBoost Algorithm

For modelling and forecasting river flows and meteorological variables, an XGBoost was chosen as it outperforms other machine-learning (ML) algorithms, such as artificial neural networks (ANN) or support vector regression (SVR), and statistical models in the case of a multilevel predictive average [62]. The process is called gradient boosting because it makes use of an algorithm to minimize losses when adding new models (Figure 3). Training is an iterative process: new trees are added that predict the residuals or errors of the previous trees, and they are subsequently combined with the previous trees to make the final prediction. XGBoost minimizes a regularized (L1 and L2) objective function that combines a convex loss function (based on the difference between predicted and target outputs) and a penalty term for model complexity.

**Figure 3.** Schematic principle of the XGBoost algorithm for the training process (the figure was prepared on the basis of XGBoost documentation [63]).

In Figure 3, *αi* and *ri* are regularization parameters and residuals computed in the *i* th tree separately. Here, *hi* is a function that is trained to predict results, and *ri* uses *X* for the *i* th tree.

To compute *αi*, the computed residuals used *ri* and estimated the following:

$$\text{arg}\min\_{\infty} = \sum\_{i=1}^{m} L\left(Y\_{i\nu} \ F\_{i-1}(X\_i) + ah\_i \left(X\_{i\nu} \ r\_i\right)\right) \tag{2}$$

where *L* (*Y*, *F* (*X*)) is a differentiable loss function. *L*(*θ*) = ∑*<sup>i</sup>* (*yi* − *y*ˆ*i*) 2 , where *L* represents the training loss function, *θ* represents model parameters, *y* represents the actual value, and *y*ˆ represents the predicted value.

Using a sampling procedure, 80% of the sample was qualified for the training set and 20% for the test set. The training of the model was performed by using eXtreme Gradient Boosting Training with the following parameters of tree booster: maximum depth of a tree = 3; number of threads = 2; maximum number of boosting iterations = 20; control the learning rate *-* = 0.3; minimum loss reduction required to make a further partition on a leaf node of the tree γ = 0.5; minimum sum of instance weight (hessian) needed in a child = 5; and evaluation metric *rmse*. For task parameters, the regression with squared loss was used as an objective function with early stopping rounds = 5. The regression with squared loss is referred to as:

$$L(\theta) = \sum\_{i} (y\_i - \mathcal{y}\_i)^2 \tag{3}$$

where *L* represents the training loss function, *θ* represents model parameters, *y* represents the actual value, *y*ˆ represents the predicted value.

#### **3. Results**

*3.1. River-Flow Dynamics in Relation to Changes in Precipitation and Air Temperature over the Decades 1961–2020*

For a graphical representation of the flow distribution during the study period, see Figure 4. From the data provided, one can follow the dynamics of flow over the decades. In the case of the Prut River at the Yaremche station, variable dynamics was evident over the study period. The smallest values of flow were recorded in the period 1960–1970, the largest in the period 2000–2010. In the last decade, a significant decrease in flow was observed at the level of the local minimum for 50 years (median = 8.0 m3/s). A similar situation was observed in the case of the Styr River (Lutsk water gauge), where the maximum flow values recorded over 2000–2010 changed into the minimum flows during the entire observation period recorded over 2010–2020. In the last decade, the median was 23.7 m3/s (with the median of 26.9 m3/s for the period 1961–2020), which is 7.1 m3/s lower than that in the 2000–2010 decade. Similarly, in the case of the Sula River (Lubny water gauge), the period of increase in the flow between 1960 and 1990 was followed by a continuous decline, with a global minimum recorded in the last decade.

The trend dynamics of river-flow amounts over the decades is shown in Figure 5. In the case of the Prut–Yaremche, an alternating trend was observed, with periods of growth (1960–1980, 1990–2010) alternating with periods of decline (1980–1990, 2010–2020). The last decade was characterized by the sharpest decline during the study period, when the flow reached the global minimum. A similar situation was observed in the case of the Styr–Lutsk, where the flow decreased by more than 25% in 2010–2020. In the case of the Sula–Lubny, the period of increase in flow (1960–1990) was followed by a period of decrease at the end of the study period, when the global flow minimum was reached in 2010–2020.

**Figure 4.** Flow time series distribution grouped by rivers (1961–2020). Note that figure descriptions refer to the nonparametric measures of central tendency—median (Quartile 1–Quartile 3).

**Figure 5.** Flow trend values grouped by rivers (1961–2020). Note that figure descriptions refer to nonparametric measures of central tendency—median (Quartile 1–Quartile 3).

The seasonal component is the same for all rivers, the only difference being its variability. The lowest variability (0.3) was observed in the case of the Styr River at Lutsk, whereas

the highest in the case of the Sula–Lubny (0.9). The dynamics of the seasonal component of river flows in the period 1961–2020 is shown in Figure S2 (in Supplementary Materials). The level of the random component was similar for all rivers. Variability in the case of rivers Prut–Yaremche and Sula–Lubny was also similar, but it was lower for the Styr River at Lutsk. In all cases, a slow decline in the average component and variability was observed throughout the study period, with local or global minima in the last decade, except the Sula–Lubny, where a wave of variation was observed. The dynamics of the random component of river flow is shown in Figure S3. In the case of the Prut River catchment, the lowest precipitation was recorded between 1980 and 1990, and it was followed by an increase in precipitation with a global maximum observed in the last decade 2010–2020 (Figure S4). In the case of the Styr River catchment, a slight variation in the dynamics of precipitation was observed: the period of increasing precipitation alternated with the period of decreasing precipitation. In the Sula River catchment, following a clear increase in precipitation in the 1970s and a corresponding stable level during the next 40 years, a decrease in precipitation below 50 mm was recorded in 2010–2020, which corresponds to the local minimum for 50 years.

The trend dynamics of precipitation amounts over the decades is shown in Figure S5 (in Supplementary Materials). The precipitation data indicate that in the case of the Prut River catchment, the declining trend period in the 1960s–1990s was replaced by an increasing trend with the global maximum in the last decade of the study period. In the case of the Styr River at Lutsk, there was a slight variation in the trend over the period with increases/decreases during one or two decades with the global average at the end of the study period. In the Sula River catchment, maximum values of the precipitation trend were observed in the period 1970–2000, which was followed by a period of a downward trend with a local minimum in the last decade of the study period. The seasonal component is the same for the Prut River and the Styr River catchments. With the same average seasonality, the Sula River catchment was characterized by the lowest variability of the three rivers. The dynamics of the seasonal component of precipitation is demonstrated in Figure S6. The mean of the random component was the same for all the rivers' catchments, the degree of variability increasing slightly in the second half of the study period, reaching local or global maxima. The lowest variability was observed for the Prut River and the highest for the Sula River catchment. The dynamics of the random component of precipitation is demonstrated in Figure S7.

For all three river catchment areas, the lowest average air temperatures were recorded in the first decade of the study period. The study period was characterized by an almostcontinuous increase in average air temperatures. The value of the increase ranged from 1.6 ◦C to 2.1 ◦C with global maxima in the last decade (Figure S8 in Supplementary Materials). The trend dynamics of air temperature over the decades is presented in Figure S9. The data show a continuous upward trend for each river's catchment from the beginning of the observation to the end of the study period. Global trend maxima were reached in the last decade. The degree of variability for all river catchments was characterized by a period of increase, with the maximum in the period 1980–1990, followed by a decrease in variability to local or global minima. The seasonal component of air temperature was the same for the Prut River and Styr River catchments. With the same average seasonality, the Prut River catchment was characterized by the lowest variability (Figure S10). The level of the random component was similar for all the rivers' catchments, the variability was also similar (Figure S11). In all cases, a slow decline in the average component and variability was observed throughout the study period, with global minima in the last decade.

#### *3.2. Forecasting Monthly River Flows and Meteorological Variables*

River-flow, precipitation, and air temperature forecasting was performed in three stages: using train (80%) and test (20%) subsets, using a model with train-test data, and training with all data.

#### 3.2.1. River-Flow Forecasting

The fitting of the data within the selected validation subgroups is presented in Figure S12 (in Supplementary Materials), where the dashed lines represent the actual data not fitted by the model. Measures of data fit in terms of prediction errors can be found in Tables S1 and S2 (in Supplementary Materials). The data show that the fitting errors are quite large. This is best seen in the normalized RMSE (NRMSE), which was calculated as the ratio of the RMSE to the global SD of flow for each river, separately. Global NRMSE values of 0.78, 0.82, and 0.61 tell us that the model does not fit the data well enough to trust such predictive results (>0.5), but it is better than the usual data randomization based on global M and SD (<1). The absence of the overfitting and underfitting (the RMSE value is similar at each fitting step for the training and test sets) of the data is evidence for judiciously selecting model parameters. A significant RMSE value at each fitting step indicates high variability in the data that cannot be reasonably well fitted by retrospective data, regardless of the selected horizon length.

Given the above limitations, the prediction of river flow for the next approach was performed by using the training set and the test set for individual rivers in periods of 12, 60, and 120 months. The aggregated prediction plots are shown in Figure 6.

**Figure 6.** 12-, 60-, and 120-step-ahead forecasts of river flows (purple line: the Prut–Yaremche, turquoise line: the Styr–Lutsk, yellow line: the Sula–Lubny).

The predictions of river flow within each horizon are shown in Table 4.


**Table 4.** Forecast by river-flow model with train-test data grouped by horizon.

Explanations: N—sample size; Mdn—median; IQR—interquartile range.

An analogous prediction was made, taking into account all the data in the form of a test set. Figure 7 shows combined multiple direct-horizon forecast models for producing a 120-month step-ahead forecast grouped by "river id".

The predictions for the entire set are essentially consistent with the predictions in Table 4. Based on the results given in Figure 7, it can be concluded that the flow level of the Prut–Yaremche will remain the same or will decrease by about 6% in the period 2021–2030 (compared with the flow level in the period 2010–2020). In the case of the Styr River at Lutsk, the flow level will decrease by about 12–14%. In contrast, the flow of the Sula–Lubny in 2021–2030 is expected to be 1.2–1.7 m3/s higher than in 2010–2020, which corresponds to a projected increase of 16–23%.

**Figure 7.** Combined plot with partial actual data (2015–2020) and predicted data (2021–2030) for the river-flow model trained with all data. Explanations: The left side of the graph shows a portion of the actual data (2015–2020) and the red vertical bar shows the boundaries of the actual data and the beginning of the forecast data. The data in purple represents the forecast data. The purple circles reflect the values with a monthly step, and the filled background represents the prediction intervals.

#### 3.2.2. Precipitation Forecasting

In the case of precipitation forecasting, the fitting of the data within the selected validation subgroups is shown in Figure S13 (in Supplementary Materials), where the dashed lines represent the actual data not fitted by the model.

The data in Tables S3 and S4 (in Supplementary Materials) show that the fitting errors of prediction are quite large. This is best seen in the normalized RMSE (NRMSE), which was calculated as the ratio of the RMSE to the global SD of precipitation for each river catchment, separately. Global NRMSE values of 0.79, 0.86, and 0.96 indicate that the model does not fit the data well enough to trust such predictive results (>0.5) but are better than the usual data randomization based on global M and SD (<1). The absence of the overfitting and underfitting of the data (the RMSE value is similar at each fitting step for the training and test sets) is evidence for judiciously selecting model parameters.

Given the above limitations, as in the case of river-flow forecasting, precipitation forecasting using the training set and the test set approach were performed for individual rivers in periods of 12, 60, and 120 months. The aggregated forecast plots are shown in Figure S14. The predictions of precipitation values within each horizon are shown in Table 5.


**Table 5.** Forecast by precipitation model with train-test data grouped by horizon.

Explanations: N—sample size, M—mean, SD—standard deviation.

An analogous prediction (as for river-flow forecast) was made, taking into account all the data in the form of a test set. Precipitation predictions for the 120-month period for each river catchment are shown in Figure 8.

**Figure 8.** Combined plot with partial actual data (2015–2020) and predicted data (2021–2030) for the precipitation model trained with all data. Explanations are the same as those for Figure 7.

Based on the results provided in Figure 8, it can be concluded that in the period 2021–2030, the amount of precipitation in the Prut River catchment will decrease by about 7–10% (compared with the precipitation level in 2010–2020). In the case of the Styr River catchment, the amount of precipitation will remain at the level of the last decade or slightly decrease (by 1–2 mm). For the Sula River basin, the forecast shows an increase in average precipitation (by 1–3 mm) compared with the last decade.

#### 3.2.3. Air Temperature Forecasting

In air temperature forecasting, the fitting of the data within the selected validation subgroups is shown in Figure S15 (in Supplementary Materials), where the dashed lines represent the actual data, almost all of which were fitted by the model.

The data in Tables S5 and S6 (in Supplementary Materials) show that the fitting errors of prediction were rather moderate. This is best seen in the normalized RMSE (NRMSE), which was calculated (similar to river flows and precipitation) as the ratio of the RMSE to the global SD of temperature for each river, separately. Global NRMSE values of 0.28–0.30 tell us that the model provides a 3.3–3.5 times better prediction of the data than in the case of the usual data randomization based on global M, SD. A low RMSE value at each fitting step indicates that the data can be reasonably well fitted by retrospective data, regardless of the horizon length chosen.

Given the above limitations, air temperature forecasting for the approach with the training set and the test set was performed for individual rivers in periods of 12, 60, and 120 months. The aggregated forecast plots are shown in Figure S16.

The predictions of air temperature values within each horizon are shown in Table 6.

An analogous prediction was made taking into account all the data in the form of a test set. Air temperature predictions for the 120-month period for each river catchment are shown in Figure 9.

Based on the results provided in Figure 9, it can be concluded that the air temperature in each individual river basin will increase by 0.1 ◦C to 0.4 ◦C. In the case of the Styr River and the Sula River catchments, a trend of further variability decrease is observed. Air temperature variability for the Prut River basin will remain unchanged in the period 2021–2030.

**Table 6.** Forecast by air temperature model with train-test data grouped by horizon.


Explanations: N—sample size, M—mean, SD—standard deviation.

**Figure 9.** Combined plot with partial actual data (2015–2020) and predicted data (2021–2030) for the air temperature model trained with all data. Explanations are the same as those for Figure 7.

#### **4. Discussion**

The results of the monthly flow forecasts made for three rivers in Ukraine show a different fit of the model to the data, which is underlined by the measures of the model's performance. To some extent, similar patterns of the long-term variability of the monthly flow were obtained. However, forecasts of flows were different, which to a greater or lesser extent is related to the multiyear variability and different forecasts of precipitation and air temperature.

#### *4.1. Long-Term Changes of Monthly River Flow under Different Precipitation and Air Temperature Conditions*

In all the studied rivers, a variable dynamic of flows was observed in the years 1961–2020, with a clearly marked decrease in monthly flows in the last decade of 2010–2020. The last decade of low flows on the Prut–Yaremche and the Styr–Lutsk rivers was preceded by their noticeable increase in the earlier 2000–2010 decade. On the other hand, in the case of the Sula–Lubny river, a continuous decrease in monthly flows had been observed since the end of the 1990s (after a 30-year period of flow increase since the beginning of the 1960s) with a global minimum in the 2010–2020 decade. These changes are confirmed by the dynamics of flow trends, showing alternating periods of increase and decrease in flows over the decades of the analysed multiannual period (Figure 5).

The variability of river flows is partly related to the trends of precipitation and air temperature in the studied catchments. The decrease in flows observed on the Prut River in 1980–1990 corresponds to the lowest precipitation recorded in the 1980s–1990s. On the other hand, the observed small increase in precipitation during 1961–2020 was accompanied by a small decrease in river flow. One of the reasons for this discrepancy was a significant increase of air temperature and a corresponding increase of evaporation from the catchment area, which is also documented by the research of Vyshnevskyi and Donich [53]. Similar tendencies of climate parameters and river-flow changes were obtained for the Styr and the Sula rivers. Despite a small increase of precipitation, water runoff tends to decrease. This decrease is bigger than that of the Prut River. In our opinion, such decrease of flow of plain rivers may be explained by two factors: an increase of water intake during warm weather (in particular due to the irrigation of farming land) and an increase of evaporation from ponds and reservoirs. As shown in the paper by Vyshnevskyi [64], the dependence between evaporation and air temperature is nonlinear: a small increase in air temperature results in a significant increase of evaporation.

For all the studied rivers, a similarity was identified in the seasonal fluctuation component and the level of the random component of monthly flows. A detailed analysis of random flow fluctuations for selected rivers has not been carried out; therefore, research on these issues should be continued in the future.

#### *4.2. Forecasting of Monthly River Flow with Precipitation and Air Temperature Forecasts*

Results of the forecasts for 2021–2030 for the studied catchments showed that the dominant pattern is as follows: a decrease in flows is accompanied by a decrease in precipitation and an increase in air temperature (the Prut and the Styr). The differences relate to the magnitude of the flow change of approximately 6% compared with the 2010–2020 decade for the mountainous Prut River (or even no change in flow) and about 14% for the Styr River, predominantly in the uplands. On the other hand, the forecasts for the Sula River provide for an increase in flows (16–23%), with an increase in precipitation and air temperature, which is what makes it different from other rivers. It may be assumed that with an increase in precipitation in the Sula River catchment area, soil retention (underground retention) will also become more important in shaping the river flows, which will contribute to levelling the differences between high spring flows and low summer–autumn flows. The predicted changes in the flows of the Prut, the Styr, and the Sula rivers correspond to the forecasts in other regions of Ukraine and Europe. Over the past decades, climate change has been observed in Europe, and it has been emphasized that these changes will be stronger in the future [65–70]. Many studies on the impact of climate variability and change on water resources in Europe, using hydrological models on a continental and global scale, predict lower flows, with increasing air temperature and lower or no significant changes in precipitation [71–74]. However, the various types of changes are not evenly distributed either in time or in space. Didovets et al. [69] analysed changes in precipitation and temperature until the end of the 20th century in the catchment areas of Ukraine and determined their impact on the water availability of main river basins by using global hydrological models, which were compared with global climate models

(GCM) in two representative concentrations pathways scenarios: RCP 2.6 and RCP 8.5. A historical period (1971–2000) and two future periods (2041–2070 and 2071–2100) were considered, showing that changes in mean annual precipitation range between −14% and +10% and changes in the mean annual river discharge range from −30% to +6%, depending on RCP. The largest decrease in the mean annual runoff was forecast for the Pripyat, the Southern Bug, and the Dniester basins, reaching up to −30% by the end of the century in RCP 8.5, which was also confirmed by the present current study for the Styr River, which is a tributary of the Pripyat River.

#### *4.3. Assessment of the Forecasts of Monthly River Flow*

In the present research (regarding the Prut, Styr, and Sula rivers), attention was drawn to a significant RMSE value at each fitting step, which indicates high variability of the data that cannot be reasonably well fitted by retrospective data, regardless of the selected horizon length. This leads to the conclusion that the monthly data frequency for flow is too low to produce good fitting results and that it is worthwhile to use the daily frequency. Comparing the obtained results with other forecasts yields an emphasis on the fact that uncertainty regarding the forecasting of hydrological characteristics is still high, regardless of the location of the rivers and the adopted research methods. The variability of the river-flow time series can be influenced by various parameters, such as air temperature, precipitation, and evaporation, which makes accurate estimations and predictions of the flow almost impossible, the more so as it is further complicated by complexity, nonstationarity, and the nonlinearity of the river-flow phenomenon [75–77].

An analysis of the homogeneity and stationarity of an average monthly flows of Ukrainian rivers (305 water gauge stations, from the beginning of the observation until 2010), carried out with the use of hydrogenetic methods [78], has shown that most of the observation series are homogeneous and stationary. The exceptions include rivers with significant anthropogenic impact (river regulations, water intakes, deforestation, etc.). Observation series with a full cycle of long-term cyclic fluctuations (dry and wet phase) are stationary, whereas other observation series are quasistationary. Many studies indicate a large diversity of predictions and forecasts, especially those obtained for large rivers (often transboundary) and a high uncertainty as to the predicted hydrological effects in the future due to the visible climate changes and their impact on the hydrological characteristics of rivers, also in Ukraine [15,78]. The importance of studies on the strength of simulated changes and the effects on river flows has also been confirmed [79–82].

A staged forecasting process used in this study takes into account the gradient enhancement approach of the XGBoost model. In the approach applied to Ukrainian rivers, a single multioutput model was proposed that simultaneously forecasts over short- and longterm horizons. The generalization performance of the model over a variety of available data sets over time was assessed, and the hyperparameters which minimize the forecast error over the relevant forecast horizons were selected. The algorithm used allowed direct multistep-ahead forecasting with multiple time series (grouped or hierarchical) inspired by Bergmeir et al. [57]. As suggested by Sanders et al. [41] and Osman et al. [83], XGBoost models use all input data directly, without any preprocessing or selection, and a proper selection of the input fields is essential to improve the performance of XGBoost models. Usually, unavoidable errors may result from introduced model limitations or the nonstationarity of hydrological processes [84]. An additional difficulty may be the quality of hydrological and meteorological data, another being the calibration of parameters [85,86]. Yang et al. [87] suggest including additional hydrological information in river-flow forecasts to improve model performance. On the other hand, research conducted by Sanders et al. [41] shows that using more information on the catchment area and data from point measurements can lead to a reduction in the performance of forecasting models. According to Aghelpour et al. [88], consideration of the climatic indices (e.g., Southern Oscillation Index, Global Mean Temperature Index, North Pacific pattern, Pacific Decadal Oscillation, North Atlantic

Oscillation) as inputs in prediction models of river flows in Iran may improve monthly flow prediction by 24% on average.

Complex river-flow prediction models typically require multiple hydrological and climatological parameters as input, many of which may not be available for some locations, and the complexity of the flow processes makes it difficult to use physical models. Many estimates and forecasts of river flows have well documented that combined models (hybrid models) may show better performance compared with stand-alone models [89–91].

#### **5. Conclusions**

Long-term changes were estimated, and a forecast of monthly river flows was made for three rivers (the Prut, the Styr, and the Sula) in Ukraine with reference to the variability and forecasts of precipitation and air temperature. The years 1961–2020 were selected as the historical period, while the forecast horizons of 12, 60, and 120 months into the future were specified for this dataset, i.e., December 2021, December 2025, and December 2030.

The forecasting of monthly river flows was carried out in three stages with the use of input data sets characteristic for the XGBoost algorithm: train (80%) and test subsets (20%), using a model with train-test data, and forecasting by training with all data for forecast horizons (December 2021, December 2025, and December 2030). The adoption of this algorithm allowed direct multistep-ahead forecasting with multiple (grouped or hierarchical) time series, which may be of importance taking into account sometimesscarce measurements of daily river flows. This is one of the main advantages and a new, unique feature of the adopted forecasting method, which is characterized by high versatility of use, regardless of the spatial location of the river catchment areas. The direct forecasting approach used included the following steps: build a single multioutput model that simultaneously forecasts over short- and long-term forecast horizons, assess the generalization performance of the model over a variety of held data sets over time, and select the hyperparameters that minimize forecast error over the relevant forecast horizons and re-train.

The selected rivers flow through different physical geographical and climatic regions and represent different hydrological systems: the Prut, a mountain river; the Styr, an upland river; and the Sula, a lowland river. These features undoubtedly influenced their different degree of sensitivity and resistance to variability of climatic conditions and possible anthropogenic influences. The conducted analysis showed some differences in the trends of changes in monthly flows in individual decades of the period 1961–2020 and similarities in terms of the seasonal component and random flow fluctuations. An emphasis was put on the last decade, significant for future forecasts, in which a decrease in flow was determined for the Prut and the Styr rivers, whereas an increase was determined for the Sula River. Two patterns were obtained in the forecasts, which partially correspond to the results of flow forecasts obtained for Ukrainian rivers with the use of other methods and models. A decrease in flow, accompanied by a decrease in precipitation and an increase in air temperature until 2030, are expected according to the forecasts for the Prut and Styr rivers, whereas for the Sula River, the forecasts assume an increase in flow, accompanied by an increase in precipitation and air temperature.

The advantages of the model used for forecasting of monthly river flows in Ukraine, i.e., modelling many time series simultaneously using one model, include the simplicity of modelling, potentially more-robust results because of pooling data across time series, and solving the "cold start" problem when few data points are available for a given time series. Unlike the recursive or iterative method of creating forecasts in multiple steps used in traditional forecasting (such as ARIMA), direct forecasting involves creating a series of distinct, horizon-specific models. Although there are several hybrid methods for creating multistep forecasts, the simple, direct forecasting method used avoids the exponentially more-difficult problem of "predicting the predictors" for forecast horizons beyond one step in the future.

The obtained results of the forecasts, although approximate to some extent, provide information on future (in the next decade) changes in monthly flows of Ukrainian rivers (the Prut, the Styr, and the Sula) occurring under various environmental and climatic conditions, and they highlight the need to extend research into assessments of the impact of anthropogenic factors on the long-term dynamics of changes in river flows. In future studies of the daily flows, precipitation and air temperature should test the effectiveness of multistage forecasting with multiple, grouped, or hierarchical time series in the forecasting of hydrological and meteorological parameters. The research results can be used for the management and rational planning of water resources in the studied regions of Ukraine. They can provide the basis for the rational use of river water resources, such as for irrigation purposes, which is certainly a great challenge in the situation of increasingly frequent and longer periods of drought.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/resources11120111/s1, Figure S1. Custom validation subsets (grey areas) with skipped data (white areas) for river flow, Figure S2. Seasonal flow values grouped by rivers (1961–2020), Figure S3. Random flow values grouped by rivers (1961–2020), Figure S4. Precipitation time series distribution grouped by rivers' catchment areas (1961–2020). Note: Figure descriptions refer to non-parametric measures-mean (standard deviation), Figure S5. Precipitation trend values grouped by rivers (1961–2020), Figure S6. Precipitation seasonal values grouped by rivers (1961–2020), Figure S7. Precipitation random values grouped by rivers (1961–2020), Figure S8. Air temperature time series distribution grouped by rivers' catchment areas (1961–2020). Note: Figure descriptions refer to non-parametric measures-mean (standard deviation), Figure S9. Air temperatur trend values grouped by rivers (1961–2020), Figure S10. Air temperatur seasonal values grouped by rivers (1961–2020), Figure S11. Air temperatur and random values grouped by rivers (1961–2020), Figure S12. River flows predicted (solid line) vs. actual (dashed line) values grouped by horizons (purple line: the Prut River-Yaremche, turquoise line: the Styr River–Lutsk, yellow line: the Sula River-Lubny), Figure S13. Precipitation predicted (solid) vs. actual (dashed) values grouped by horizons (purple line: Prut River catchment. turquoise line: Styr River catchment. yellow line: Sula River catchment), Figure S14. 12, 60 and 120 step-ahead forecasts of precipitation (purple line: the Prut River catchment, turquoise. line: the Styr River catchment, yellow line: the Sula River catchment), Figure S15. Air temperature predicted (solid) vs. actual (dashed) values grouped by horizons (purple line: PrutRiver catchment. turquoise line: Styr River catchment. yellow line: Sula River catchment), Figure S16. 12, 60 and 120 step-ahead forecasts of air temperature (purple line: the Prut River catchment, turquoise line: the Styr River catchment, yellow line: the Sula river catchment). Table S1. Forecast error by river flow forecast horizon, Table S2. Global river flow forecast error, Table S3. Forecast error by precipitation forecast horizon and river catchment, Table S4. Global precipitation forecast error grouped by river catchment, Table S5. Forecast error by air temperature forecast horizon and river catchment, Table S6. Global air temperature forecast error grouped by river catchment.

**Author Contributions:** Conceptualization, R.G. and V.V.; methodology, R.G.; software, R.G.; data curation, V.V.; writing—original draft preparation, R.G. and V.V.; investigation, R.G.; validation, R.G.; visualization, R.G.; writing—reviewing and editing, R.G. and V.V.; discussion, R.G. and V.V.; supervision, R.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The source of data derives from the measurement and observation data of the Hydrometeorological Service in Ukraine.

**Acknowledgments:** The authors thank their universities for their support. The authors would also like to acknowledge the Hydrometeorological Service in Ukraine for the release of the input database.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Shmueli, G. To explain or to predict? *Stat. Sci.* **2010**, *25*, 289–310. [CrossRef]


**Renata Graf 1,\*, Tomasz Kolerski <sup>2</sup> and Senlin Zhu <sup>3</sup>**


**Abstract:** Forecasting ice phenomena in river systems is of great importance because these phenomena are a fundamental part of the hydrological regime. Due to the stochasticity of ice phenomena, their prediction is a difficult process, especially when data sets are sparse or incomplete. In this study, two machine learning models—Multilayer Perceptron Neural Network (MLPNN) and Extreme Gradient Boosting (XGBoost)—were developed to predict ice phenomena in the Warta River in Poland in a temperate climate zone. Observational data from eight river gauges during the period 1983–2013 were used. The performance of the model was evaluated using four model fit measures. The results showed that the choice of input variables influenced the accuracy of the developed models. The most important predictors were the nature of phenomena on the day before an observation, as well as water and air temperatures; river flow and water level were less important for predicting the formation of ice phenomena. The modeling results showed that both MLPNN and XGBoost provided promising results for the prediction of ice phenomena. The research results of the present study could also be useful for predicting ice phenomena in other regions.

**Keywords:** river freezing; Multilayer Perceptron Neural Network (MLPNN); Extreme Gradient Boosting (XGBoost); predictor variables; balanced accuracy; Poland

#### **1. Introduction**

Prediction of ice phenomena in rivers is an important element of hydrological regime analysis [1] and the assessment of the risk of ice jam type floods [2]. The changing thermal conditions of river waters during the winter season and the nature of river ice may significantly change the hydro-ecological and socio-economic aspects of the functioning of the river ecosystem.

Due to the stochastic nature of ice phenomena, their prediction is difficult, especially when data sets for rivers are sparse or incomplete. An additional complication is the scale of the event (local and regional scales) and the influence of numerous factors on the process of river freezing, e.g., meteorological (e.g., air temperature, solar radiation, wind velocity) [3,4], hydrological (e.g., flow rate, inflow and outflow conditions) [5,6], the complexity of interactions between hydroclimatic factors [7,8], hydraulic (e.g., trough crosssection geometry, river bathymetry, water table drop) [9] and thermodynamic factors (e.g., water temperature and thermal conductivity) [10,11]. Relations between river freezing and features of the hydrological regime, including flow, water state, and water temperature, are usually complex and non-linear, and are also spatially heterogeneous due to the variability of environmental conditions. In addition to the process that determines the number of occurrences of a given phenomenon, there is also a dichotomous process determining whether it has a chance of occurring in a given period [12]. This task is further complicated

**Citation:** Graf, R.; Kolerski, T.; Zhu, S. Predicting Ice Phenomena in a River Using the Artificial Neural Network and Extreme Gradient Boosting. *Resources* **2022**, *11*, 12. https:// doi.org/10.3390/resources11020012

Academic Editor: Diego Copetti

Received: 22 November 2021 Accepted: 24 January 2022 Published: 26 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

by the fact that ice phenomena occur in three phases: freezing of the river (first symptom of ice), permanent ice cover, and the disappearance phase when an ice floe is formed and related phenomena appear, such as ice jams, which often lead to winter floods. However, the full freezing cycle is not always recorded for rivers.

The analysis of time series relating to ice phenomena allows for the determination of the frequency and duration of their occurrence and the tendency of changes over time, and also for an assessment of the ice phases, which provides a good background for the characteristics of the freezing process in many regions. However, it is not sufficient for their prediction and forecasting [13]. Although Shulyakovskii [14] has developed a manual for forecasting the freezing of inland rivers and lakes, there are few studies related to this topic, especially works dealing with the prediction of ice phenomena at various stages of their occurrence. The problem most frequently discussed is the prediction of ice jams on rivers and their consequences in the form of ice jam floods. The theoretical model of river ice jams was developed by Uzuner and Kennedy [15]. Existing forecasts of ice extent are most often based on the location of the 0 ◦C isotherm [16]. Good results in this regard have also been obtained from observations of river ice ranges carried out with the use of satellites. Remote sensing is useful for the monitoring of ice characteristics, such as different types and thicknesses of ice or ice cover, and for tracking the progress of the breakup of ice jams, which can help predict the location and timing of ice blockages [17,18]. However, the results of field studies and analyses of satellite images do not always provide accurate data for forecasting ice and scenarios of changes in ice dynamics [19].

Prediction models for ice phenomena are usually limited to the empirical or the stochastic due to the difficulties in applying deterministic models. The methods used to predict ice phenomena (e.g., ice jams) include empirical single-variable threshold analyses, logistic regression [2,20], and discriminant function analysis [21]. Many numerical models have been developed to simulate ice formation on rivers [22,23]. According to Wang et al. [24] and Beltaos [25], a better understanding of physical processes has increased the possibility of developing more accurate numerical models of ice jams and ice jam floods in rivers, e.g., the public-domain river-ice RIVICE model [10], the DynaRice model, a two-dimensional coupled hydrodynamic and ice dynamic model [23], and hydraulic models [19]. An interesting ice jam flood forecasting system that considers requirements for the real-time predictions of water, ice, and sediment transport, was developed for the lower Odra River [26]. The prediction of ice phenomena was also carried out using teleconnection indices, as presented by Sutyrina [27] in relation to spring ice phenomena in lakes and reservoirs (including for Lake Baikal).

In the prediction of ice phenomena, machine learning methods are used less frequently, although they have already been utilized widely in forecasting time series of hydrometeorological data [28–30]. Artificial neural networks (ANNs) have been used to forecast freezing conditions in rivers [31,32] and predict ice jams [4,33]. For example, Chokmani et al. [34] estimated the thickness of ice using artificial neural networks (ANNs), while Hu et al. [35] predicted the disappearance of ice phenomena using hybrid artificial neural networks. Furthermore, fuzzy logic systems have provided favorable results as regards ice phenomena forecasting and its effects [20,36]. For example, Zhao [3] predicted the breakup date of flood ice using a wavelet neural network (WNN) model. Whereas Yan and Ding [37] proposed a predictive model of ice formation based on a dynamic fuzzy neural network (D-FNN) in combination with a particle swarm optimization (PSO) algorithm. The significant advantages of artificial neural networks over standard statistical classification methods consist in their ability to adapt to data of different formats and configurations [32,34]. Ensemble machine learning methodologies, including resampling methods (bagging, boosting, and dagging), model averaging, and stacking, are used in the solving of problems related to simulation and prediction in hydrology [38,39].

The increase in the amount of hydrological and meteorological data makes it more difficult not only to select the methods for their analysis but also to choose predictive and prognostic models so as to maintain both their legibility and accuracy [40]. For the integrated management of an aquatic ecosystem, it is necessary to determine how the thermal–ice regime of the river will develop and change in the future, considering global climate change and local conditions in particular [41]. The identification of the most important hydrological and thermal variables influencing the course of ice phenomena in a river may result in more accurate forecasts of the freezing process.

The main goal of the present study is to predict ice phenomena in a river with the use of the Multilayer Perceptron Neural Network (MLPNN) and Extreme Gradient Boosting (XGBoost) algorithms, which belong to the group of machine learning methods. MLPNN is one of the most widely used ANN models in the field of hydrology [4,7,8]. According to Zounemat-Kermani et al. [42], the boosting methods (e.g., boosting, AdaBoost, and Extreme Gradient Boosting) are becoming more and more effective for modeling and forecasting water quality, runoff, sediment transport, groundwater, flooding, and drought. One of the advantages of XGBoost compared to neural networks is the ability to assess the importance of predictors in the model, and in this study, by employing the XGBoost model, we can assess the dominant factors controlling the dynamics of ice phenomena in the studied river. The objective of the present research was to show the predictability of the selected models and explain spatial differences in terms of the predictors: air temperature (Ta), water temperature (Tw), water level (H), and river flow (Q), as well as the ice phenomenon of the previous day and of the month of occurrence of the phenomenon. The predictions will be carried out using the example of the Warta River in Poland (Central Europe), which is a river of great economic significance and considerable natural value. The results of the study are important for determining the range of intensification of thermal and hydrological ice phenomena variables and the conditions under which their reduction will occur.

#### **2. Study Area**

The Warta River is a tributary of the Odra River and the third longest (808 km) river in Poland (Figure 1). Its catchment area (area 54,500 km2) is characterized by a significant diversity of topography and terrain and climatic and hydrological conditions [43]. Within the Warta Water Region there are three main types of relief: old-glacial in the southern part, young-glacial in the northern and central parts, and upland, south of Wielun.

The catchment belongs to nine out of 28 climatic regions designated in Poland by Wo´s [44]. The average annual air temperature ranges from 7.5 ◦C in the north to 8.5 ◦C in the west. In the coldest month—January—the average temperature ranges from −1.2 ◦C (in the west) to −2.5 ◦C (in the southeast). Annual rainfall totals in the study area are diverse and range from 520 mm in the Kujawy region (in the northeast) to 675 mm in the south. A regional differentiation of features of the hydrological regime has been observed along the analyzed section of the Warta [45]—from a medium-developed (the upper and lower course of the river) to a highly-developed (along the section from Nowa Wies (Nowa Wie´s Podgórna to Poznan) nival regime (Figure 1). The average dates of appearance of ice phenomena on the Warta River, as well as the dates of their disappearance, vary. Research by Graf et al. [46] for the 1991–2010 observation series showed that the earliest ice phenomena occurred in the third decade of December (about 45% of the total number of observations) and the latest in the first decade of January. The disappearance of ice phenomena is usually observed from the end of January to the end of March [47], while about 30% of observations are made in the third decade of February. Most days with ice phenomena on the Warta River are in January (41% of observations), and the most common form of ice is frazil ice (46% of phenomena) and ice cover (30%).

**Figure 1.** Study area and the locations of water gauges and meteorological stations of the Institute of Meteorology and Water Management—National Research Institute (IMGW-PIB, Warsaw, Poland).

#### **3. Materials and Methods**

Predictions of ice phenomena and their numerical descriptions were performed based on daily data on the number of occurrences (the number of days on which the phenomena were observed) and the nature of ice phenomena, and on air temperature (Ta), water temperature (Tw), water levels (H), and river flow (Q) for the years 1983–2013 from the Central Database of Historical Data of the Institute of Meteorology and Water Management— National Research Institute (IMGW-PIB) in Warsaw, Poland (Figure 1). The observation series includes data for the period after 1980, when changes in water temperature in rivers and further consequences, including the lower incidence of ice phenomena, were revealed in response to the sudden climate change associated with changes in CRS (climate regime shift). The regime shift of the late 1980s is a well-documented example of CRS in Poland [48].

Use was made of data from eight water gauges on the Warta River (Bobry, Sieradz, Uniejow, Nowa Wies (Nowa Wie´s Podgórna), Srem, Poznan, Skwierzyna, and Gorzow Wielkopolski) and seven meteorological stations (Wielun, Sieradz, Koło, Słupca, Kornik, Poznan, and Gorzow Wielkopolski). Data have been presented in relation to the hydrological year, which in Poland lasts from 1 November until 31 October.

#### *3.1. Classification of Ice Phenomena*

The full ice cycle of the river includes forms of ice phenomena observed within the IMGW-PIB water gauge network: frazil ice, border ice, border ice and frazil ice, frazil ice jam, ice cover, ice floes, ice floes and border ice, ice floes and frazil ice, and ice jams. For modeling and predicting ice phenomena, these were grouped into three basic categories: (1) river freeze-up, (2) stable ice cover, and (3) breakup of ice cover—the disappearance of ice (Table 1). The joining into classes is not random, and indeed follows from the order in which ice phenomena appear on the river depending on the thermal and hydrological conditions of the winter season. Each observed ice phenomenon was assigned to the date of occurrence (month and year).


**Table 1.** Classification and grouping of ice phenomena.

In each of the mentioned phases, characteristic fluvial processes occur as a result of the appearance of various forms of ice. Strictly defined forms, such as frazil ice jams or ice jams, are ephemeral forms. According to data from the IMGW-PIB, in the analyzed period, the Warta River had only several days with frazil ice jams and ice jams, which accounted for 0.1% of all observed ice phenomena. Frazil ice jam embolism occurred on the Warta River in Skwierzyna (three days with jams) and Poznan (one day with jams). Ice jams on the Warta occurred in Srem (three days with jams) and Uniejow (two days with jams).

It should be emphasized that ice jams are not a common event in the Warta River. They are among the characteristic features of the river's morphology, which makes the reaches of the river susceptible to the formation of frazil ice [46]. In addition, climatic change serves to significantly decrease the intensity of ice jams by increasing the temperature in the vicinity of supercooled water and thereby prohibiting the formation of ice jams. Furthermore, the Warta River is strongly impacted by anthropogenic activity. Regarding the features of its morphology, the Warta River consists of various bed slopes, from mild to steep, which follow each other in a way that affects the pattern of ice formation. Surface ice is observed mainly over the milder sloped beds. Less surface ice is observed over the steeper slopes, while more suspended frazil particles are present. These ice particles may accumulate under the cover of the following flat sections, forming hanging dams at the inlet of the low-sloped sector. There is no specific data for the Warta River regarding hanging dams that would allow us to distinguish between the ice cover itself and frazil depositions with a greater degree of certainty.

#### *3.2. Data Preparation*

Predictions of ice phenomena were performed based on daily data on the number of occurrences (number of days with the phenomena) and the nature of ice phenomena and on air temperature (Ta), water temperature (Tw), water levels (H), and river flow (Q), as well as ice phenomena of the previous day (the day before occurrences of ice phenomena from classes 1, 2, or 3, or from class "none") and the month of individual phenomena (six months of the hydrological winter half-year XI-IV). The choice of input variables was not accidental. The hydroclimatic factors and thermal conditions are important predictors for the process of ice phenomena formation.

To improve the predictability of the tested models and accelerate the process of simulation convergence (in particular as regards artificial neural networks (ANNs)), before inputting the data into Multilayer Perceptron Neural Network (MLPNN) the variables were normalized by converting their values into standardized values (so-called Z-scores) [29]:

$$z\_{l} = \frac{\mathfrak{x}\_{l} - \underline{\mathfrak{x}}}{sd\_{\mathfrak{x}}} \tag{1}$$

where *xi* is the *i*-th value of *x*, *x* is the mean of *x*, and *sdx* is the standard deviation of *x*. The resulting variable *z* has a mean of zero and a standard deviation of one while retaining all the properties of the original variable. The following variables were transformed: Ta, Tw, H, Q, the day of the month (mon.), and the year (Y).

Additionally, ice phenomena of the previous day, encoded in four columns with the one-hot method (zero for the absence of a given phenomenon, and one when it occurred), and the month of individual phenomena, also encoded with the one-hot method, were introduced into the models. Encoding ice phenomena using the one-hot method with the addition of labels allows the assignment of your own characteristics to specific phenomena, showing the similarities between phenomena or the features that make them different. The ice phenomena data used are treated as categorical variables (also called nominal variables), that is, they represent the types of data that can be broken down into groups. In the tested example, three classes were distinguished (Table 1). However, the categories cannot be ordered from highest to lowest. In the classification methods these variables—as target variables (the ones we want to predict)—are usually converted to numerical form using one-hot coding.

#### *3.3. Descriptive Statistics of the Frequency of Ice Phenomena*

The research methodology included several stages. The development of predictive models for ice phenomena on the river was preceded by a statistical description of ice phenomena and their changes in the studied period.

The statistical description of ice phenomena included an analysis of the frequency of ice phenomena in the set of analyzed data, which was determined separately for each measuring station and assuming the classification of ice phenomena into three classes, as presented in Table 1. The next stage concerned the analysis of ice phenomena as sequential phenomena. For this purpose, cross-tables were made comparing ice phenomena from the current day with phenomena from the previous day. In the last step, the relationships between the classes of ice phenomena and air temperature, as well as water temperature, water level, and river flow, were analyzed. For this purpose, box and violin plots were made for the distribution of these parameters for each class of ice phenomenon. The violin plot is a combination of a box plot and a density plot, thus showing more details of data distribution, especially the kernel density distribution [49]. As a result, the problem of overlapping the traditional density plot, which is difficult to identify, is eliminated. Wider sections of the graph signify the higher probability of occurrence of certain values, while narrower sections denote lesser probability. According to Hintze and Nelson [50], the violin plot is used to visualize quantitative and qualitative data, including those that do not conform to the normal distribution, and to define the data structure. Like box plots, violin plots are used to present a comparison of variable distribution (or sample distribution) across different categories.

#### *3.4. Prediction Models*

Figure 2 depicts the stages of research activities in brief. The current research was carried out in three stages in total. To begin with, the data that had been cleaned, standardized, and adapted to the needs was referred to as prepared data. The second step was to use the R tools to test model predictions using both XGBoost and MLPNN methods. In the

prediction of ice phenomena, the following formula was used (according to the choice of the MLPNN and XGboost algorithms):

$$\begin{array}{c} \text{ice\\_0 + ice\\_1 + ice\\_2 + ice\\_3 \sim} \\ \text{Ta} + \text{Tw} + \text{Q} + \text{H} + \text{D} + \text{Y} + \\ \text{day\\_before0} + \text{day\\_before1} + \text{day\\_before2} + \text{day\\_before3} + \\ \text{mo1} + \text{mo2} + \text{mo3} + \text{mo4} + \text{mo5} + \text{mo6} \end{array} \tag{2}$$

where ice\_0 means no ice phenomena, ice\_1–3 is the classes of ice phenomena (classes 1–3 adopted on the basis of the classification and grouping of ice phenomena presented in Table 1), Ta is air temperature, Tw is water temperature, Q is river flow, H is water level, D is day of the month, Y is day of the year, day\_before0 is day before the day with no ice phenomena (class "none"), day\_before1–3 means the day before occurrences of ice phenomena from classes 1, 2, 3 and mo1–6 means months in the winter half-year (November– April, according to the hydrological year).

**Figure 2.** Stages of research activities.

The training and test sets were created using the stratified sampling algorithm, with the year and month variables functioning as layers. The process of determining datasets is detailed in the description Evaluating the Predictions. The confusion matrices were formed on the basis of the second stage of the activity and were then used as inputs for the third stage, which involved evaluating the performance of the XGBoost and MLPNN methods.

#### 3.4.1. The Multilayer Perceptron Neural Network (MLPNN)

The most commonly used type of neural network method is the multi-layered perception method. In this method, the signal is passed to a one-way loop-free input-to-output network. Neither neuron acts on itself. This architecture is referred to as feed-forward, and consists of multiple inputs, hidden layers, and an output, as shown in Figure 3.

The first model used to predict ice phenomena was the Multilayer Perceptron Neural Network (MLPNN), which included an input layer, one hidden layer, and an output layer, and is one of the most widely used ANN models in the field of hydrology [4,7,8]. The input layer, which comprises the predictors, does not perform any calculations. The hidden layer is made of artificial neurons. A single hidden neuron 'collects' activations from each neuron of the input layer and calculates the weighted sum of the input variables. Each hidden layer neuron is connected to each input layer neuron. The hidden layer neurons then perform

a non-linear transformation of the weighted sums using an activating function and pass the results to the output layer, which in this application is represented by ice phenomena. A neural network of this type with an output variable *Y* and containing n neurons in the hidden layer can be expressed as follows [51]:

$$Y = f\_2 \left[ \sum\_{j=1}^{n} w\_{jk} \left[ f\_1 \left( \sum\_{j=1}^{n} \mathbf{x}\_i w\_{ij} + \delta\_j \right) \right] + \delta\_0 \right] \tag{3}$$

where *xi* is the value of the input variable *i, wij* is the weight (synapse) between the input variable *i* and the hidden neuron *j*, *δ<sup>j</sup>* is the bias of the hidden neuron *j*, *f*<sup>1</sup> is the sigmoidal function constituting the activation function for hidden neurons, *wjk* is the synapse between the hidden neuron *j* and the output neuron *k* (here *k* = 4), *f*<sup>2</sup> is also the activation sigmoid function, and *δ*<sup>0</sup> is the bias of the output layer neuron. The use of the sigmoidal function as an activation function for neurons of the output layer ensured that the predictions would be obtained from the model.

**Figure 3.** Feed-forward multilayer perceptron architecture.

To estimate the weights and biases, the neural net package [52] and implemented elastic back propagation [53] were used. Cross entropy was used as a function of cost. Models with three, four, five, and six neurons in the hidden layer were calculated for each station.

#### 3.4.2. The Extreme Gradient Boosting (XGBoost) Model

The second model tested was the Extreme Gradient Boosting (XGBoost) implemented by Chen et al. [54]—also in the form of the XGBoost library for the R platform. The gradient boosting machine is a team learning technique based on decision trees. A decision tree generates an output variable estimate based on optimized predictor thresholds that divide the data into multiple groups. The gradient boosting algorithm in each subsequent step aims to reduce the prediction error of the previous step. Technically, in each subsequent step the algorithm estimates the parameters of the model whose purpose is to predict the residuals (prediction errors) of the model estimated in the previous step. The objective function (*J*) in round *t* (step *t*) is given by Equation (4) [54]:

$$J^{(t)} = \sum\_{i=1}^{n} l(y\_i \quad \hat{y}\_{\\_i}) + \sum\_{k=1}^{K} \Omega(f\_k) \tag{4}$$

where: *l* is the training loss, Ω is regulations, *fk* is the function of the K–tree. In this study, *yi* is the observed ice phenomena and *y*ˆ*<sup>i</sup>* is the obtained final prediction value.

In the present study, decision trees with a maximum depth of five nodes were used. Formally, a tree is any consistent acyclic graph, i.e., a graph that does not contain cycles. The multi-class log loss function was used as the cost function. Predictions of *Y*(*t*) from the model for iteration *t* are obtained from Equation (5) [39,54]:

$$Y^{(t)} = \sum\_{k=1}^{t} f\_k(X) = Y^{(t-1)} + f\_t(X) \tag{5}$$

where *X* is the predictor or set of predictors and *fk* is the function that returns the predicted values of the predictors. The second part of the equation shows explicitly that the algorithm prediction in the iteration *t* is the sum of predictions from the *t* − 1 iteration and the new predictions from the *t* iteration. In XGBoost, the function *fk* consists of classification and regression trees that enable the modeling of arbitrary nonlinear relations and the prediction of variables of any nature (Figure 4).

**Figure 4.** A general architecture for XGBoost.

One of the advantages of XGBoost compared to neural networks is the ability to assess the importance of predictors in the model. The importance of a predictor for regression and classification trees in the gradient boosting algorithm is defined as the profit that the predictor contributes to the entire model by using it to create successive branches of the tree. In this study, by employing the XGBoost model, we can assess the dominant factors controlling the dynamics of ice phenomena in the studied river.

#### *3.5. Evaluating the Predictions*

To assess the predictive power of the tested models, cross-validation and four goodnessof-fit metrics were used. Cross-validation was performed by training the models on the available data (training data) and then calculating predictions and goodness-of-fit metrics for the data on which the algorithms were not trained (test data). The XGBoost model was taught on 70% of the training set, and the prediction model was tested on 30% of the test set. The ANN model was taught on the first 50% of the sample, and the prediction model was tested on the remaining 50%. The test and training sets were created using the stratified sampling algorithm. The year and month variables were used in the form of layers. This was done specifically so that the training and test sets had a comparable number of observations within each year and month included in the analysis. Such divisions are in line with the general practice of evaluating machine learning algorithms [51].

The test and training sets were created using the stratified sampling algorithm, with the year and month variables functioning as layers. As a result, the test and random sets had a comparable number of observations within each year and month.

Four metrics were used as goodness-of-fit metrics, calculated separately for each class of ice phenomena: sensitivity, specificity, precision, and weighted validity [55]. For ease of interpretation of these statistics, consider the following cross tables (Table 2), where the letters A–D represent the counts:

**Table 2.** Goodness-of-fit metrics.


Explanation: A—TP (True-Positive), the number of true positive predictions, i.e., correctly classified examples from the selected class, B—FP (False-Positive), the number of false-positive predictions, i.e., examples incorrectly assigned to the selected class when in fact they do not belong to it, the so-called false alarm, C—FN (False-Negative), the number of false-negative predictions, i.e., misclassified examples from this class, i.e., a negative decision while the example is positive (the so-called error of miss), D—TN (True-Negative), the number of truly negative predictions, i.e., examples correctly not assigned to the selected class (correct rejection).

The statistics used are defined by the formulas:

Sensitivity = A/(A + C) (6)

$$\text{Specificity} = \text{D} / (\text{D} + \text{B}) \tag{7}$$

$$\text{Precision} = \text{A} / (\text{A} + \text{B}) \tag{8}$$

Balanced Accuracy = (Sensitivity + Specificity)/2 (9)

Sensitivity TPR (True-Positive Rate) is a measure of "reach" (coverage, "reaching") that indicates the percentage of the positive class that has been covered by a positive prediction [56]. Specificity TNR (True-Negative Rate) is a measure of "coverage" that indicates the percentage of the negative class being covered by the negative prediction. Theoretically, Sensitivity (TPR) and Specificity (TNR) are independent measures, however. in practice increasing sensitivity often leads to a decrease in specificity [55]. Precision, referred to as the Positive Predictive Value (PPV), is a measure of precision that indicates how confidently we can trust positive predictions, i.e., the percentage of positive predictions that are positive. The confidence interval for the three distinguished measures is built based on the Clopper–Pearson method for a single proportion. Accuracy is the proportion of correct predictions with a set of test data. It is the ratio of the number of correct predictions to the total number of input samples. In turn, Balanced Accuracy is the arithmetic mean of the recall for each class. The closer the value is to 1, the better the prediction. However, exactly 1 indicates a problem that may be typically labeled as over-fitting. For highly unbalanced classification problems, as in the case of the analyzed data, balanced accuracy is particularly useful, because this statistic depends on both the level of correct prediction of a phenomenon and the level of prediction of the absence of a phenomenon.

Data analyses and operations were performed using the R 4.02 statistical environment [57]. The analyses and the necessary data restructuring, as well as the visualization of the data and the results of the analyses, were performed using the basic functions of the R environment and dedicated libraries for a given type of algorithm. The libraries used are cited in the corresponding analysis.

#### **4. Results**

#### *4.1. Probability of Occurrence of Ice Phenomena*

The frequency of ice phenomena on the Warta River in the analyzed period has been presented in Table 3. At the majority of measuring stations, ice phenomena from class 1 were observed on slightly more than 10% of days, while the frequency of occurrence of phenomena from class 2 varies from about 1.5% to over 8%. Ice phenomena from class 3 (breakup of ice cover—disappearance of freezing) were the least frequently observed. At each measuring station, this class was observed on less than 1% of days.


**Table 3.** The frequency of ice phenomena.

\* No—no ice phenomena.

The probability of occurrence of ice phenomena in specific months of the year (in the cold semester of the hydrological year) has been presented in Figure 5. The probability of occurrence of ice phenomena from class 1 is highest in the months of December and January. Class 2 events are most likely to occur in January and February, whereas the greatest probability of the breakup of ice cover and disappearance of freezing (class 3) is associated with the month of January; with February in Sieradz, and with March in Poznan.

The results of the analysis of ice phenomena as sequential phenomena have been presented in Figure 6. The cross tables compare ice phenomena from the current day with the phenomena of the previous day. It was noted that each class was most often preceded by a phenomenon from its class. Additionally, there was often no ice at all at the river stations the day before the occurrence of ice phenomena from class 1. In a small percentage of days, ice phenomena from class 1 preceded class 2 events. Class 3 occurrences were regularly preceded by phenomena from classes 1 and 2 (Figure 6).

#### *4.2. The Relationship between Ice Phenomena and Hydrological Conditions and Thermal Variables*

The assessment of the relationship between various classes of ice phenomena and thermal conditions and hydrological factors has been presented in the form of violin plots of the distribution of these parameters for each class of the phenomenon. Figure 7 shows the differentiation of the variables with respect to the water gauges. For each water gauge station on the Warta River, the differentiation of the occurrence of ice phenomena in relation to air (Ta) and water (Tw) temperatures as well as water level (H) and flow (Q) was presented. As drawn, the graphs indicate certain regularities of occurrence of ice phenomena on the river.

The phenomena from the first stage of ice (border ice, frazil ice) are characteristic of the conditions of poor cooling of the water and mild flow, i.e., for the months of November and December. Although frazil ice requires a significant subcooling of the water and an effective dissipation of the heat of solidification, it can form particularly abundantly during strong, cold winds, even if the air temperature drop is insignificant (even at a few degrees below 0 ◦C). The analysis of the data showed that the ice phenomena from the first phase occur even at a water temperature of the Warta River of 0.2–0.8 ◦C, and ice cover is maintained at a water temperature of 0.2 ◦C and at negative air temperatures, which is understandable. In the case of lowland rivers, which also include the Warta River, the ice cover expansion phase, due to low flows and falls, lasts the longest, and its formation is favored by the persistence of negative air temperatures for a long time [46]. The period of ice cover disappearance as a result of an increase in air temperature occurs on the river in stages, as a result of which an ice floe is created that moves downstream (ice procession). The flow of ice floes in the river usually accelerates the cracking of the ice cover caused by the rise in the water level in the spring.

**Figure 5.** Probability of occurrence of ice phenomena (classes 1–3 and none) as a function of the month for the water gauge stations on the Warta River. Note: Water gauge stations are labeled in the order (**a**–**h**), according to their location on the river (from upper to lower course).

**Figure 6.** Probability of the order of occurrence of ice phenomena classes (none, C1, C2, C3) as a function of the ice phenomenon of the previous day for the water gauge stations on the Warta River. Note: Water gauge stations are labeled in the order (**a**–**h**), according to their location on the river (from upper to lower course).

**Figure 7.** Distributions of the relationship between classes of ice phenomena (none, C1, C2, C3) and air and water temperatures (Ta, Tw), the water level (H), and river flow (Q).

The distribution of ice phenomena concerning water temperature has a distinct character. In this case, the distribution for phenomena classes 1–3 is unimodal and has one high "peak" at very low water temperatures, which indicates the typical regularity of the occurrence of the first ice on the river. Considering the ice cover, the distribution partly takes the form of a slanting distribution with a long tail, which can be seen in the graph for the water gauges of Uniejow, Nowa Wies, and Srem (Figure 7).

As regards the relationship between ice phenomena and air temperature, distribution becomes more diverse depending on the class of the phenomenon and the location of the observation post. For class 1, the distribution is predominantly unimodal. For the majority of measuring stations, distribution is asymmetric and has features of skewed distribution. At the Sieradz and Skwierzyna stations, this distribution shows a tendency to bimodality, which would suggest the presence of two characteristic periods of air temperature and thus favor an increase in the probability of occurrence of ice phenomena from the first phase in these locations. In the case of class 2 (permanent ice cover), distributions at almost all stations are unimodal with a clear skew towards very low air temperatures, which strongly suggests that the probability of ice cover is related to the accumulation

of days with negative air temperature. The exception is the Poznan water gauge, for which distribution has features similar to the bimodal distribution (Figure 7). For class 3 (breakdown of the ice cover), the distribution has features typical of unimodal distribution and is focused on an air temperature ranging from 0 ◦C to a few degrees above zero. This form of relationship is typical of most water gauges on the Warta River. Finally, as regards the water gauges in Nowa Wies and Skwierzyna, the skewness of the distribution increases, and this points to an increase in outliers.

The distribution of ice phenomena from class 1 (formation of ice phenomena) concerning the water level displayed predominantly bimodal features (Figure 7). In the case of the Bobry and Sieradz stations, the distribution shows features of asymmetry and develops as a skewed distribution. In the case of distributions with two or more mods, the widest sections of the violin diagram indicate the greatest probability of observing ice phenomena on the river at a low and medium water level. However, additional periods with a specific water level on the Warta River (states above the average) at which ice phenomena will occur under favorable river thermal conditions are not excluded. The bimodal distribution indicates that the distribution of ice phenomena in this relationship is unstable or very variable. Distribution displays similar features in the case of the relation of the ice cover (class 2) to the state of the water, which is also bimodal at most stations (Figure 7). The distribution shows similar features as regards the relation between the ice cover (class 2) and the state of the water, which, too, is bimodal at the majority of water gauges. The exceptions here are Sieradz and Poznan, for which a unimodal distribution with a specific skewness has been identified. For ice phenomena from class 3, the relationship with the water level shows different types of distributions: unimodal (Poznan and Gorzow Wlkp.) and biomodal (Nowa Wies and Skwierzyna). In the case of Sieradz, the distribution is flat, while in Gorzow Wlkp. it is strongly skewed.

The distribution of the relationship between ice phenomena from class 1 and river flow is bimodal at Nowa Wies, Poznan, and Skwierzyna and unimodal at other stations (Figure 7). In the case of Uniejow, Srem, and Gorzow Wlkp., distribution is also strongly skewed. The greatest probability of occurrence of ice phenomena in the initial period of the Warta River's freezing is associated with the low flow of the river. In the case of permanent ice cover (class 2), the distribution is unimodal at all water gauges, except for Nowa Wies, where it exhibits features of bimodality. This means that the distribution of ice in this relationship is relatively stable along the entire river. However, as regards ice phenomena from class 3, the distribution at certain stations has unimodal (Uniejow, Srem, and Gorzow Wlkp.) or bimodal (Sieradz, Nowa Wies, and Skwierzyna) features.

#### *4.3. Predicting Ice Phenomena*

The results of predictive modeling have been presented for three sections of the Warta River: the upper course (Bobry, Sieradz, and Uniejow water gauges)—Table 4; the middle course (Nowa Wies, Srem, and Poznan water gauges)—Table 5; and the lower course (Skwierzyna and Gorzow Wlkp. water gauges)—Table 6. In most of the analyzed instances, the predictive power of the tested models was comparable, and the differences in metrics between models were inconsiderable.

In the upper section of the Warta River (Bobry station), the MLPNN with four hidden units (NN4) was the best among the models, as indicated by the highest values of "balanced accuracy" (BA) statistics for ice phenomena from class 2 (BA = 0.971), and for the "no ice phenomena" class (BA = 0.933), and the second-highest value of statistics for class 1 (BA = 0.913) in the test set (Table 4). The XGBoost model predicted ice phenomena to a comparable extent. It exhibited a similar "balanced accuracy" profile, but one slightly infe-rior to NN models. Class 3 was too small in terms of abundance for the model to success-fully learn the relationship between the class and the predictors in this dataset. For the Sieradz station, it is difficult to identify the model with the highest predictive power (Table 4). The XGBoost model and the NN3–NN5 models successfully predicted each class of ice phenomena. From the NN models, the model with four hidden units

was the most sensitive to the rarely occurring class 3 (balanced accuracy BA = 0.76), while predictions for class 1 were less accurate (BA = 0.834). Among all the models used in the work, the XGBoost model best predicted ice phenomena from class 1 (BA = 0.923) but demonstrated the weakest prediction of phenomena from classes 2 (BA = 0.954) and 3 (BA = 0.64) (Table 4). The predictive power of the tested models for the Uniejów station was different depending on the class of phenomena (Table 3). At this location, the XGBoost model achieved the highest values of balanced accuracy in the test set for ice phenomena from classes 1 (BA = 0.951) and 2 (BA = 0.984). At the same time, the NN5 model was the only one to correctly predict ice phenomena from class 3 (BA = 0.998).



\* No means no ice phenomena.


**Table 5.** Results of predictive modeling of ice phenomena for the Warta River (middle course).

\* No means no ice phenomena.


**Table 6.** Results of predictive modeling of ice phenomena for the Warta River (lower course).

\* No means no ice phenomena.

In the middle course of the Warta River, for the Nowa Wies water gauge, the NN5 model turned out to be the best at predicting ice phenomena (Table 5). This model performed well for ice phenomena from classes 1 (BA = 0.906) and 2 (BA = 0.949), comparable with other models, while at the same time being the most sensitive for class 3 (BA = 0.623). However, the best performance in predicting class 1 events was achieved by the NN6 model (BA = 0.936), and the best performance for class 2 by the NN3 model (BA = 0.965). The XGBoost model achieved similar performance to the NN5 model as regards the prediction of phenomena from class 2 (BA = 0.945). In the case of the Srem water gauge, it is difficult to indicate the best model (Table 5). A neural network model NN5 best predicted the ice phenomena from class 3 (BA = 0.749). The NN6 model showed the best prediction for class 2 (BA = 0.993), and class 1 events were best predicted by XGBoost (BA = 0.946). Nevertheless, it is the neural network model with five hidden units (NN5) that seems to have the most balanced prediction profile for all classes of ice phenomena. For the Poznan water gauge (Table 5), the predictive power of the tested models was comparable. The NN3 model can be viewed as the best for predicting ice phenomena at this location because it predicted classes 1 (BA = 0.933) and 2 (BA = 0.958) best and was the third most effective in

predicting the level of class 3 events (BA = 0.706). At the same time, the NN5 model turned out to be the best for ice phenomena from class 3 (BA = 0.759).

For the lower course of the Warta River, in the Skwierzyna profile, the predictive power of the tested models was comparable (Table 6). The NN3 model appears to present the most balanced predictive profile. This model best predicted the occurrence of ice phenomena from class 1 (BA = 0.958) and also the absence of river freezing (BA = 0.973). The best prediction of ice phenomena from class 2 was achieved by the NN6 model (BA = 0.982), with the results of the XGBoost model being comparable (BA = 0.98). The NN3 model also displayed good predictability of phenomena from classes 2 (BA = 0.965) and 3 (BA = 0.692). Its class 3 prediction performance is comparable to that of the NN4 model, for which BA = 0.694 was determined. For the Gorzow Wlkp. water gauge, one of the better predictive models for ice phenomena was the NN4 model (BA = 0.932 for class 1, BA = 0.98 for class 2) (Table 6). The NN5 model predicted classes 1 and 2 comparably to the other models, and at the same time was the most sensitive in terms of predicting ice phenomena from class 3 (BA = 0.665). The XGBoost model predicted the phenomena from group 2 best (BA = 0.987), similarly to the NN3 model (BA = 0.983).

#### 4.3.1. Spatial Differences in Model Performance

Among the NN models used, the best predictions were given by the NN5 (eight-fold confirmation of the best prediction) and NN4 models (seven-fold) (Table 7). The XGBoost model also has high predictive power, and the model turned out to be the best in predicting ice phenomena from classes 1 and 2. In three cases, its performance was comparable with those of the NN models. The phenomena from the initial stage of freezing (class 1) were best predicted by the XGBoost model. On the other hand, the disintegration of the ice cover and accompanying ice phenomena were best predicted by the NN5 model (at five water gauge stations). No dependence of the models' performance on the location of water gauges (Table 7) was observed, although as regards predictions of ice phenomena in the upper section of the Warta River (Bobry, Sieradz, Uniejow stations), the XGBoost model and the NN4 and NN5 models proved to be superior.


**Table 7.** Models with the best prediction of ice phenomena on the Warta River.

\* No results—no results from the learned relations between the class of the phenomenon and predictors.

Ice phenomena predictions for the river along its middle section (stations in Nowa Wies, Srem, and Poznan) were made most reliably by the XGBoost and NN5 models (Nowa Wies and Srem) and the NN3–NN5 models (Poznan) (Table 7). For the prediction of ice phenomena along the lower section of the Warta, superior performance was demonstrated by the NN models, taking into account the lower efficiency of the XGBoost model.

The most difficult prediction was that for ice phenomena in the decay phase and the formation of ice floes and, consequently, ice jams. Due to the lowest frequency of observations, there were problems with their prediction in the case of the Bobry station. In this case, no results were obtained from the relations determined between the class of ice phenomena and the predictors.

#### 4.3.2. Evaluation of the Importance of Predictors in the Models

The use of XGBoost, as opposed to ANNs, made it possible to assess the importance of predictors in the model. The selected predictor variables were ranked according to the normalized reduction in model error, also known as "variable importance". Figure 8 shows the most important predictor variables in the final model: water and air temperature, hydrological conditions (water level and river flow), and data for the "day before", month, and year. The results of this analysis indicate that for each measuring station the most important predictor of ice phenomena is the type of ice phenomenon the day before the identification of a given event, with water temperature and air temperature coming next. In the case of the stations in Uniejow and Srem, water temperature is the second most important predictor of the occurrence of ice phenomena.

**Figure 8.** *Cont*.

**Figure 8.** Relative importance of predictors in the XGBoost model (profit values).

These results suggest that when looking for a balance between the complexity of the model and its predictive power, the two most important predictors for the occurrence of ice phenomena on the Warta River should be taken into account, i.e., the nature of the ice phenomenon on the day preceding the observation (especially for class 2 or class 1 events), and water temperature.

#### **5. Discussion**

#### *5.1. Selection of Predictors as Input Variables*

The predictive modeling of ice phenomena carried out on the example of the Warta River showed that the prediction of their occurrence in different phases and spatial locations gives different results. In this case, the prediction was a difficult process, mainly due to the complexity of interactions between hydroclimatic factors and thermal conditions that contribute to the occurrence of freezing.

In the research conducted on the Warta River, an important assumption was the selection of input variables that affect the accuracy of predictions in the neural network models and XGBoost. A set of daily data were used, these including thermal and hydrological variables, the type of ice phenomenon (group of phenomena) on the day preceding water gauge observations, and the month of their occurrence. The premises confirming the correctness of their choice are the results of studies of the ice regimes of rivers in Poland, including those conducted on the Vistula River [13,58], Oder River [59], on the rivers of the Baltic coastal zone [60,61], Bug River [62], and Warta River [46,63,64]. The selection of input variables significantly affects the performance of ANN models [7], however, it is often arbitrary [8].

#### *5.2. The Most Important Predictor Variables in the Final Model*

The results of the predictive models that we developed for the Warta River showed that all the input parameters (predictors) that were taken into account had some significance for the formation of ice phenomena from different classification groups. However, under the thermal conditions established for the reference period (research period 1984–2013), hydrological parameters—river flow and water level—were less important for the process of ice phenomena formation. The research established that ice phenomena occurred irregularly and periodically in the studied period and that the structure of freezing along the river course was diversified. The phenomena from class 1 were predominant, i.e., from the freezing phase of the river, represented by frazil ice and border ice, which is now a typical feature of the ice regime of most rivers in Poland [47].

The results of modeling confirmed that the most important predictors in the analyzed case were the nature of the phenomenon on the day preceding the observation (most often class 2 or class 1), as well as water temperature, and then air temperature (Figure 8). Graf [12] examined the dependencies of the trends of ice phenomena in the Note´c River, in western Poland (a tributary of the Warta River), on air and water temperature using regression models for count data and the Zero-Inflated Negative Binomial Model; results

showed that the temperature values are the best predictors. In some locations, however, the model predicting the number of ice phenomena—taking into account the relationship with temperature—turned out to be statistically insignificant. Graf and Tomczyk [11] determined that for the Note´c River, a faster increase in accumulated sequences of negative air temperature contributes to an increase in the probability of a permanent ice cover, and the average degree day increase by one degree increases the chance of ice cover on the river in the range of 1.5–6.0% in different water gauges.

The period of intense changes in thermal conditions in the Warta catchment area, e.g., the cold period or sudden spring river supply, can be represented by changes in the models. However, it is also visible in the types of distributions illustrating the relationship between the classes of ice phenomena on the Warta River and hydrological factors and thermal conditions, which has been presented in the violin plots. The violin plots show diverse and complicated relations resulting from the differing variability of hydroclimatic factors and thermal conditions, which determines the nature of the distribution (Figure 5). Conditions conducive to the emergence of ice phenomena are not always the same in every location on the river, which is the result of local conditions, including, e.g., channel morphology and the influence of anthropic pressure.

#### *5.3. The Performance of Predictive Models*

The models performed promisingly in predicting the occurrence of ice phenomena on the Warta River, and this—in addition to their low demand for computational data resources, speed of operation, and ease of use—makes them particularly attractive. Further, it was found that the ANN approach served its purpose. By using more advanced and specialized network architectures (NN5, NN6), the ability to learn and predict the non-linear behavior of ice phenomena was increased for classes 2 and 3, which were characterized by a lower frequency of occurrence in the Warta River.

In ANN predictive modeling, the use of the sigmoidal function as an activation function for output layer neurons ensured that predictions obtained with the model would be the probabilities of occurrence of a given class, since this function maps real values to the range 0–1. As a result, models with 3, 4, 5, and 6 neurons in the hidden layer were developed, and this made it possible to compare their performance in predicting ice phenomena. Guo et al. [4], basing their work on the ANN theory, used the sigmoidal function as the activation function in the hidden layer for the forecasting of ice jams during river ice breakup. Their results were promising, as they predicted the annual occurrence of ice blockages with an accuracy of 85%, while the projected decay date with the projected ten-day period showed a maximum error of two days.

Concerning the prediction of class 1 and class 2 phenomena in the Warta River (permanent ice cover) and their non-occurrence, ANNs require further improvement, although present results indicate that they are comparable to the XGboost algorithm for predicting group 2 phenomena, i.e., permanent ice cover. The performance of the NN models and the XGBoost algorithm is also comparable for the different water gauge locations on the Warta River, although an overall better fit of XGBoost and NN4/NN5 models was observed for the upper course of the river; XGBoost and NN5/NN3 were most successful for the middle course, while NN models predominated for the lower course. The results of a comparison of both types of models in terms of their suitability for predicting ice phenomena on the Warta River showed a high accuracy of prediction for the XGBoost method, which has not been used on a larger scale in this regard so far.

XGBoost models variable interactions and handles the multi-linearity common to ecological datasets seamlessly [65]. Moreover, XGBoost works faster than many other gradient-increasing algorithms due to the regularization factor and the parallel computing functionality. One of the advantages of this method is its resistance to outliers, which eliminates the need to supplement missing data, and thus, in the case of the Warta River, it allows an increase in the efficiency of the prediction of a given phenomenon, even

when eliminating outliers. The obtained results were considered satisfactory, which was confirmed by four model fit measures.

The comparison of results obtained for the Warta River with the ANN models used to predict ice phenomena on other rivers shows their considerable similarity. Most of the predictive and prognostic models developed confirm that the results of ice condition forecasts made with the use of ANNs are satisfactory and consistent with the measured data [37,66,67]. Moreover, the high accuracy of forecasts is indicated, which takes into account factors influencing the formation and disappearance of ice phenomena. Too many simplifications, as made in some models, may lower their prognostic accuracy and limit their usefulness for other rivers [36]. According to Massie [33], neural network classifiers, just like in the case of empirical methods, are most likely location-specific, but it is possible to transfer ANN models to other locations with minimal modifications. However, there are still no solutions for the prediction of phenomena in individual phases of their occurrence using the XGBoost algorithm. In Poland to date, models from the ANN group and the XGBoost algorithm have not been used to predict ice phenomena.

A review of the literature shows that numerous parameters are needed to support models developed for forecasting ice phenomena, most commonly ice jams and the resultant floods, but obtaining this data is sometimes difficult or even impossible. Despite the progress made in forecasting ice processes on rivers, this field still has great research potential; however, it also requires comprehensive observations, the collection and testing of data from stationary measurement networks, and direct field studies [8].

#### **6. Conclusions**

In the present study, MLPNN and XGBoost models were developed to forecast ice phenomena on the Warta River in Poland. The results obtained lead to the following conclusions:


The results of the research conducted here have important implications for forecasting ice phenomena, specifically as regards the application of XGBoost. Preliminary results seem to indicate that XGBoost, as an ensemble machine learning model, works well as a forecasting tool in hydrological research. Though the MLPNN and XGBoost models performed competently, there is still scope for further improvement through additional studies and the construction of hybrid models. Other factors influencing the occurrence of ice phenomena on rivers that would additionally help to improve the accuracy of these models should also be looked at (e.g., channel morphology, the accumulated degree days of frost and thaw, and the rates of change in water level and flow during the freeze-up

and breakup periods). Since the present results concern only one river, future research will focus on applying models to rivers in different geographic locations and hydrological regimes to more accurately test the suitability and effectiveness of models.

**Author Contributions:** Conceptualization, R.G.; funding acquisition, R.G.; methodology, R.G.; project, R.G.; administration, R.G.; investigation, R.G.; software, R.G.; formal analysis, R.G.; validation, R.G.; visualization, R.G.; writing—original draft preparation, R.G; writing—review and editing, R.G., T.K., and S.Z.; supervision, R.G., T.K., and S.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding from an internal grant "Research subsidy for scientific research activities of the Faculty of Geographical and Geological Sciences of the Adam Mickiewicz University in Poznan, Poland.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Hydrological and meteorological measurement data are available directly from https://hydro.imgw.pl (accessed on 10 July 2021) and https://danepubliczne.imgw.pl (accessed on 15 July 2021) as operational data. The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** The authors thank the Faculty of Geographical and Geological Sciences, Adam Mickiewicz University, Poznan, Poland for their support. The authors would also like to acknowledge the Institute of Meteorology and Water Management—National Research Institute (IMWM-NRI, Warsaw, Poland) for the release of the input database.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Changes in Terrestrial Evaporation across Poland over the Past Four Decades Dominated by Increases in Summer Months**

**Urszula Somorowska**

Department of Hydrology, Faculty of Geography and Regional Studies, University of Warsaw, Krakowskie Przedmie´scie 30, 00-927 Warsaw, Poland; usomorow@uw.edu.pl

**Abstract:** Given the importance of terrestrial evaporation (ET) for the water cycle, a fundamental understanding of the water quantity involved in this process is required. As recent observations reveal a widespread ET intensification across the world, it is important to evaluate regional ET variability. The specific objectives of this study are the following: (1) to assess annual and monthly ET trends across Poland, and (2) to reveal seasons and regions with significant ET changes. This study uses the ET estimates acquired from the Global Land Evaporation Amsterdam Model (GLEAM) dataset allowing for multi-year analysis (1980–2020). The Mann–Kendall test and the Sen's slope were applied to estimate the significance and magnitude of the trends. The results show that a rising temperature, along with small precipitation increase, led to the accelerated ET of 1.36 mm/y. This was revealed by increased transpiration and interception loss not compensated by a decrease in bare soil evaporation and sublimation. The wide-spread higher water consumption especially occurred during the summer months of June, July, and August. Comparing the two subperiods of 1980–2020, it was found that in 2007–2020, the annual ET increased by 7% compared to the reference period of 1980–2006. These results can serve as an important reference for formulating a water resources management strategy in Poland.

**Keywords:** terrestrial evaporation; components; Global Land Evaporation Amsterdam Model (GLEAM); increasing trends; spatial–temporal pattern; evaporative ratio; Poland

#### **1. Introduction**

Terrestrial evaporation (ET), alternatively called land surface evaporation or evapotranspiration [1], is an important component of the global water cycle. It consists of biophysical (transpiration from vegetation) and physical (evaporation from the interception, bare soil, and open water) water fluxes. The contribution of these different fluxes to the total amount of evaporated water depends on both the climate controlling the atmospheric water demand, and the land surface features, especially vegetation characteristics, influencing the energy balance at the land surface and determining the volume of interception [2]. At a global scale, the largest share in terrestrial evaporation is transpiration (Et; 59%), followed by vegetation and floor interception (Ei; 31%), soil moisture evaporation (Es; 6%), and, lastly, open water evaporation (Ew; 4%) [3]. As changes in terrestrial evaporation can lead to either warming or cooling of the land surface [4], acquiring a better understanding of ET changes is of high priority for ongoing research. Additionally, knowledge of the temporal changes of ET is necessary for accurately quantifying global and regional water budgets and for a better understanding of the hydrological interactions between the land and the atmosphere [5]. Revealing the current ET trends might contribute to the prediction of the Earth's runoff changes, and to the introduction of water management strategies.

Globally, ET has shown a significant upward trend over the 32-year period of 1982–2013, mainly driven by vegetation greening and rising atmosphere moisture demand [6]. It has been mainly caused by increases in transpiration from vegetation and the vaporization of intercepted rainfall from vegetation [7]. An increase in global terrestrial ET was also

**Citation:** Somorowska, U. Changes in Terrestrial Evaporation across Poland over the Past Four Decades Dominated by Increases in Summer Months. *Resources* **2022**, *11*, 6. https://doi.org/10.3390/ resources11010006

Academic Editor: Diego Copetti

Received: 7 December 2021 Accepted: 9 January 2022 Published: 12 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

reported for the period 1982–2011, as estimated by different remote sensing-based physical models, machine-learning algorithms, and land surface models [8]. Despite the overall increasing trend, decreasing ET trends were also observed in shorter periods of time. For example, it was reported that global annual evapotranspiration declined due to the limited moisture supply from 1998 to 2008 [9]. Additionally, a long-term (1950–2017) relative decline in evapotranspiration, accompanied by increasing runoff, was found in 27% of the global land areas, which was explained by a reduction in surface conductance [10]. Apart from long-term changes observed globally, the strongly regional and temporal differentiation of ET trends are evidenced, as proven for central Europe and central North America [11], and for the north-eastern United States [12]. For the period 1980–2010, largescale increases in ET were observed in south-eastern China, while decreases in ET occurred over the north-east [13]. Other studies conducted over the 71-year period from 1948 to 2018 confirmed that ET exhibited an increasing trend over almost 90% of the territory of China [14].

In most of Europe, ET increased in response to land use (mainly large-scale reforestation and afforestation) and climate change, except in the Iberian Peninsula and some other parts of the Mediterranean where negative ET trends were found [15]. For example, Great Britain is becoming warmer and wetter, through which increases in precipitation and air temperature are driving increases in runoff and evapotranspiration, as proven for 1961–2015 [16]. In Germany, significant trends were observed in transpiration and evapotranspiration in the period of 1961–2019 [17]. Over European Russia, positive trends in the annual values of potential evapotranspiration were revealed for 1966–2017, the distribution of which has a strong, patchy character [18]. As projected for most of Europe, widespread and relatively uniform ET increases of around 75–125 mm are expected by the 2050s [19]. Overall, the trends in evapotranspiration exhibit wide regional differentiation, depending on varying climate and human impacts. The changes are linked to the main climatic drivers of ET that comprise precipitation, net solar radiation, air temperature, vapor pressure deficit, and wind speed. Transpiration, which is the main component of ET in forested landscapes [7,20], is dependent not only on the climate conditions, but also on the soil moisture regime, the physiological features of vegetation, and the duration of the growing season [21,22].

In Poland, only a few studies have focused on changes in ET from a long-term perspective. Based on data from 18 weather stations across Poland, an increase in the growing season reference evapotranspiration was detected in the period from 1971 to 2010 [23]. Moreover, a significant increase in the mean daily values of the summer reference ET in Poland in the same multiyear period was observed, with the potential to occur more frequently in the future [24]. Projections of changes of areal ET in the Wielkopolska Region (central-western part of Poland) indicate that the regional average increase in the annual ET is predicted to equal 45 mm, with the highest changes occurring during the winter, when comparing the control period of 1961–1990 to the projection horizon of 2061–2090 [25]. In the north-western part of Poland, significant increasing trends in the potential ET were observed during the period of 1952–2018, as investigated based on data acquired from six weather stations [26]. In this area, particularly high increases were detected during the spring and summer months, explained by a significant increase in air temperature and a decrease in relative humidity, marking the end of the 20th century and the start of the 21st century.

While the above-mentioned studies investigated the actual, potential, or reference ET trends from selected weather stations' data across Poland, none identified countryscale spatial–temporal patterns of ET derived from the satellite-based or reanalyzed data covering the whole territory of Poland. Hence, undertaking such research can bridge the gap and improve the understanding of the ongoing ET changes. This might provide insight into the hydrological implications that result from the ET variations helping to manage the water resources at the country and regional scales. This paper uses well-validated, satellite-based Global Land Evaporation Amsterdam Model (GLEAM) data that comprise

spatially continuous estimates of terrestrial evaporation and its main components [27–29]. For the first time, ET is quantified throughout the recent 41-year period (1980–2020), and the spatial–temporal trends are estimated across the country. New insights into the ET variability and its changes are acquired by considering the grid cell and country scale patterns. Thus far, such an approach has not been given detailed consideration with respect to the territory of Poland. Apart from investigations focusing on specific regions and based on point estimates of ET, quantitative multi-year studies incorporating spatial–temporal characteristics have not been reported.

The key scientific questions that need to be answered in this study are the following:


It is hypothesized that the annual and monthly ET has remarkably increased in the past four decades, and that extremely dry and warm years have led to a high evaporation ratio, defined as the ratio between the actual terrestrial evaporation and precipitation. It is also hypothesized that the ET trends are regionally differentiated, with the highest trend rates occurring in the summer months.

#### **2. Materials and Methods**

#### *2.1. Basic Geographical Characteristics of Poland*

Poland is located in the central part of Europe. The area of the country stretches from 54.83◦ N to 49.00◦ N and from 14.18◦ W to 24.15◦ W. The land area of the country covers an area of 311,888 km2. While the average elevation is approximately 170 m a.s.l., the relief is characterized by significant diversity. The lowest point is located in Zuławy Wi´ ˙ slane (−1.8 m a.s.l.) at the Baltic coast in northern Poland, and the highest peaks are in the Tatras Mountains (Rysy 2499 m a.s.l.) in the south. Lowlands (0–300 m a.s.l.) dominate and cover 91% of the country's area. Highlands (300–500 m above sea level) occupy 6%, while mountains (above 500 m above sea level) occupy only 3%. As much as 99.7% of the country's area constitutes the drainage basin of the Baltic Sea, which consists of the Vistula River Basin (53.7%), Oder River Basin (33.9%), and the river basins draining directly to the Baltic Sea (12.1%). A large part of the country (from the 14th to approximately 20–22 meridians) is influenced by a humid, temperate (Cfb) climate [30]. The rest of the area to the east and the mountainous parts to the south are classified as a snowy, humid climate (Dfb) with a warm summer. The highest annual precipitation total of approximately 1600 mm/y occurs in the Tatra Mountains, while the lowest, below 500 mm/y, concern the central, lowland part of the country. The average air temperature in the months from April to October ranges from 9.3 ◦C in the mountains to 14.9 ◦C in the midwest lowland, with a country-wide average of 13.8 ◦C. The land use is dominated by agricultural land (60%) and forests (33%), while the rest of the territory is covered by artificial surfaces (6%) and other categories (1%) (Figure 1).

**Figure 1.** (**a**) Geographic location of the study area. (**b**) Land use according to CLC 2012 acquired from https://land.copernicus.eu, accessed on 7 June 2021.

#### *2.2. Data*

In this study, the satellite-based GLEAM (https://www.gleam.eu, accessed on 20 March 2021) dataset was used for detecting the magnitude and trends in terrestrial evaporation for 1980–2020. Driven by remote sensing data, GLEAM provides estimates of different components of land evaporation at a resolution of 0.25◦ × 0.25◦, including transpiration (Et), vaporization of intercepted rainfall from vegetation (Ei), evaporation from the bare soil (Eb), snow sublimation (Es), and open-water evaporation (Ew). Version 3.5a is derived from the reanalysis of net radiation, air temperature, satellite and gauge-based precipitation, satellite-based vegetation optical depth (VOD), ESA-CCI (European Space Agency Climate Change Initiative) soil moisture, and snow water equivalent [27,28]. Changes in land cover are derived from the MEaSUREs Vegetation Continuous Fields dataset [31]. In this study, the monthly values of terrestrial evaporation were acquired and extracted for the territory of Poland. Moreover, the annual values of evaporation components (Et, Ei, Eb, Es, and Ew) were also acquired within the same spatial domain consisting 1056 grid cells. It is worth mentioning that spatial–temporal ET patterns might be captured by different regional and global ET products [32]. The GLEAM ET dataset belongs to the sophisticated land surface model group, which has an acceptable accuracy when compared to the benchmark ET values [33]. It has been recently used as a diagnostic dataset to investigate the global ET trends [34], to identify best-performing evaporation datasets [35], and to improve the structure of a simple conceptual rainfall–runoff model [36]. The key distinguishing features of these data are the multi-year observation period (1980–2020), the availability of five pieces of terrestrial evaporation component data, and the long-term satellite surface soil moisture and phenology observations assimilated into the GLEAM model.

For the analysis of the climate conditions, the air temperature and precipitation dataset was employed for 1980–2020 from version 23.1e of the station-based E-OBS gridded dataset (https://www.ecad.eu, accessed on 9 October 2021), available from the European Climate Assessment and Dataset Project [37]. It comprised daily values of precipitation (P) and air temperature (T) acquired from a regular latitude/longitude grid of 0.25◦ × 0.25◦. From the daily gridded values, the country scale, annual air temperature, and precipitation estimates were calculated. Due to E-OBS uncertainties caused by gauge measurement errors, lack of wind corrections, and possible station relocations [38], the areal precipitation averages are considered to be underestimated [39,40]. To overcome the underestimation of the precipitation data available, a correction factor was applied to the annual precipitation totals derived from the E-OBS data. The correction factor results from the linear relationship between the E-OBS data (PE-OBS) and the G2DC-PL+ data (PG2DC+), a gridded 2 km daily climate dataset for the Polish territory, which is corrected for snowfall and rainfall under-catch and spans until 2019 [41,42]. The relationship is represented by the equation PG2DC+ = 1.2 × PE-OBS

established for the period 1980–2019. The value of the Pearson's correlation coefficient is 0.996, which confirms the strong linear relationship between the two variables. Thus, the PE-OBS-corrected value of the corrected E-OBS precipitation is estimated from the following equation: PE-OBS-corrected = 1.2 × PE-OBS, where annual values of PE-OBS are acquired for the period 1980–2020. The corrected values of annual precipitation totals were used to calculate the values of the evaporation ratio, defined as a ratio of the actual evaporation over a given area to the precipitation falling on that area [43]. Then, the annual precipitation and air temperature time series were tested for the presence of trend, as explained in Section 2.3.

#### *2.3. Methods*

The non-parametric Mann–Kendall (MK) test [44,45] was applied for detecting trends in time series of terrestrial evaporation, precipitation, and air temperature data. It is widely applied test to analyze the hydro-meteorological time series [46]. Here, it was used to test the null hypothesis of no trend, H0, against the alternative hypothesis, H1, where there is an increasing or decreasing monotonic trend. The test statistic *S* is defined as

$$S = \sum\_{k=1}^{n-1} \text{sgn}(\mathbf{x}\_{\bar{j}} - \mathbf{x}\_k) \tag{1}$$

where *xj* and *xk* are the annual values in years *j* and *k*, *j* > *k*, respectively, n is the length of the time series, and

$$\text{sgn}(\mathbf{x}\_{j} - \mathbf{x}\_{k}) = \begin{cases} 1 & \text{if } \begin{pmatrix} \mathbf{x}\_{j} - \mathbf{x}\_{k} \end{pmatrix} > 0 \\ 0 & \text{if } \begin{pmatrix} \mathbf{x}\_{j} - \mathbf{x}\_{k} \end{pmatrix} = 0 \\ -1 & \text{if } \begin{pmatrix} \mathbf{x}\_{j} - \mathbf{x}\_{k} \end{pmatrix} < 0 \end{cases} \tag{2}$$

The *S* statistic is approximately normally distributed when *n* ≥ 8, with a mean of 0 and the variance of statistics *S*, *σ*2, given as

$$
\sigma^2 = \left[ n(n-1)(2n+5) \right] / 18 \tag{3}
$$

The standardized test statistic *Z* is computed by

$$Z = \begin{cases} \begin{array}{l} (S-1)/\sigma \end{array} & \text{if } S > 0\\ 0 & \text{if } S = 0\\ (S+1)/\sigma & \text{if } S < 0 \end{array} \tag{4}$$

For the country scale, annual time series of precipitation, air temperature, terrestrial evaporation, and its components, the trends were tested at a significance level of *α* = 0.05. For *α* = 0.05, the null hypothesis is accepted when −1.960 ≤ *Z* ≤ 1.960 (no significant trend), while it is rejected when *Z* < −1.960 (significant decreasing trend) or when *Z* > 1.960 (significant increasing trend). Thus, the trend is significant if the null hypothesis cannot be accepted. For the monthly terrestrial evaporation time series in a multi-year period 1980–2020, the significance level *α* = 0.1 was assumed. For the assumed *α* = 0.1, the null hypothesis is accepted when −1.645 ≤ *Z* ≤ 1.645 (no significant trend), while it is rejected when *Z* < −1.645 (significant decreasing trend) or when *Z* > 1.645 (significant increasing trend). For the computation of the magnitude (slope) of an existing trend (as change per year), the directional coefficient *β* expressed by the Theil–Sen estimator [47,48] is calculated by the formula:

$$\beta = \operatorname{Median}\left(\left(\mathbf{x}\_{j} - \mathbf{x}\_{k}\right) / \left(j - k\right)\right) \tag{5}$$

A positive value of *β* signals an increasing trend, and a negative value of *β* indicates a decreasing trend. In case the change is not statistically significant, but shows an inclination, it is called a tendency. The Climate Data Toolbox (CDT) for MATLAB [49] was used to calculate the MK test statistic (*Z*) and the *p*-values.

For the shift-type changes, the terrestrial evaporation time series were partitioned into two subseries by minimizing the sum of the residual (squared) error of each subset from its local mean, and finally returning the index, which, in this case, is the year in which the change occurs. The change point was identified for the country-scale time series of annual ET sums by using the MATLAB function "findchangepts". In the analysis performed here, a complete 41-element ET time series was divided into subperiods of differing length, consisting of 27 and 14 records covering the years 1980–2006 and 2007–2020, respectively. Finally, the Kruskal–Wallis test was applied to test for statistically significant differences between the subseries in two selected subperiods [50,51] by using HYDROSPECT software [52]. This test is a non-parametric test that compares mean ranks (i.e., medians). For this test, the null hypothesis is that the subseries medians are equal, versus the alternative that there is a difference between them. Under the null hypothesis of equal subperiod means, the statistic follows the Chi-square distribution. The Kruskal–Wallis test statistic is calculated as

$$H = \frac{12}{n(n+1) - T} \sum\_{i=1}^{N} n\_i (m\_i - m)^2 \tag{6}$$

$$T = \frac{1}{n-1} \sum t\_j^3 - t\_j \tag{7}$$

where *n* is the series length, *m* the global mean of ranks, *N* denotes the number of subperiods, *ni* is the number of values in the *i*-th subperiod, mi is the mean rank for the *i*-th subperiod, *T* is the "tie correction", and *tj* denotes the number of ties in subsequent tie groups [50]. The value of *T* for *t* = 1 is equal to 0. The Kruskall–Wallis test was also applied to the 41-element monthly ET time series partitioned into two 27- and 14-element subseries, as a split series of annual ET totals.

#### **3. Results**

#### *3.1. Country Scale, Inter-Annual Terrestrial Evaporation over the Period 1980–2020*

Figure 2 shows the course of the country scale of the annual sum of terrestrial evaporation comprised of its components in the analyzed period 1980–2020. It is within the range of 409–491 mm, and the multi-year average (1980–2020) equals 455 mm. The largest share in ET is Et (78%), and the second quantitatively important component is Ei (17%). The other contributions, Eb, Es, and Ew, have much smaller shares, i.e., 2%, 1%, and 2%, respectively. The annual ET shows an increasing trend of 1.36 mm/y in response to the slightly increasing tendency in precipitation (P) and the significant increase in air temperature (Figure 3a). The increasing ET is the result of the statistically significant increases in Et, Ei, and Ew, which outweigh the decreasing trend in Eb and the decreasing tendency in Es (Figure 3a,b). Noteworthy are the high values of Es in 1995 and 2006, which presumably were caused by particularly long-lasting snow cover, usually occurring from mid-December to the end of April. The multi-year rate of change of Et, Ei, Eb, Es, and Ew is equal to 0.84, 0.61, −0.04, −0.08, and 0.03 mm/y, respectively.

A shift-type temporal change of annual ET occurred in 2007, as detected by the change point analysis (Figure 3c). The Kruskal–Wallis test confirmed the significance of the difference between the median values of ET in the two subperiods. It shifted from 449 mm in the subperiod of 1980–2006 to 480 mm in the subperiod of 2007–2020. Additionally, the mean annual ET changed from 444 mm to 476 mm. There was also a significant shifttype change in the components of terrestrial evaporation, and this concerns transpiration, interception loss, and open-water evaporation. The average evaporation ratio (ET/P) in the period before the change was 65%, while after the change, it increased by 3% to a value of 68%. However, in particularly wet years, it was much lower, while in dry years, it reached very high values. For example, in the wet years of 2010 and 2017, it was equal to 50% and 56%, respectively, while in dry years, it was far above the long-term average, reaching 95% in 1982, 84% in 2015, and 87% in 2018.

**Figure 2.** Annual terrestrial evaporation (ET) over Poland in 1980–2020, composed of transpiration (Et), interception loss (Ei), open-water evaporation (Ew), snow sublimation (Es), and bare soil evaporation (Eb).

Thus, atmospheric conditions, represented here by temperature and precipitation, control extreme ET values affecting the range of anomalies. In this work, precipitation and temperature are considered as selected ET drivers, although it should be noted that other factors, such as wind, solar radiation, and air humidity, can also play an important role. Except for climatological–meteorological conditions, other features such as hydrogeological, topographical, and physiological characteristics can also be critical.

#### *3.2. Country-Scale, Monthly Terrestrial Evaporation over the Multi-Year Period 1980–2020*

Figure 4 shows the course of the country-scale, monthly ET in the analyzed period 1980–2020. Statistically significant increasing trends occurred in summer and early fall, from June to October. In winter, an increase took place in January (0.08 mm/y). In the remaining months, the changes were not significant. The highest increases occurred in June (0.37 mm/y), July (0.30 mm/y), and August (0.21 mm/y). Much smaller changes concerned the months of September (0.13 mm/y) and October (0.05 mm/y).

An annual cycle is exhibited by ET, with the highest values in the summer season (June–August) and the lowest values in the winter season (November–February), as shown in Figure 5. Comparing the two subperiods, the second subperiod (2007–2020) is marked by a clear increase in ET in June (from 73 mm to 81 mm), July (from 77 mm to 84 mm), and August (from 63 mm to 68 mm) (Figure 5b). In these months, the increase in ET reached on average of 9%. The shift type changes in ET were statistically significant at 6 = 0.05 in January and in all months from June to September, as determined by the Kruskal–Wallis test. The same is true for ET totals in the months from April to October.

#### *3.3. Monthly and Annual Spatial Patterns of Terrestrial Evaporation*

Figure 6 shows the mean monthly spatial patterns of ET for 1980–2020. In December, January, and February, the values were the lowest. From March on, throughout April and May, they gradually increased, reaching the highest values in June, July, and August. Then, starting from September, they gradually decreased, reaching the lowest values in the winter months. Throughout the whole year, the highest monthly ET occurred in the south, and it concerned the mountainous, forested area. The area with relatively high ET values was also located in the western and northwestern part of Poland, where there were dense forest complexes. The lowest monthly ET were characteristic for the central zone, stretching from west to east.

**Figure 3.** (**a**) Changes in annual precipitation (P), air temperature (T), terrestrial evaporation (ET), and transpiration (Et) in the multi-year period 1980–2020. (**b**) Changes in evaporation from interception (Ei), evaporation from bare soil (Eb), sublimation (Es), and evaporation from water (Ew) in the multiyear period 1980–2020. (**c**) Terrestrial evaporation (ET) in two subperiods (1980–2006 and 2007–2020), separated by a change point in 2007. The presence of a trend is determined at the significance level of *α* = 0.05.

To check if there was a trend in monthly ET, 41-element time series were prepared for each grid cell for each month. Then, the Mann–Kendall test was applied. The results are presented in Figure 7 and Table 1. A statistically significant increasing trend occurred in many grid cells in June–September, covering 96% of the country's territory in June (Figure 7f), 89% in July (Figure 7g), 84% in August, and 75% in September. An increasing ET trend was also observed in 59% of the territory of Poland in January (Figure 7a).

**Figure 4.** Monthly time series of terrestrial evaporation (ET) over Poland in 1980–2020. Statistically significant trends are indicated by filled markers. The presence of a trend is determined at a significance level of *α* = 0.05.

**Figure 5.** (**a**) Monthly terrestrial evaporation (ET). (**b**) Mean monthly terrestrial evaporation (ET) for subperiod 1 (1980–2006) and subperiod 2 (2007–2020). The ET values are averaged over the territory of Poland.

The trend rate of ET over the 41-year period was the highest in June, July, and August (Figure 8). In June, 63% of the territory experienced a trend rate of 0.20–0.39 mm/y, 29% had a trend rate of 0.40–0.59 mm/y, 2% was marked by a highest trend rate of 0.60–0.79 mm/y, and 2% was marked by a rate of 0.01–0.19 mm/y (Figure 8f). In July, the trend rate of 0.20–0.39 mm/y dominated as well, covering 57% of the territory, while the rate of 0.01–0.19 mm/y (either trend or tendency) was observed in 22% of the territory, a trend rate of 0.40–0.59 mm/y concerned 17% of the territory, and 0.60–0.79 mm/y concerned 2% of the territory. The remaining 2% of the territory experienced a slight decreasing tendency of −0.19–0 mm/y. The contribution of particular trend rate classes changed in August; 46% of the country had a trend or tendency rate of 0.01–0.19 mm/y, while the rates of 0.20–0.39 mm/y, 0.40–0.59 mm/y, and −0.19–0.00 mm/y concerned 43%, 9%, and 2% of the country, respectively. While the rising ET trend dominated in the summer, a slight downward trend or tendency was observed across the country during the spring months.

This applied to the months of March, April, and May when the changes were in the range of −0.19–0 mm/y, and they covered 68%, 48%, and 42% of the country territory, respectively.

**Figure 6.** Spatial patterns of mean monthly ET during 1980–2020 for: (**a**) January, (**b**) February, (**c**) March, (**d**) April, (**e**) May, (**f**) June, (**g**) July, (**h**) August, (**i**) September, (**j**) October, (**k**) November, and (**l**) December.

**Figure 7.** Mann–Kendall test results of ET trend detection over Poland at significance levels of 6 0.001, 0.01, 0.05, and 0.1 in the period 1980–2020 for: (**a**) January, (**b**) February, (**c**) March, (**d**) April, (**e**) May, (**f**) June, (**g**) July, (**h**) August, (**i**) September, (**j**) October, (**k**) November, and (**l**) December. Red and orange grid cells indicate a decreasing trend, light green and blue grid cells show an increasing trend, whereas yellow and yellowish green indicate no trend detected. Four grid cells with no data are shown in white.

**Table 1.** Percent of the country area with decreasing or increasing ET trends, or without changes in ET detected. The presence of a trend is determined at a significance level of *α* = 0.1.



**Table 1.** *Cont.*

**Figure 8.** Trend rate of monthly terrestrial evaporation (ET) in the multi-year period 1980–2020 for: (**a**) January, (**b**) February, (**c**) March, (**d**) April, (**e**) May, (**f**) June, (**g**) July, (**h**) August, (**i**) September, (**j**) October, (**k**) November, and (**l**) December.

Finally, the spatial patterns of the annual ET were elaborated (Figure 9). Overall, the highest ET values occurred in the southern mountainous part of the country (Figure 9a),

determined by the relatively high precipitation totals and the land use dominated by forests. The land strips in the north-west and west were also marked by relatively high ET values, associated with the presence of dense forest complexes in this area. Much lower ET occurred in the central part of the country, but also in the south-western, north-eastern, and eastern parts of the country, where agricultural fields have a large share in the land use. The comparison of ET in the two selected subperiods showed that subperiod 2 (2007–2020) had significantly higher ET values than subperiod 1 (1980–2006) (Figure 9b,c). The minimum ET values observed across the country were 35 mm higher, and the maximum values increased by 63 mm. The annual ET difference in the range of 0.1–20 mm/y was detected in the zone spreading from the south-west to the north-east, while much higher differences in the south reached, locally, up to 115 mm (Figure 9e). Only five grid cells were characterized by an ET decrease of up to 3% (Figure 9f), of which only one cell had a statistically significant decrease (Figure 9d). A total of 87% of the country experienced an increase of 1–10% ET, while a 11–22% increase was characteristic for the rest of the area (Figure 9f). Overall, the statistically significant increase in annual ET concerned 90% of the country's territory (Figure 9d).

**Figure 9.** (**a**) Annual terrestrial evaporation (ET) in the multi-year period 1980–2020. (**b**) Annual terrestrial evaporation (ET) in the subperiod 1980–2006. (**c**) Annual terrestrial evaporation (ET) in the subperiod 2007–2020. (**d**) Trend rate of annual ET in the multi-year period 1980–2020. (**e**) Anomalies in annual ET detected as a difference between ET (2007–2020) and ET (1980-2006). (**f**) Anomalies in annual ET detected as a ratio of ET (2007–2020) to ET (1980–2006).

#### **4. Discussion**

This study demonstrated that terrestrial evaporation has increased both at the annual and monthly time scale. This finding supports the hypothesis that terrestrial evaporation has remarkably increased in the past four decades. The ET increase was caused by increases in transpiration, interception loss, and open-water evaporation, partially counteracted by bare soil evaporation and sublimation decreases. The increasing trends in transpiration (0.87 mm/y) and interception loss (0.68 mm/y) have the largest contribution to the ET

changes. They are presumably driven by increases in the vegetation leaf area index, dominated by greening [6]. In Poland, greening has recently been confirmed by increasing NDVI trends, showing that 59% of the country has been greening continuously since the 1980s and that it has a widespread character [53]. This process affected over 44% of 1980s agricultural land and 87% of 1980s non-agricultural land. In total, this concerned nearly 60% of the country's area. Thus, the intensified greening all over the country coincides with the increase in terrestrial evaporation detected in this study; 81% of the territory is marked by an ET increase of 0.1–2 mm/y (trend or tendency), while 17% has a change of 2.1–4.9 mm/y. Only 2% of the territory has a decreasing tendency. To conclude, most of Poland has experienced an increasing ET trend in the last four decades, and the evaporation ratio was extraordinarily high in dry and warm years, such as 1982, 2015, and 2018.

The follow-up hypothesis is also confirmed; ET trends are regionally differentiated, with the highest multi-year trend rates occurring in June, July, and August, equaling, on average, 0.35 mm/y, 0.30 mm/y, and 0.22 mm/y, respectively. These are signals of an accelerated, more intense water cycle. The main implication of this finding is that increasing ET rates might have a great impact on the other components of the water cycle. If the observed ET trends continue, the associated changes in water budget components might present challenges for water managers at country and regional scales.

In this study, a well-validated satellite-based GLEAM dataset was used for analyzing ET dynamics for 1980–2020. To the best of the author's knowledge, this is the first study for the Polish territory that gives country-scale estimates of terrestrial evaporation over the last four decades. This enabled the analysis of total terrestrial evaporation and its components in the multi-year period. This study explained the relative contribution of each ET component, revealing the dominant role of transpiration and interception loss in the terrestrial evaporation. The mean annual ET of 455 mm/y is comparable to the ET of 442 mm/y, which was estimated for the whole of Europe [7]. However, the contribution of transpiration to the ET differs (69% for Europe and 78% for Poland). Additionally, clear differences concern the role of interception loss and bare soil evaporation. In this study, interception loss (77 mm/y) contributed in the amount of 17% to the ET (455 mm/y), which is almost double the amount estimated for Europe (9%). In turn, the bare soil evaporation contribution was found to be only 2%, while for Europe, it was found to be 21%. Such divergent results might be due to differences in model-dependent partitioning approaches, datasets forcing the ET models, the differences in climatic and land surface characteristics, and vegetation morphological attributes. The implication is that the relationships between particular components of ET should be further explored in a future study. It is worth noting that the multi-year ET average (455 mm) determined in this study is relatively consistent with earlier estimates for the territory of Poland; it is only slightly higher than estimates for the Oder River basin (454 mm) and for the Vistula River basin (432 mm) reported in 1950s [54]. The discrepancy is very low and may be due to the difference in source data (derived from the water budget) and a different period (1920–1940). Much higher ET values (553 mm) were signaled for 1961–1990 in the Wielkopolska region located in central and western Poland [25]. Such contrasting results are probably due to differences in applied methods and different sets of variables forcing the ET models from which areal ET values are derived. Moreover, this region cannot be treated as representative for the whole territory of Poland.

#### **5. Conclusions**

To conclude, this study evaluated the magnitude and multi-year trends in terrestrial evaporation and its components across Poland. Benefitting from the novel satellite observations assimilated to GLEAM, together with the reanalysis data used as the model forcing, it was found that most of Poland experienced significant terrestrial evaporation increases from 1980 to 2020. The warmer climate, along with a small precipitation increase, led to increased vegetation activity. This was revealed by increased transpiration and interception loss not compensated by a decrease in bare soil evaporation and sublimation. The enhanced vegetation activity was manifested by the wide-spread higher water consumption, especially during the summer months of June, July, and August. These monthly increases contributed to the annual changes in terrestrial evaporation of 1.36 mm/y. Comparing the two subperiods of 1980–2020, it was found that in the subperiod 2007–2020, the annual evaporation increased by 7% compared to the reference subperiod of 1980–2006. Further study on the quantification of terrestrial evaporation is required to explain its impacts on the changing water budget structure. This can serve as a reference for formulating a water resources management strategy in Poland.

**Funding:** This research was funded by the Faculty of Geography and Regional Studies, University of Warsaw, Poland, grant no. SWIB55/2021, and by the University of Warsaw, Poland, "Excellence Initiative–Research University" program, Action IV.3.1, grant no. BOB-661-519/2021.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the data used during this study are available at locations cited in the acknowledgements section.

**Acknowledgments:** The author thanks the editors and reviewers for their insightful comments and suggestions. The author acknowledges the GLEAM dataset (version 3.5a), which is available at https://www.gleam.eu (accessed on 20 March 2021). The author also acknowledges the E-OBS dataset from the ECA&D project (http://www.ecad.eu, accessed on 9 October 2021); version 23.1e of the dataset was downloaded from the Copernicus Climate Change Service (https://surfobs. climate.copernicus.eu/dataaccess/access\_eobs.php) (accessed on 9 October 2021). The CORINE Land Cover was downloaded from the Copernicus Land Monitoring Service (https://land.copernicus.eu), European Union, European Environment Agency (EEA) (accessed on 7 June 2021). The author also acknowledges the G2DC-PL+ dataset openly available from the 4TU Centre for Research Data (https: //doi.org/10.4121/uuid:a3bed3b8-e22a-4b68-8d75-7b87109c9feb, accessed on 16 September 2021).

**Conflicts of Interest:** The author declares no conflict of interest. The funder had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


## *Article* **Economic Indicators in Water and Wastewater Sector Contributing to a Circular Economy (CE)**

**Marzena Smol 1,\* and Renata Koneczna 1,2**


**Abstract:** Protection and sustainable management of water was indicated as one of the strategic tasks in the process of transformation towards a circular economy (CE) in the European Union (EU), therefore, the water and wastewater sector plays an important role in this process. At the same time, the European Commission (EC) strongly underlined the importance of the possibility to assess the transformation process toward the CE, and developed a set of CE indicators that are available on the Eurostat website. However, these indicators have limited ability to assess the transformation progress in the water and wastewater sector. This paper presents a set of indicators for assessing the economic progress of transformation towards the CE in this sector. The proposed economic CE indicators were grouped into the following actions of the CE model in the water and wastewater sector: reduction, reclamation (removal), reuse, recycling, recovery and landfilling. The selection of specific indicators was based on a systematic review of the literature presenting economic indicators developed by international organisations and researchers (covering different thematic areas, scopes and potential applications). The selected economic CE indicators were assigned to three groups of the cash flow: income (revenues, expenses), costs, and investment financing. The proposed CE indicators can be used by water supply and sewage companies (i.e., supplying water to the public and wastewater treatment plants, and companies that use water in their production processes) to assess the level of the transformation toward the CE at a microeconomic level. An important aspect of future application and usage of the proposed set of CE economic indicators is the collection and processing of data needed for their reporting. The proposed set of CE economic indicators refers to information that are reported by the companies to prove its revenues, costs and investment outlays, and are collected by companies anyway. The proposed set of economic CE indicators is flexible, allowing the adaptation of indicators and areas of interest to maintain effectiveness throughout the transition period from linear to the CE model.

**Keywords:** circular economy (CE); monitoring; indicators; economic indicators; water; wastewater

#### **1. Introduction**

A circular economy (CE) is defined as a regenerative system [1] where the value of materials, products and resources is maintained as long as possible in the economy and the production of waste is minimized [2]. The CE enables more efficient use of available resources, but also promotes a more sustainable management of waste. The integrated initiatives along the entire life cycle of raw materials [3], from extraction to the circular final processing are more and more often and successfully introduced in various industries [4]. It has to be underlined that the CE refers not only to raw materials ( such as animal, vegetable or mineral), but also to water [5], which is an irreplaceable resource with life-giving property for nature, people and the economy [6].

Water resources are currently under unprecedented pressure in most countries [7]. The problem of water stress (reaching the level above 70%) occurs mainly in the regions of

**Citation:** Smol, M.; Koneczna, R. Economic Indicators in Water and Wastewater Sector Contributing to a Circular Economy (CE). *Resources* **2021**, *10*, 129. https://doi.org/ 10.3390/resources10120129

Academic Editor: Diego Copetti

Received: 2 November 2021 Accepted: 13 December 2021 Published: 16 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the world such as Northern Africa, Middle East, Western, Central and Southern Asia [8]. However, this problem also affects the European Union (EU), as the water scarcity was estimated to have affected at least 17% of the EU territory and at least 11% of the European population [9]. Moreover, next to water scarcity, an important issue in water management in Europe is water pollution, from industry and agriculture. To minimize the effects of anthropogenic use of water, both in agriculture and in industry, the European Commission (EC) announced further initiatives focused on water management in the second CE Action Plan [10]. This new Action Plan is one of the main blocks of the European Green Deal (EGD)-new agenda for sustainable growth of the EU. The main goal of the EGD is to achieve climate neutrality in Europe by 2050 [11] by turning climate and environmental challenges into new opportunities across all policy areas, and ensuring that the transition is fair and inclusive. In the coming years, the EC plans to facilitate water reuse and efficiency in both industrial processes and agriculture [10].

The protection and sustainable management of water and water-based waste (as wastewater, sewage sludge or sewage sludge ash) are indicated as one of the strategic tasks in the transformation process towards CE [12]. In the White Paper "Water and the Circular Economy" [13], common characteristics, ideas and approaches between the CE initiatives being implemented by organizations and Water System Management were identified. The three key dimensions of water use were grouped into the three themes of: (i) water as service (consumptive use, production use, process use), (ii) water as source of energy (kinetic, thermal, bio-thermal), and (iii) water as carrier (nutrients, chemicals, minerals). The main areas of the transformation of water and wastewater sector to the CE model have been also indicated by the International Water Association (IWA) [14]. The IWA also proposed three pathways to support the water utility leaders in boosting their progress towards CE: (i) water pathways (upstream investments, rainwater harvesting, water recycling for non-potable use, water reuse for agriculture/aquaculture/industry), (ii) materials pathways (resource efficiency, used water sludge and products for agriculture, bioplastics, fertilizer and other materials), and (iii) energy pathways (energy saving, energy reduction and recovery, biosolids to energy production, renewable energy). In turn, the EC assumed that in applying the CE main principles—reduce, reuse and recycle—in the water and wastewater sector will accelerate the process of transition to the CE model in the EU. However, at the moment it is not possible to determine whether the subsequent actions are bringing the intended effects, whether environmental, social or economic. This is due to the lack of a dedicated CE monitoring framework for the water and wastewater sector, which would take into account indicators and measures allowing the assessment of the level of transformation towards CE in this sector. In 2018, the EC proposed the CE monitoring framework with ten CE indicators grouped in four thematic areas: production and consumption, waste management, secondary raw materials, competitiveness and innovation [15]. However, the potential application of these indicators in water and wastewater sector is limited, and does not evaluate all sector elements. Despite the fact that the EC underlined that the monitoring such important areas as production and consumption is essential for understanding progress towards the CE, the presented data does not take into account the water usage. Moreover, apart from the environmental benefits (resulting from the protection of water resources) and social benefits (securing drinking water supplies), the economic benefits of taking measures to implement CE in the water and wastewater sector should be also demonstrated [16]. In this area, the EC indicated that water savings in all sectors in the EU could lead up to 5% of reduced total primary energy consumption, which bring economic benefits for individual players [4]. To encourage companies to implement CE measures, which could generate greater value and commercial opportunity [13], economic indicators that allow for the assessment of the level of transformation towards the CE in different sectors should be identified. Therefore, the objective of the current paper was to present an inventory of economic indicators that can be used for the evaluation of the progress toward CE in the water and wastewater sector. In the previous years, several indicators have been proposed to access water-related aspects in the economy; however, their goals, scope and potential application for the assessment of level of transformation toward CE in water and wastewater sector must be analyzed and evaluated from the point of view of the possibility of their monitoring at the microeconomic and macroeconomic levels.

#### **2. Materials and Methods**

The current research includes three steps of collecting and processing data. The research framework is presented on Figure 1.

**Figure 1.** The research framework.

In the first step, a detailed analysis of the published research was conducted with the use of the desk research method. This state of the art analysis was based on the review approaches used in [17,18] to conduct searching and eligibility screening of available literature while retaining the procedural scope of analysis, and ensuring that that the review process is objective, repeatable and. The objective of this step of research was to review the indicators (economic, social, technological and environmental) from national and international organizations. The analyzed indicators regarded different aspects related to the CE and sustainability. Moreover, from the list of identified CE-related indicators, the specific indicators that can be used in the water and wastewater sector were also analyzed. The following data sources were analyzed: international and European official documents related to the water and wastewater sector and circular economy, published in EUR-Lex (eur-lex.europa.eu) and the official webpage of the EC (ec.europa.eu). The analysis also included the review of the statistical documents at the international and European levels (Organization for Economic Cooperation and Development—OECD, United Nations—UN, World Bank, European Investment Bank, Eurostat, European Environment Agency—EEA) and selected reviewed articles available in the scientific databases Elsevier Scopus, Elsevier Science Direct, Google Scholar and the Multidisciplinary Digital Publishing Institute (MDPI) database [19,20]. The selection of the articles was conducted based on the identified keywords "circular economy", "CE", "economic" "indicator", "index", "measurement", "assessment", "water", and "wastewater". The results of this step of research are presented in Section 3.1.

In the second step of the research, identified indicators, measures and indices have been analyzed and grouped according to the CE model for the water and wastewater sector, developed under the MonGOS project [14]. The objective of this step of research was to propose a set of economic indicators that could be used in water and wastewater

management. At the beginning, the economic indicators were selected from the list of indicators analyzed in the first step of the research. Social and environmental indicators have been rejected as they are not the subject of the current research. Then, these indicators that are directly or indirectly related to the measurement of economic efficiency were analyzed and grouped in the following actions of the CE model in the water and wastewater sector: *reduction*, *reclamation* (*removal*), *reuse*, *recycling*, *recovery* and *rethink* [14]. The concept of the CE model framework in the water and wastewater sector is presented in Table 1. Finally, new economic indicators have been proposed for each action of this model. The round table discussion which included the consultation in the group of internal and external experts was used for this purpose. The research group consisted of 6 expertsthree representatives of enterprises operating in the water and wastewater sector, and three specialists (scientists) dealing with economics and environmental technologies. The criterion for the selection of experts from enterprises was to have at least a master's degree and a minimum of 10 years of experience in a managerial position in a company using water or/and dealing with wastewater treatment. In turn, the criterion for selecting scientific experts was to have a doctoral degree in the economic and environmental sciences, taking into account experience in the water and wastewater sector. The indicators were analyzed and discussed during three meetings with these experts: (1) consultation online with industry experts, (2) consultation online with scientific experts; (3) consultation online during joint meeting of the MonGOS project. The final results of this step of the research was the list of economic indicators that can be used for the evaluation of the level of the transformation toward the CE in water and wastewater sector. The results of this step of research are presented in Section 3.2.

**Table 1.** CE model in the water and wastewater sector (own based on [14]).


The third step of the research includes the discussion of the possibility of the usage of the identified economic indicators at the microeconomic levels. The synthesis method, that is, formulating generalizations based on recognized partial theorems, was used to interpret the obtained results. Moreover, the desk research was used to compare obtained results from the perspective of previous studies and other authors. The findings and their implications, as well as future research directions, are discussed in Section 4.

#### **3. Results**

#### *3.1. Inventory and Classification of CE Indicators*

In recent years, several CE indicators have been developed by various organizations. The inventory of the groups of circularity indicators in the documents of international organizations is presented in Table 2. As a result, 10 documents proposing or discussing CE indicators at the macro (European and international monitoring framework) and micro (products, services and companies) levels were presented. There are also CE indicators proposed by individual authors. Therefore, to the listed documents, selected scientific papers presenting specific CE indicators were also analyzed (Table 3). In total, 742 indicators (provided by international organizations) and summary indicators and indexes (provided by other authors) were analyzed.


**2.** Groups of circularity indicators in the documents of international organizations


**Table 2.** *Cont*.


**Table 3.** Circularity indicators in scientific papers.

There are various methods of grouping the indicators into specific thematic groups, for example, according to the perspective of sustainability—economic, social and environmental indicators proposed by the Global Reporting Initiative (GRI) [21] or due to individual 17 Sustainable Development Goals (SDGs) proposed by United Nations (UN) [22]. A holistic picture of the level of transformation towards a CE in European countries is indicated by the CE monitoring framework, developed by the EC in 2018 [15]. It proposes four groups of indicators, divided according to the key areas of CE implementation in the EU, such as production and consumption, waste management, secondary raw materials and competitiveness and innovation. Other organizations group indicators according to specific industries, environmental problems, individual elements of the life cycle [23], or components of the environment [24]. Depending on the potential application, there are different goals with different scope with regard to the proposed specific indicators.

As can be seen from Table 2, there are various classifications of indicators related to the features of a circular economy, focusing on the assessment, improvement, monitoring and communication on the results of the CE [25]. However, there are no official or recognized indicators, methods or tools for measuring a company's performance in the transition from a linear economic model to a more sustainable one, and there are no tools to support and track this transition [26]. Indeed, most of the CE indicators are in their early stages of development [27], and they cannot capture the overall performance of circular products and services [25]. However, many existing indicators can help measure performance in several areas (micro, meso, macro) that contribute directly or indirectly to the circular economy [28].

There are also interesting indicators developed by the individual authors. The summary of these indicators is presented in Table 3, taking into account the level of their measurement (micro, meso, macro). Most of these indicators can be grouped to the perspective of sustainability-economic, social and environmental. Most of the presented indicators focus on individual stages or all stages of life cycle, such as the Circularity Measurement Toolkit (CMT). They mainly measure the efficiency of materials use.

From the list of analyzed indicators, the CE-related indicators that can be used directly or indirectly for the circularity analysis in the water and wastewater management sector were selected. The conducted inventory shows that direct CE indicatos published by international organizations mainly focus on single indices for the water and wastewater sector not rotating the entire CE model. For example, the UN presented 231 indicators, including 11 water and 75 economic indicators, while the Organisation for Economic Cooperation and Development (OECD) published 153 indicators, including 17 water and 56 economic indicators.

CE-related indicators for the water and wastewater sector were published by the EC—water exploitation index, water productivity, price of water scarcit and water use (calculated as water abstraction minus returned water), the European Environment Agency (EEA)—water exploitation index, water productivity. Additionally, the World Business Council for Sustainable Development (WBCSD) presented the water circularity index (%), and the OECD presented a water productivity index. The European Investment Bank

(EIB) proposed circular value recovery models: reuse/recycling of wastewater. However, despite the proposed CE indicators, there are no statistical data that would allow them to be calculated and reported. The proposed indicators refer to both the micro and macro levels of water management. However, none of the organizations proposed a set of indicators that could assess the level of the CE transformation in the water and wastewater sector.

On the other hand, the indicators published in scientific papers only take into account individual aspects of the resources (incl. water) management, as cost analysis without analysis of possible revenue for CE. The economy aspect is combined with others, e.g., social or environmental. Only the Eco-costs Value Ratio (EVR) is an indicator of the economic dimension. The authors propose to calculate all environmental effects in monetary terms based on the costs that should be incurred to reduce environmental pollution and materials depletion to "no effect level" [39]. However, the EVR does not take into account the income that a company or household may gain as a result of involvement in CE.

The revision of the available CE-related indicators shows that there is a lack of a set of economic indicators that could be used for measurement of the level of transformation towards CE in the water and wastewater sector. Therefore, the next sections provide a proposal for a set of CE -related economic indicators in mentioned sector of the economy.

#### *3.2. CE Economic Indicators in the Water and Wastewater Sector*

As part of the MonGOS project, the assumptions for the CE model in water and wastewater management was proposed and published in [14]. The assumptions for this model were developed on the basis of the "xR" models in waste management, as well as the EU waste hierarchy [42]. The CE model in water and wastewater management has been classified into groups of activities that fit into the assumptions of the CE, i.e., *reduction*, *reclamation* (*removal*), *reuse*, *recycling*, *recovery* and *rethink* [14]. In the current section, the economic indicators have been proposed for each action of this model. The selected economic indicators were assigned to three groups of the cash flow (income revenues, expenses – costs, and investment financing) and to the specific actions of the of the CE model in water and wastewater sector (*reclamation*/*removal*), *reuse*, *recycling*, *recovery)* and landfilling. The proposed indicators are dedicated to water supply and sewage companies, i.e., supplying water to the public and wastewater treatment plants (WWTPs) and companies that use water in their production processes. The proposed CE economic indicators in the water and wastewater sector are presented in Table 4.

The first element of the model (*reduce*) includes revenues from less water consumption. It creates two levels of value added-lower charges for water abstraction from the water supply, but also less pressure on the environment (lower water consumption and then less amount of generated wastewater that need to be treated). Investments are related to the cost of equipment for the optimization of water usage.

In the second element of the CE model (*reclamation*/*removal)*, the revenues come from the possible sale of water, the provision of wastewater collection and treatment services, to ensure the continuity of collective water supply of adequate quality and quantity, and collective wastewater disposal. The expenses are related to costs of water production and wastewater treatment services. The investments are related to equipment for water purification and wastewater treatment. In addition, there are also costs of water intake and abstraction, operation, maintenance and expansion of the water supply and sewage.

In the next level of the CE model (*reuse*), the revenues come from sales of non-drinking water, lower wastewater treatment services fees and reduced water abstraction from the waterworks. The expenses include the costs of non-drinking water production (costs of non-consumer water recovery). There are also costs of electricity production, water consumption, external services and employee remuneration, as well as investments in equipment for the wastewater treatment and water recovery. The economic added value of implementation of water reuse activities are lower annual water bills, an increase in additional business entities benefiting from improved wastewater treatment, and less waste landfilled. Analogous economic indicators are proposed for the next element of the CE

model (*recycling*). However, in this case, due to the need for additional water purification (for human consumption), more energy and material costs will be incurred. On the other hand, higher revenues from sales of drinking water are expected.


**Table 4.** Proposed CE economic indicators in the water and wastewater sector.

In the next element of the CE model (*recovery*), the revenues are related not only to lower wastewater treatment services fees, but also the ability to sell electricity, fertilizers and other materials. Here, the expenses include costs of electricity, fertilizers and materials production, while investment are related to equipment needed for recovery of those resources.

It should be noted that the implementation of technological and organizational solutions in each of element of the presented CE levels can reduce the amount of landfilled waste and thus bring an economic benefit. In addition, emissions and environmental charges for landfilling are reduced. The proposed economic CE indicators selected can allow companies to identify both positive and negative trends in their activities. The measurement of these indicators can provide information about which actions are effective or created negative effects or regression at the CE economic level. The implementation of CE-related solutions is expected to increase in coming years, due to the tightening legal regulations on water and water-based waste, and EC recommendations for the transformation toward the CE model.

An important aspect of the future application of proposed indicators is the collection and processing of data needed for the reporting of these indicators. The proposed CE economic indicators are related to the information that is reported by the companies to demonstrate their revenues, costs and investment outlays. Therefore, the collecting of this economic data should not face significant barriers in the individual units.

#### **4. Discussion**

The transformation process towards the CE requires more rational use of resources and waste management practices in all sectors of the economy [43]. This also applies to the water and wastewater sector and its key elements, i.e., water, sewage, sewage sludge, other waste, and by-products arising from water purification and wastewater treatment [14]. In practice, the implementation of the CE assumptions in various sectors of the economy is often supported by the use of various methods of rational management of raw materials, products and resources, as well as sustainable waste management [44].

The process of transformation towards the circular economy in the water and wastewater sector [45] requires the involvement of all stakeholder groups, both experts working for innovative and pro-ecological solutions, and the society, which should reduce water waste in households. In addition, the implementation of circular economy principles in water and wastewater management is important for enterprises dependent on water (e.g., the cosmetics industry) and wastewater treatment plants, because their environmentally conscious decisions regarding the implementation of sustainable and circular solutions in the management of water and wastewater may accelerate the transformation process towards the CE. in the given country.

In the recent years, some significant progress has been made in the area of the assessment of circularity of products [46–48], companies [34,49,50], and regions [51,52]. The CE indicators are created to assess the progress of transformation towards CE at the micro, meso and macro levels. They are an important element of new business models for CE, which are systematically improved and introduced into the activities of enterprises, including those operating in the water and wastewater sector. The EC clearly indicates that the transition to the CE model brings economic benefits for those involved in the transformation process [10]. In order to assess the economic benefits of implementing CE solutions, specific financial data on the operation of the company and introduced changes must be reported. Therefore, the indicators developed in this research can be used by the enterprises to assess the level of transformation towards CE in the water and wastewater sector. The proposed indicators refer to the CE model for the sector developed in the MonGOS project [14], thus providing a broader perspective of the sector. In addition, they require reporting of information that is collected by the enterprise anyway, therefore its collection and processing should not pose a significant challenge to individual entities.

The added value of the presented economic indicators is the possibility of their application in various enterprises operating in the wastewater sector, i.e., supplying water to the public, WWTPs, and other companies (public and private) that use water in their production processes and create an innovations for the water and wastewater sector (e.g., Schwander [45], Veolia [53] or Outotec [54]). The proposed economic indicators could measure the CE-related activities in water and wastewater management, as minimization of water consumption, water and wastewater treatment, water reuse (for non-consumption purposes), water recycling (for consumption purposes) and the recovery of water, energy and raw materials produced in the water and wastewater treatment processes. The implementation of those activities may bring significant environmental benefits, resulting from the reduction of water consumption and the reduction of the impact of wastewater discharge on the quality of the aquatic environment. In the following years, further technological progress and new investments should be expected to reduce the consumption of water, raw materials and energy, in line with the CE model. Therefore, the proposed economic indicators can be widely applied in the sector and can complement new business models for CE.

Moreover, the usage of the proposed set of economic indicators at multiple levels would facilitate policy development, measuring economic performance, sector benchmarking, and improving business investment decisions. Such a framework should provide meaningful answers to decision-makers questions covering all relevant dimensions of the CE transition: resource consumption and material flows, economic parameters, financial flows and policy effectiveness. The presented set of economic indicators is flexible, allowing the adaptation of indicators and areas of interest to maintain effectiveness throughout the transition period. Further work on testing the developed indicators in the individual companies is undertaken as part of the MonGOS project.

#### **5. Conclusions**

Water and wastewater management is an important part of the CE model. The circular management of resources (as water, nutrients, energy) could generate financial benefits for companies operating in the water and wastewater sector. The financial benefits could be related to lower water consumption, the sale of drinking water, electricity, fertilizers, and lower fees for wastewater treatment services and taxes for waste landfilling. However, there are also unavoidable costs of the CE-related activities as investment outlays in new infrastructure. In the long run, there should be a reduction in average total costs that could bring further economic benefits.

The implementation of the CE solutions requires estimating the costs of such activities, both in terms of possible expenses and revenues. Therefore, the proposed set of economic CE indicators can be used by the companies not only to assess the transformation level but also to plan the future CE-related investments.

**Author Contributions:** Conceptualization, M.S. and R.K.; methodology, M.S.; investigation, M.S. and R.K.; resources, M.S. and R.K.; data curation, M.S. and R.K.; writing—original draft preparation, M.S. and R.K.; writing—review and editing, M.S. and R.K.; visualization, M.S. and R.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Polish National Agency for Academic Exchange (NAWA) as the part of the project "Monitoring of water and sewage management in the context of the implementation of the circular economy assumptions" (MonGOS, mon-gos.eu).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

