Modelling of climate change impact on flow conditions in the lowland anastomosing river

Site description

The Upper Narew is the sub-catchment of the largest Polish river basin—the Vistula (Fig. 1). It drains an area of 4,231 km² (of which 27% belong to Belarus) and is characterized by a flat relief with an average elevation of 152 m a.s.l. The prevailing type of soils in the catchment are sands and loamy sands, whereas heavy, impervious soils are rare. The land cover in the Upper Narew catchment is predominantly forested (43.6%) and agricultural (41%), whereas wetlands and grasslands occupy 16% of the area. The catchment is located in a temperate climatic zone with moderately warm summers (mean temperature in July 18 °C), cool winters (mean temperature in January −2 °C), and an annual average precipitation total of ca. 600 mm.

The Upper River Narew stretch, which was investigated in this study (Fig. 1C), is characterized by low-gradient (0.0002 m/m) and anastomosing pattern. It is situated in the Narew National Park (NNP) for which Natura 2000 sites (under both Bird and Habitat directives) were established. Within the NNP (6,810 ha), the river is characterized by a complex network of small interconnected channels (Gradziński et al., 2003). Vegetation cover within the NNP is dominated by early growth sedge and reed communities. In-stream vegetation densely overgrows channels during the summer season, but also survives the dormancy season to some extent, influencing the hydraulics of channels for the whole year. The flow regime (based on discharge data series from the Suraż gauging station spanning from 1951 to 2018) is characteristic for lowland snow/rainfall-fed rivers with clearly dominating spring flooding, sourced by snowmelt occurring most commonly in April. Summer floods caused by heavy rain events are significantly lover. The average yearly discharge for the analyzed period is 15.4 m³/s, whereas extreme events (i.e., minimum and maximum flows) equal 1.5 m³/s and 250 m³/s, respectively. The number of days with flooding is diverse (from 0 to 119; 37 on average) and depends majorly on the snow cover thickness in winter season.

This particular stretch has been selected for the analysis, as according to Marcinkowski, Grabowski & Okruszko (2017), the channel extinction process was recognized to be in progress there for the last few years. It is worth noting that in current conditions during the growing season on average, nearly 98% of the flow is concentrated in the main channel (right-hand side reach on Fig. 1C) and only 2% is distributed to left-hand side anabranch. Such extremely low flows and unevenly distributed discharge create favorable conditions for vegetation encroachment and channel sedimentation, which in turn lead to the gradual closure of the side channels, as has been observed in the field for the last few years. Therefore, it has to be highlighted that if current hydrological conditions, especially low flows during the vegetation period, remain stable in the coming decades, the anastomosing system will most likely shift into a single-channel planform. The extinction process is in progress now and it might be assumed that it will be continued unless the low flows increase.

Integrated modeling framework

The IMF applied in this study consists of a catchment scale SWAT model and a river reach scale HEC-RAS model. The Hydrological SWAT model driven by climate change projections was used for daily discharge calculations for the future time horizons: 2021–2050 (near future—NF), and 2071–2100 (far future—FF). For a baseline scenario, which determines a reference point for the projected climate changes, historical (1971–2000) bias-corrected climate model data was used. In the study, flow condition changes were assessed at two different spatial scales: (1) at catchment scale mean monthly flows were calculated (for all SWAT sub-basins) showing the general impact of climate change on water resources in the area and (2) at reach scale daily hydrographs (implemented in 1D HEC-RAS hydraulic model as an upstream boundary condition and groundwater inflow from subcatchment) were calculated from the SWAT project. This way SWAT-based simulations of daily flow could be directly transferred and routed in the HEC-RAS stream network (Fig. 1C) as the flow hydrograph input, which allowed assessing the impact of climate change on hydraulic conditions in the anastomosing section of the River Narew. The projected flow hydrographs of models’ ensemble as a SWAT model output can be found in the raw Supplemental Material (flow simulation files). As low flows prevail in summer and autumn months, triggering anabranch loss (Marcinkowski, Grabowski & Okruszko, 2017), HEC-RAS is only applied during these periods. The resulting daily discharge was next reach-averaged and analyzed for flow velocity. Additionally, to assess not only the magnitude but also the duration of low flow conditions, hydrographs derived from SWAT were statistically examined. Given that the channel loss process is triggered by low flow conditions, a threshold value has been defined based on the daily discharge observations from Suraż gauging stations (1952–2017), expressed as Q₈₀ (with 80% exceedance probability). This particular exceedance probability has been selected for two reasons: (1) according to Smakhtin (2001), flows within the range of 70–99% time exceedance are usually most widely used as design low flows, (2) collected discharge measurements used for HEC-RAS model calibration correspond to Q₈₀. Taking into account the Q₈₀ flow (7 m³/s), for each climate model in the baseline, NF and FF horizons, the number of days lower than Q₈₀ were accumulated. Additionally, the mean value of flow of consecutive days below threshold was calculated. The analysis scheme is presented in Fig. 2.

SWAT

SWAT is a process-based, continuous-time model which simulates the movement of water, sediment and nutrients on a catchment scale (Arnold et al., 1998). It is a comprehensive tool suitable for investigating the interaction between climate, land use and water quantity. It enables simulation of long-term impacts of climate changes on water, sediment and nutrient loads in catchments with varied topography, land use, soils and management conditions. Spatially, the catchment area is divided into sub-basins and further into Hydrologic Response Units (HRU), which are the smallest spatially explicit units. Water balance as well as sediment and nutrient loads for each HRU are calculated at the land phase. Runoff is further aggregated to the sub-basin level and routed through the stream network to the main outlet in order to obtain the total runoff for the river basin. SWAT has also a simplified groundwater flow component which partitions groundwater into two aquifer systems: a shallow, unconfined aquifer which contributes return flow to streams and a deep, confined aquifer which contributes return flow to streams outside the watershed. In this study, we build upon the existing SWAT model of the Upper Narew catchment (Marcinkowski et al., 2016).

Model setup, calibration and validation

Delineation of the catchment based on the 10-m resolution DEM resulted in division of the catchment into 243 sub-basins. The intersection of land cover map, soil map and slope classes resulted in creation of 4,509 HRUs. Daily precipitation and air temperature (minimum and maximum) data (1951–2013) were acquired from 5 km resolution gridded, interpolated using kriging techniques, dataset (CPLFD-GDPT5) based on meteorological observations coming from the Institute of Meteorology and Water Management (IMGW-PIB; Polish stations) (Berezowski et al., 2016). The use of interpolated climate data in the SWAT model, as reported by Szcześniak & Piniewski (2015) increased the model performance for a case study in Poland.

The calibration phase was conducted in SWAT-CUP using the SUFI-2 algorithm (Sequential Uncertainty Fitting Procedure Version 2) where the Kling–Gupta efficiency (KGE) was used as an objective function (Gupta et al., 2009). In the calibration and validation of the daily discharge, 10 flow gauges (data acquired from IMGW-PIB) were used. The calibration period was set to 1976–1985, and the validation period was 1986–1991, which was dictated by the highest number of observation records available in the catchment (Marcinkowski et al., 2016). Marcinkowski et al. (2016) reported satisfactory values of goodness-of-fit measures across for all gauging stations (median KGE above 0.7). Figure 3 depicts the observed vs simulated daily discharge in the Suraż gauging station (cf. Fig. 1), which is at immediate upstream to the HEC-RAS model created for the anastomosing section of the Narew River.

HEC-RAS

HEC-RAS is a well-known hydrodynamic model for one and two-dimensional hydraulic calculations and is often used in calculating both steady and unsteady, gradually varied flow (Dysarz et al., 2019; Horritt & Bates, 2002). In addition, the HEC-RAS system solves movable boundary sediment transport computations and water quality analysis. In 1D hydrodynamic component model solves the full 1D St. Venant equations for unsteady open channel flow with the finite difference method, using a four point implicit box scheme (Brunner, 2016). Manning’s roughness coefficients are used to describe the flow resistance and were calibrated for the analyzed river branches. HEC-RAS performs hydraulic calculations for a full network of natural and constructed channels. Depending on the type of reach junction (flow split or flow rejoin) in a river network, one equation of continuity of flow or equation continuity of stage is used as an interior boundary equation.

The catchment scale hydrological response to climate change

The projected mean annual temperature in the Upper Narew catchment is expected to increase by approximately 1.2 °C (ensemble range 0.8–1.6 °C) in NF and 2.1 °C (ensemble range 1.5–2.6 °C) in FF following the RCP4.5 according to the ensemble mean. Comparing the monthly variation, the highest change is expected to occur in November–March (2.4 °C, ensemble range 2.1–2.6 °C) and lowest in April–October (1.7 °C, ensemble range 1.5–1.9 °C) (Marcinkowski et al., 2017).

The annual total precipitation is projected to increase by 5.6% in the NF and by 9.5% in the FF, on the annual basis. The seasonal patterns show a relatively high increase in winter and spring, and a mixed response (both increase and decrease cases) in the summer and autumn. In the FF, the spring precipitation increase is distinctly higher than in other seasons, exceeding 20% (Marcinkowski et al., 2017) (Fig. 4).

The SWAT model simulations indicate that the median of projected changes in average water discharge, on an annual basis and catchment scale, is expected to increase by 11% and 25% in the NF and FF respectively. Monthly patterns showed the most pronounced increase occurring in December–March, with the ensemble median reaching 17% and 35% for NF and FF respectively. In the April–November period, a significantly lower increase in discharge is observed, reaching 7% and 23% for the NF and FF, respectively (Fig. 5). Additionally, more extreme and outlying values (both low and high) are noted for the future time horizons, compared to historical period. High flows might be attributed to the occurrence of extremely heavy rain events in the projected rainfall patterns, and low flows to an increased air temperature triggering the increase of evapotranspiration.

HEC-RAS model performance

A simulated split flow at the channel junction was compared with observations in the NR and NL reaches. In the calibration process for a discharge of 7.4 m³/s, the error of the flow distribution is 0.7% and 11% for section NR and NL respectively. The average error of the water table in the river reaches varies from 0.03 m to 0.05 m. The mean error of simulated water velocity is 6%, 7% and 14% for the ND, NR and NL segments respectively (Table 1). The model accuracy discussed above was obtained for the following set of Manning’s coefficients: n_NU = 0.028; n_NR = 0.028; n_NL = 0.1; n_ND = 0.035 m^−1/3/s. These values are in the range of the roughness coefficients used in other models for the analyzed section of the Narew River (Marcinkowski, Kiczko & Okruszko, 2018). Verification of the developed model was based on measurements of the water level and water velocity performed at a discharge of 10.5 m³/s. In this case, the errors for the flow distribution are 0.9% and 9% for the sections NR and NL respectively. The average error of the water table varies from 0.04 m to 0.07 m. The mean error of simulated water velocity is 7%, 3% and 9% for the ND, NR and NL river reaches respectively (Table 2). The obtained errors are considered acceptable to allow the use of the model to simulate the unsteady flow during the vegetation period. Then, the developed hydrodynamic model was used to simulate the flow for the long historical period (1973–2000) and climate change scenarios for the years 2070–2100. In these simulations, the upper boundary condition in the form of a flow hydrograph was calculated in the SWAT hydrological model.

Reach scale changes in flow condition

At the reach scale the flow condition changes were calculated based on the SWAT model simulations for the particular sub-basin, which spatially overlaps with the HEC-RAS hydraulic model. As indicated in Fig. 6, the duration of low flows changes significantly when comparing the historical period and future horizons. Considering the annual values (Fig. 6B), a gradual decrease from 151 (baseline) to 135 (NF) and 123 days (FF) is observed. However, the change is not uniform across the year and varies notably between each month (Fig. 6A). In general, for November–July, a clear drop in the duration of low flows in observed reaching 30% and 50% (ensemble median) for NF and FF respectively. However, for August–October, a slight increase in the duration of low flows is noted (8%) for the NF. Such a trend subsides in FF where again the decrease is observed, but significantly lower than for winter, reaching only 20%. The analysis of mean flow values of consecutive days below threshold indicate a gradual increase at the annual level from the historical period to future horizons (Fig. 6D). However, the change is not uniform across the year. In general, a higher increase is observed for FF compared to NF, and the highest increase occurs for February–May in both time horizons. What is noteworthy, the mean value of flow for days below Q₈₀ seems to stay at almost the same level in NF for August–October (Fig. 6C).

Further analysis included assessment of low flow magnitude changes, which for the summer and autumn seasons was modeled in HEC-RAS using daily flow hydrographs derived from SWAT. Reach-averaged discharge and flow velocity is presented separately for the main channel and the anabranch. Simulations derived from SWAT indicate that both the discharge rate and mean flow velocity slightly increase in NF (3% for ensemble median) and FF (5% for ensemble median) (Fig. 7). Based on the T-student test, NF changes in discharge and flow velocity were not statistically significant (p value of 0.16 and 0.15 respectively) and the FF changes were statistically significant (p value of 0.0001 and 0.001 respectively).

Anastomosing river perspectives under climate change

The conservation of anastomosing rivers worldwide, although extremely important nowadays, has not been addressed in many studies. The success of their conservation depends on finding and steering the appropriate key factors, which differ from stream to stream and site to site (Verdonschot & Nijboer, 2002). Marcinkowski, Grabowski & Okruszko (2017) pointed out the low flow to be the main stressor triggering the loss of anastomosing channels in the NNP. In this study, flow conditions were analyzed in terms of duration and magnitude by using a cascade of modeling driven by climate change projections. It was due to the fact that Regional Climate Model (RCM) data operates at coarse scale, whereas in river conservation programs, particular river reaches of fine scale are targeted. The results were presented at catchment (using SWAT) and reach (using HEC-RAS) scale on an annual and monthly basis. Catchment-averaged mean annual duration of low flows (Q below Q₈₀) was projected to decrease and mean flow to increase in both future time horizons. Conclusions drawn on an annual basis results suggest that hydrological conditions might improve in future. However, more detailed monthly analysis indicated that for August–October, when according to Marcinkowski, Grabowski & Okruszko (2017) the loss of anastomoses is the most intense, duration of flows below Q₈₀ increases and the mean flow in NF is projected to remain at the same level as it is currently. The projected increase in precipitation and flow is distributed mainly to the winter and spring seasons, when vegetation is scarce and the potential of overgrowing the anastomoses very low. In consequence, near future hydrological conditions are not projected to improve in the most sensitive (in terms of anastomoses loss) period of the year.

The reach-scale HEC-RAS simulations of water flow conditions in the anastomosing sections of the river Narew indicate a very low flow velocity (0.09 m/s–ensemble median) in the anabranch, which remains the same in NF and only slightly increases in FF (0.1 m/s—ensemble median). Such hydraulic conditions lead to in-stream vegetation development. Typically the abundance of macrophytes are stimulated at low velocities and their growth is restricted at higher velocities (Marcinkowski, Grabowski & Okruszko, 2017; Marcinkowski, Kiczko & Okruszko, 2018). This is due to the fact that water velocity affects photosynthesis through the regulation of the availability of dissolved substances (Madsen et al., 2001). It has been reported that the photosynthetic and nutrient uptake rates of in-stream macrophytes are positively correlated with low water velocities (0–0.1 m/s) (Marshall & Westlake, 1990; Madsen et al., 2001). The prospect of a further reduction of river flow velocity in NF horizon suggest further degradation on the anastomosing system in coming decades.

One of the major issues in freshwater ecosystems conservation is that adaptive management must be applied on time to be useful. Therefore, the delayed detectability of a climate change signal might cause the implemented practices to be ineffective or bring economic and ecological losses. It is well known that the climate models themselves, although continuously improved, are to some extent uncertain (Ito et al., 2020). However, reasonable conservation management may consider the use of model projections to guide future interventions. In the anastomosing Narew case, the simulations at reach scale, which were run taking climate projections and process-based catchment modeling into consideration, indicate that most possibly by the end of 2050 (NF) hydrological conditions and consequently the sediment transport potential, will not improve. This means that by that time most anastomosing channels in the NNP will be aggraded and fully colonized by vegetation. Marcinkowski, Kiczko & Okruszko (2019) conducted a model-based assessment of the hydraulic efficiency of channel vegetation removal, dredging and water damming, in increasing stream capacity and sediment transport potential in the anastomosing section of the NNP. The results of this paper indicated that summer flow conditions in NF will be comparable to the current state, which means that the flow boundary conditions used in the model of Marcinkowski, Kiczko & Okruszko (2019) are applicable to NF flow conditions. Consequently, their simulations of stream capacity and sediment transport potential might be used in this study to accurately reflect future hydraulic conditions and potential efficiency of conservation measures. The projected hydraulic conditions suggest that the least invasive in-stream vegetation removal in the future will not increase the sediment transport potential preventing river bed aggradation and consequently, stream colonization by plants. A sufficient increase in the stream flow capacity might be achieved only by implementing high invasive channel dredging and water damming, provided that they are implemented relatively soon. Such alterations of hydraulic properties in anastomoses might result in a significant increase in sediment transport potential, which ensures the stability of restoration eliminating the main cause of channels extinction that is, excessive sediment deposition and plants colonization.

Study limitations

In this study, we used coupled hydrological and hydraulic models to test the impact of climate change on flow conditions in the anastomosing section of the NNP. The designed methodological approach as well as tools and the datasets used are burdened with some limitations. First of all, it has to be mentioned that low flows, which are of special concern in this study, are mainly triggered by groundwater, whereas both the SWAT and HEC-RAS models do not address this flow component in sufficient detail. To reduce the errors while transferring the projected discharge values from SWAT to HEC-RAS beside the daily stream hydrograph, a groundwater flow (generated by groundwater recharge) from the direct sub-catchment area of the anastomosing section, was added to the HEC-RAS domain. According to Bates et al. (2008), hydrological changes to European rivers caused by climate change are strongly related to the stream flow type. Following their findings for catchments dominated by groundwater flow (Narew case), a decline in groundwater recharge is expected which poses a greater risk of groundwater and surface water drought in these catchments. Therefore, the simplistic approach for groundwater simulation presented in SWAT and HEC-RAS, which might give unreliable results in terms of the magnitude of groundwater flow contribution to the channel, is rather consistent with the overall direction presented by Bates et al. (2008).