Assessing Impacts of Land Use and Land Cover (LULC) Change on Stream Flow and Runoff in Rur Basin, Germany

Abstract

Understanding the impact of land use/land cover (LULC) change on hydrology is the key to sustainable water resource management. In this study, we used the Soil and Water Assessment Tool (SWAT) to evaluate the impact of LULC change on the runoff in the Rur basin, Germany. The SWAT model was calibrated against the observed data of stream flow and runoff at three sites (Stah, Linnich, and Monschau) between 2000 and 2010 and validated between 2011 and 2015. The performance of the hydrological model was assessed by using statistical parameters such as the coefficient of determination (R²), p-value, r-value, and percentage bias (PBAIS). Our analysis reveals that the average R² values for model calibration and validation were 0.68 and 0.67 (n = 3), respectively. The impacts of three change scenarios on stream runoff were assessed by replacing the partial forest with urban settlements, agricultural land, and grasslands compared to the 2006 LULC map. The SWAT model captured, overall, the spatio-temporal patterns and effects of LULC change on the stream runoffs despite the heterogeneous runoff responses related to the variable impacts of the different LULC. The results show that LULC change from deciduous forest to urban settlements, agricultural land, or grasslands increased the overall basin runoff by 43%, 14%, and 4%, respectively.

1. Introduction

An increasing population, climate change, and human activities have led to a progressive decrease in per capita water availability. Vörösmarty et al. [1] indicated that the changes in the world population and economic activities will influence the water supply demand relationship much more strongly than climate change over the period from 2000 to 2025. It is estimated that water use has increased 6-fold over the past century and is rising by 1% per year [2]. Anthropogenic activities such as deforestation, agricultural development, and urban development lead to change in LULC, which causes environmental change and has major implications on earth–atmosphere interactions and sustainable development [3,4]. Globally, humans are altering the landscape at an unprecedented rate and many regions have experienced a massive LULC change. However, there is a significant regional variability in the LULC change trends [5]. For instance, there is a decreasing trend in the forested area in the tropical region, while an increasing trend is observed in the temperate region [6,7]. During the 1980–2000 period, more than 55% of the new agricultural land across the tropics came at the expense of intact forests, and another 28% came from disturbed forests [8]. In the period 2010–2015, the forest area declined by 5.5 million hectare per year across the tropics and increased by 2.2 million hectare per year in the temperate region [9].

LULC changes have impacts on the hydrological processes (evapotranspiration, runoff, leaching, infiltration, etc.) and the hydrologic response [10]. Hence, the reliable quantitative estimation and calculation of various hydrological variables such as stream flows, evapotranspiration (ET), ground water level, and spatio-temporal fluctuations under variable basin conditions are very important for water resource management. Many studies have been conducted to quantify the impact of LULC change on stream flow. Zhang et al. [11] reviewed 312 basins worldwide. Their results showed that deforestation due to various reasons (e.g., urbanization, agricultural expansion, wildfires, and insect infestation) caused an increase in annual stream flow for more than 80% of the basins, whereas afforestation had relatively inconsistent effect on the stream flow. Generally, stream flow increases with deforestation whereas it decreases with afforestation [12,13,14,15,16,17]. However, the stream flow response to LULC change can differ significantly from the general trend especially in large basins (>1000 km²), and hence requires closer scrutiny [18,19]. In addition to the LULC, climate change also needs to be accounted for in the hydrological models to understand the hydrological regime of the region.

Key research methodologies to understand LULC change effects on hydrological processes include paired catchment experiments and hydrological modelling. In paired catchment experiments, two similar catchments are chosen for comparison (generally very close to each other). One is made the control catchment and the other various treatments such as afforestation, deforestation, and vegetation change are applied. Brown et al. [20] performed a review of the use of paired catchment studies for determining the changes in water yield due to variation in vegetation at various timescales. They found that the period taken to reach a new balance under permanent land-use change varies considerably. Deforestation treatments reach a new balance more quickly than afforestation ones. Bosch and Hewlett [12] reviewed 94 catchment experiments that monitored the effect of vegetation change on water resources. They found that in most of the cases, water yield is increased as forest cover is reduced, but in some occasional cases of the sparsely vegetated land, an increased forest cover reduces the water yield. It should be noted that there is high variability in the hydrologic response to such deforestation. Sahin and Hall [21] performed a meta-analysis of 145 catchment experiments using fuzzy linear regression to estimate the change in water yield due to the change in different vegetation types such as coniferous and hardwood. They found that a 10% reduction in the cover of conifer-type forest led to an increase in the water yield by about 20–25 mm, while that in eucalyptus type forest increased by only 6 mm. Brath et al. [22] showed that LULC change can affect the peak flow regime and flood frequency in the Samoggia River basin, Italy. Now, hydrological observatories have been set up to study, in detail, the causes of human influence on runoff volumes and water qualities in some basins, such as the Rur basin [23]. However, long-term catchment experiments, although an important method, take more time and are expensive to set up. Particularly, runoff behavior takes a long time to adapt to the vegetation change in catchment experiments [24].

Hydrological models offer a more viable option for quantifying water balance and evaluating land management practices, and play an important role in scaling up field observations and predicting the impacts of LULC changes on basin hydrology in the river basins [25]. Many basin models have been developed in the last three decades, such as Chemicals, Runoff, and Erosion from Agricultural Management Systems (CREAMS) [26]; Erosion Productivity Impact Calculator (EPIC) [27,28]; Agricultural Non-Point Source model (AGNPS) [29]; Soil and Water Assessment tool (SWAT) [30,31,32]; and Hydrologic Simulation Program–Fortran (HSPF) [33]. These basin models are applied for modelling runoff and soil loss prediction, water quality, land-use change effect assessment, and climate change impacts.

Among the foregoing models, the SWAT is physically a semi-distributed model for analyzing the impact of land management practices on water, sediment, and agricultural chemical yields in basins [30,34]. Numerous studies have been conducted using SWAT to evaluate the effect of LULC and basin management scenarios on various aspects of the hydrological cycle. Nie et al. [35] used SWAT to quantify the effect of LULC changes on hydrological components at the sub-basin scale in the upper San Pedro basin. Their results showed that the urbanization was the strongest contributor to the increase in surface runoff and water yield in the study period. In contrast, the replacement of desert scrub/grassland with mesquite was the strongest contributor to the decreased base flow/percolation and contributed to the increased ET. The SWAT model was used for simulating the stream and sediment runoff of Beressa basin in Ethiopia [36]. Their results showed that the basin experienced a significant increase in the stream and sediment runoff due to an increase in farmland and settlement areas. SWAT have been used for simulating stream flow [37], phosphorus export [38], pollutants [39,40], microorganisms [41], trends in frequencies of extreme precipitation and floods [42], greenhouse gas emissions [43,44,45], and water quality [32,46,47] in several catchments in Europe. Several studies showed that SWAT is not very sensible for land-use change, besides ET, simply because the corresponding change in soil parameters is not captured adequately [48,49,50].

Although the SWAT has been used for the effects of LULC changes on sediments, it has not considered the changes in soil properties and slope when LULC changes in the SWAT. The runoff response to the LULC change is a complex function of the slope, soil, vegetation, sensitivity of the landscape to runoff, basin management operations, etc. The Meuse River basin, a highly urbanized basin, is experiencing changes in its hydrological regime due to climate change and LULC change impacts [51]. The Rur River, a tributary of the Meuse River, supports a population of over 5 million by securing a water supply [52]. Despite the importance of LULC changes to the basins, no studies have been reported to quantitatively compare the performance of the widely used SWAT in the Rur basin. Particularly, LULC change could change other hydrological conditions of a basin, such as slope, ponds, and Riparian zones. These changes have not been represented explicitly in the hydrological response units of the SWAT when LULC changes. The influence of LULC change on runoffs can be masked due to multiple sources of uncertainties in the hydrological modelling in a large basin. Therefore, a careful examination in estimating the impacts of LULC changes is needed for managing water resources and evaluating the impacts of LULC changes on water resources. Under some extreme scenarios of LULC changes, the impacts of LULC changes on stream flow and runoffs of the Rur River need to be evaluated. Such a study has not been accomplished before in the Rur basin, and this area is characterized by strong land-use change, e.g., urbanization, reforestation, opencast mines, recultivation, etc. Finally, the Rur basin can be seen as representative of other areas with similar developments.

The main objective of this study is to evaluate the impacts of LULC changes on the stream flow of the Rur River and assess the feasibility of the SWAT model for this basin by performing sensitivity analysis, calibration, and validation. We will design several specific scenarios of LULC changes to represent extreme LULC changes, primarily concerning deforestation, urbanization, and agricultural expansion. We will qualify the response of stream flow and runoff to different LULC change scenarios. Hence, this study will be helpful for the basin managers and will open up new dimensions for various other applications of SWAT in the Rur basin.

2. Materials and Methods

2.1. Study Area

The Rur basin, covering 2354 km², is located largely in North Rhine-Westphalia (Germany) with small parts also in Belgium (6.7%) and the Netherlands (4.6%) [23]. It lies between the coordinate locations around 6° E–6°40′ N E and 50°10′–50°30′ N (Figure 1) and shows variability in topography, land use, soil type, and meteorological parameters, which are more pronounced when comparing northern and southern parts. The Rur basin and its sub-basins are shown in Figure S1. The elevation and mean annual temperature of the region approximately decline from 680 m to 30 m and 8.5–10.5 °C to 7.0–9.0 °C, respectively, as we move from north to south. Further, the study area is characterized by flat lowland in the northern part which is part of the German–Belgian loess belt consisting of unconsolidated rock deposits of mainly tertiary and quaternary periods. In the northern part, the major soils are Cumulic Anthrosols near the drainage lines and Haplic Luvisols, both with silt loam textures [53]. Soils with a loamy sand texture (Fimic Anthrosols and Dystric Cambisols) are located in the northernmost part of the loess plain (Figure S2). The soil in this region is very productive. Consequently, the dominant land use in the northern part is cropland, with the main crops being winter wheat, sugar beet, and maize (Figure S3). Other important land use includes open-cast lignite mines and urban settlements. This lowland region receives between 650 mm and 900 mm of precipitation annually, of which between 500 and 600 mm/year evaporates. The southern part of the Rur basin belongs to the low mountain range Eifel, which mainly consists of consolidated rocks of the Rhenish Massif from the Paleozoic and Mesozoic. The hard rock terrain, combined with mountainous slopes and very shallow soil layers, causes low infiltration rates, so that the area is at risk of flooding. Further, low temperature in winters cause low ET. Consequently, strong seasonality with a strong runoff response to rainfall events is observed [54]. Mean annual precipitation in the southern region ranges from 1200 mm/yr (windward side of High Fen Mountains) to 700 mm/yr (leeward side). Coniferous forests, deciduous forests, and pastures are dominant parts of the southern region [23]. Coniferous forests (17.6%), deciduous forests (23.8%), and pastures (15.6%) are dominant parts of the southern region. The main crops are winter wheat (WWHT) (13%), sugar beet (SGBT) (12.6%), and maize (CORN) (5.9%). Other important land uses include open-cast lignite mines and urban settlements (5.6%).

2.2. SWAT Model Description

SWAT is a small watershed to river basin-scale model that considers both upland and stream processes and was developed originally by the United States Department of Agricultural Research Service [2,55,56]. In SWAT, a basin is divided into multiple sub-watersheds or sub-basins to incorporate spatial heterogeneity in terms of land use land cover, soil type, and topography. Each sub-basin is divided into areas of homogenous land-use, soil, and slope; these are called hydrologic response units (HRUs) [57]. SWAT can simulate a number of different physical processes, such as snowfall and melt, vadose zone processes (infiltration, evaporation, plant uptake, lateral flows, and percolation), and groundwater flows in a basin [57]. The stream flow and runoffs are predicted separately for each HRU and routed, using a variable channel routing method [58,59], to obtain the total stream flow and runoff at the outlet of a basin. The SWAT model uses the Natural Resource Conservation Service (NRCS) curve number method [60] to relate a calculated runoff curve number (CN) to direct runoff, after taking into account initial abstraction losses and infiltration rates for estimating surface runoff (Q_surf). The fundamental hydrology of a basin in SWAT is based on the following water balance Equation (1), based on mass balance, which calculates the change in soil water content (SW_t):

$$S{W}_t=S{W}_0+{{\displaystyle \sum }}_{i=1}^t\left(R_{day}-Q_{surf}-E_a{W}_{sweep}-W_{gw}\right)$$

(1)

SW_t = final water content in (mm);

SW₀ = initial water content in (mm);

t = time in (days);

R_day = precipitation amount on specific days i (mm);

Q_surf = runoff amount on specific days i (mm);

E_a = evapotranspiration amount on day i (mm);

W_sweep = the amount of water percolated into the vadose zones on a day i (mm);

W_gw = return amount of flow on a day i (mm).

Soil Conservation Service (SCS)–Curve Number (CN) approach is an empirical model developed to estimate the rainfall–runoff relationships based on the land use, soil type, and soil hydrologic conditions and slope. CN values range from 0 to 100. A higher CN (e.g., 100 for pavement) means a higher potential of the surface to create runoff, in contrast to lower CN values (e.g., near 0 values for dry soils on very less slop terrain). The equation used for estimating the surface runoff is as follows:

$$Q_{surf}=\frac{{\left(R_{day}-I_a\right)}^2}{R_{day}-I_a+S}$$

(2)

Q_surf is the surface runoff. $I_a$ is the initial abstraction and S is the retention parameter. The retention parameter varies spatially due to changes in soils and land use; changes in soil water content are defined as

$$S=25.4\left(\frac{1000}{CN}-10\right)$$

(3)

The retention parameter is defined as (2) where CN = curve number for the day. The initial abstraction $I_a$ is commonly approximated as 0.2S. Therefore, Equation (2) becomes

$$Q_{surf}=\frac{{\left(R_{day}-0.2S\right)}^2}{R_{day}-0.8S}$$

(4)

The runoff will only occur when $R_{day}$ > $I_a$. The graphical solution of Equation (3) for different curve number values is also well documented in [60].

2.3. Model Input and Setup

Basin delineation is used to divide the basin into smaller sub-basins, called sub-basins with the stream reaches and directions. The digital elevation model is used when the basin delineation and flow directions (upstream and downstream) are decided. Sub-basin numbers are generated and are used as references throughout the SWAT simulation. On the stream network, sub-basin outlets are created, where water accumulates in the sub-basin to become stream flow. In this study, 27 sub-basins were created. The stream network for the Rur basin was delineated using the digital elevation model in the ArcSWAT 2012 software. Figure 2 shows a flowchart of the SWAT model’s implementation.

2.3.1. Land Use and Soil Definition

Land Use definition is performed where the SWAT-supported land-use classes are mapped to the prepared land-use map of the area. In this study, land use data were downloaded from the CORINE land cover data (CLC) website for the year 2006 and categorized into 13 classes (Figure S3) using the ArcSWAT software.

The SWAT model requires different soil textural and physicochemical properties such as soil texture, available water content, hydraulic conductivity, bulk density, and organic carbon content for different layers of each soil type (Neitsch et al., 2011) [61]. The soil data were downloaded from the Food and Agricultural Organization (FAO), US website (Table S1). The soil raster file for the study area was clipped from the global soil raster database downloaded from the FAO website. For soil definition, the soil raster database was downloaded from the Food and Agricultural Organization (FAO) website and a soil map for the area was prepared using the SWAT-supported soil classification. SWAT uses the NRCS CN method to estimate the runoff. Curve numbers are calculated for each combination of land use and soil.

2.3.2. Digital Elevation Model (DEM), Slope, and Hydrological Response Units (HRUs)

The spatially distributed data (GIS input) needed for the Arc SWAT interface include the Digital Elevation Model (DEM), soil data, land use, and stream network layers. Topography was defined by a DEM that describes the elevation of any point in each given area at a specific spatial resolution. The DEM was used to delineate the basin and estimate the sub-basin parameters, such as the slope gradient, slope length of the terrain, and the stream network characteristics such as channel slope, length, and width. Five slope classes were made for this study, as shown in Figure S4. These classes were used in the HRU definition. The slope is an important parameter in estimating the flow rates in the reaches.

The HRUs are the unique combinations of land use, soil, and slope. After the land-use, soil, and slope definition in the SWAT, they are overlaid to create HRUs to account for the spatial heterogeneity of the basin. However, the level of complexity can be adjusted by providing thresholds (in percentage) for the land use, soil, or slope during HRU definition. The greater the value of the threshold, the lesser the number of HRUs. In this study, thresholds of 0%, 5%, and 10% were provided for the land use, soil, and slope, respectively. A total of 758 HRUs were generated in this study. Creating HRUs increases accuracy and complexity by providing a more detailed description of the water balance. In these circumstances, the number of HRUs can be reduced by choosing appropriate thresholds for land use, soil, and slope during HRU definition [61]. Runoff calculations are performed separately for each HRU and then routed to obtain the total runoff for the basin.

2.3.3. Meteorological Data

Weather data include precipitation, temperature, relative humidity, solar radiation, and wind speed. Of these, precipitation, temperature, and relative humidity are the most important data for hydrological modelling. The weather input data are prepared using observed weather data (temperature, relative humidity, and precipitation) in a specific format supported by the SWAT. The locations of three weather stations are listed in Table S2. SWAT provides a weather generator that uses the average weather statistics for interpolations of the missing data between weather stations. For this study, a daily time, and a warm-up period of 6 years were chosen. The simulation was run, and the different output files for reaches, sub-basins, and HRUs were written. SWAT computes hydrological components such as runoff, ET, sediments, nutrients, and other parameters that can be obtained from the TxInOut folder. Weather data and river discharge data were also used for calibration purposes and for the prediction of stream flow. Table S3 shows the sources of different input data with their description.

2.4. Model Sensitivity Analysis, Calibration and Validation, and Performance

2.4.1. Model Sensitivity Analysis

There are two types of sensitivity analysis ‘one at a time’ and ‘global’. In “one at a time” sensitivity analysis, one parameter is selected at a time and its sensitivity is assessed. Spruill et al. (2000) [62] gave a detailed description of the sensitivity analysis results for a ‘one at a time’ selected parameter for daily stream flow, for a catchment in Kentucky, USA. However, the sensitivity of a parameter often depends on the value of some other parameter because the parameters are not necessarily independent. In this study, we have used global sensitivity analysis, analyzing all the parameters for their sensitivity simultaneously. Model parameters depict the catchment properties and the processes in the catchment. A total of 18 model parameters were identified to affect surface runoff from the preliminary literature survey. For instance, Kannan et al. [63] and White et al. [64] described sensitive parameters for surface runoff. These parameters were checked for sensitivity, and, finally, 13 parameters were selected.

Sensitivity is assessed by measuring the change in the output variable with respect to the change in the input parameter. A greater change in the output variable corresponds a to greater sensitivity. Model sensitivity analysis is an important step during the calibration phase for choosing the most sensitive parameters from a bigger pool of potential parameters affecting the output variable. For sensitivity analysis, SWAT-CUP 2012 software was used, incorporating sequential uncertainty fitting (SUFI-2) optimization algorithm, whereby an objective function (g) was chosen (R² in this study) and initial uncertainty ranges were set for each parameter (in Equation (5)). The upper and lower limit defining the uncertainty band of the parameter was a subset of the absolute ranges provided by SWAT. SUFI-2 uses Latin hypercube sampling [65] to generate as many combinations of parameter values as there are simulations (100 simulations in our case). These parameter values were regressed against the chosen objective function (equation) and a t-test was applied for assessing parameter sensitivities and used for calibration and validation processes.

$$b_{j,min}\le b_j\le b_{j,max}$$

(5)

where j = 1, m.

$$\mathrm{g}=\mathsf{\alpha }+{{\displaystyle \sum }}_{i=1}^m\left({\beta }_i{b}_i\right)$$

(6)

where m is the number of parameters.

2.4.2. Model Calibration and Validation

Calibration is the process where the parameter ranges are adjusted in such a way that a good fit is obtained between the observed and simulated time series data for the period chosen as the calibration period. A different time period is chosen as a validation period (having no overlap with the calibration period), with the same parameter values and the number of simulations used as in the final calibrated model. The goodness-of-fit is assessed based on the criteria that maximum observed data points are bracketed in the 95% prediction uncertainty (95 PPU) band (signifying uncertainty in the predicted values) [66]. The fraction of observed data lying in the 95 PPU band is denoted by the p-factor (value of p = 1 being ideal), whereas the thickness of the 95 PPU band (denoted by r-factor) should be as small as possible (value of r close to 0 being ideal) [66]. In the present study, R², p value, r value, and percentage bias (PBIAS) are used as performance statistics.

After obtaining and processing input data, they were used to calibrate the SWAT model from 2000 to 2010 and then to validate the model from 2011 to 2015. The LULC used in the calibration and validation is taken from CLC 2006. The model calibration procedure is developed based on optimization techniques [67] with the assumption that an optimal set of parameters exists in the model to describe the surface runoff in the Rur basin. The results from these permutations suggested that the set of optimal values of model parameters allowed the model to optimally describe the basin hydrology, in the sense of having the least error between the simulated and the observed surface runoff. Data from three stations were used for calibration and validation following the sequential calibration technique.

2.5. Model Performance Evaluation

The performance of the hydrological model was assessed by using statistical parameters such as the coefficient of determination ($R^2$), p value, r value, and PBIAS.

2.5.1. Coefficient of Determination (R²)

The degree of correlation between the observed and the simulated values is measured by the $R^2$ value. $R^2$ ranges from 0 to 1. An $R^2$ close to 1 means that observed and simulated data points are in good agreement, whereas a value close to 0 means no agreement. Although $R^2$ has been widely used in the literature for estimating the “goodness-of-fit” of the hydrological models, it is recommended to supplement the $R^2$ with other performance statistics because of its sensitivity to outliers and its inability to discern the additive and proportional differences between simulated and observed data points [68]. The expression for $R^2$ is

$$R^2={\left[\frac{{{\displaystyle \sum }}_{i=1}^n\left(O_i-O_{avr}\right)\left(P_i-P_{avr}\right)}{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-O_{avr}\right)}^2{{\displaystyle \sum }}_{n=1}^{\infty }{\left(P_i-P_{avr}\right)}^2}\right]}^2$$

(7)

where $R^2$ is the coefficient of determination, $O_i$ the ith observed value, $O_{avr}$ is the average observed value of the entire period in consideration, $P_i$ the ith simulated (modeled) value, $P_{avr}$ is the average of the simulated values of the entire study period.

2.5.2. p Value and r Value

The goodness-of-fit is assessed based on the criteria that maximum observed data points are bracketed in the 95% prediction uncertainty (95 PPU), calculated at the 2.5 and 97.5 percentiles of the cumulative distribution of the simulated variables to represent the uncertainty in the predicted values. The fraction of observed data lying in the 95 PPU band is denoted by the p-factor (value of p = 1 being ideal). Further, the thickness of the 95 PPU band (denoted by r-factor) should be as small as possible (value of r close to 0 being ideal). A detailed description of p and r values can be found in [66,69].

2.5.3. PBIAS Value

Percentage bias (PBIAS) measures the average tendency of the simulated data to be underestimated or overestimated. During calibration, one tries to achieve lower values of PBIAS. A positive PBIAS value means underestimation whereas a negative PBIAS means overestimation. A value of PBIAS within the range of ±25% is considered to be satisfactory for the runoff simulations [70]. The mathematical expression for the PBIAS value is

$$\mathrm{PBIAS}=100\times \frac{{{\displaystyle \sum }}_{i=1}^n\left(Q_{m,i}-Q_{s,i}\right)}{{{\displaystyle \sum }}_{i=1}^n\left(Q_{m,i}\right)}$$

(8)

where Q stands for the hydrological variable, m stands for measured, and s stands for simulated. The above expression is minimized if PBIAS is chosen as the objective function in the SWAT.

For successful calibration and model performance assessment, Moriasi et al. [70] suggested performance statistics such as R², NS, PBIAS, RSR, etc., and their ranges. However, it should be noted that these ranges are for a monthly time step. Generally, the values of the test statistics deteriorate as the time resolution increases to a daily time step. In this study, simulation was performed at the daily time step and the ranges in [70] are used only as the reference.

2.6. Scenario Generation and Impacts

Three hypothetical scenarios were constructed to assess the impact of the LULC on stream runoff due to forest area change. The average value of percentage change in the daily runoff was calculated using the following formula:

$$Q_{average of daily percentage change}=\frac{{{\displaystyle \sum }}_{i=1}^n\frac{\left(Q_{scenario, sim, i}-Q_{base level, sim, i}\right)\times 100}{Q_{base level, sim, i}}}{n}$$

(9)

Variability in the percentage change in the daily runoff values was calculated by taking the standard deviation of the percentage change in the daily runoff using the formula

$$Q_{variability }=Stdev \left(Q_{percentage change in runoff for day i}\right) 1\le i \le n$$

(10)

$$Q_{percentage change in runoff for ith day}=\frac{\left(Q_{scenario, sim, i}-Q_{base level, sim, i}\right)\times 100}{Q_{base level, sim, i}}$$

(11)

The percentage change in the average long-term (16 years) runoff was calculated using the formula

$$Q_{percentage change in average runoff}=\frac{Q_{average, scenario,sim, i}-Q_{average,base, sim, i}}{Q_{average,base, sim, i}}\times 100$$

(12)

$$Q_{average, scenario,sim, i}=\frac{{{\displaystyle \sum }}_{i=1}^n{Q}_{scenario, sim, i}}{\mathrm{n}}$$

(13)

$$Q_{average, base level,sim, i}=\frac{{{\displaystyle \sum }}_{i=1}^n{Q}_{base level, sim, i}}{\mathrm{n}}$$

(14)

3. Results and Discussion

3.1. Sensitivity Analysis

To identify sensitive parameters, 18 parameters were chosen that affected the surface runoff, reviewing previously calibrated SWAT parameters and documentation from the SWAT manual [61]. Global sensitivity analysis was performed using SWAT-CUP to identify the 13 most sensitive parameters. In this method, the smaller the p-value and greater the absolute value of t-stat, the more sensitive the parameter [71]. Table S4 shows the parameters used for different stations and their ranks at the end of the calibration process for (a) Monschau, (b) Linnich, and (c) Stah, respectively. In order to calibrate the model, the parameter ranges were adjusted to minimize the difference between simulated and observed surface runoff values until an acceptable agreement was attained.

Groundwater delay (GW_delay), curve number (CN2), and base flow factor (ALPHA_BF) were the three most sensitive parameters for all three stations. GW_delay is the time lag between water exiting the soil profile and recharging the aquifer. It depends on the depth of the water table and the hydraulic properties of vadose and ground water zones [61]. Sangrey et al. [72] provides a mathematical relationship between shallow aquifer recharge and GW_delay. Soil evaporation compensation factor (ESCO), average slope length (SLS_SUBBSN), and HRU_SLP (average slope steepness) are the lower sensitivity parameters.

3.2. Model Calibration and Validation

Niraula et al. [73] concluded that the results for the relative change in the stream runoff for the un-calibrated and the single-outlet-calibrated models can vary significantly from the spatially-calibrated models. Hence, spatially-calibrated models should be used to assess the relative change in the stream runoff due to land-use change. Model calibration and validation were carried out for three stations, Stah, Linnich, and Monschau, with the selected parameters after sensitivity analysis using a sequential uncertainty fitting (SUFI-2) algorithm in SWAT-CUP [66]. The model was calibrated for a period of 11 years (2000 to 2010) with a 6-year warm-up period (1994 to 1999) at a daily time step. The model was recalibrated multiple times by changing the lower and upper bounds of the parameter so as to obtain an acceptable agreement between the simulated and observed data. The model was automatically iterated 100 times to obtain the final simulation for the calibration period. The results are given in

Table 1. Performance statistics for the simulation at various stations.

Station	Year	Period	Evaluation of Statistics
			p	r	R2	PBIAS
Monschau (sub-basin 21)	2000–2010	Calibration	0.81	0.84	0.60	11.3
	2011–2015	Validation	0.78	0.52	0.67	30.5
Linnich (sub-basin 6)	2000–2010	Calibration	0.67	1.42	0.71	−4.4
	2011–2015	Validation	0.63	1.12	0.66	−17.3
Stah (sub basin 1)	2000–2010	Calibration	0.53	1.8	0.73	−8.6
	2011–2015	Validation	0.48	1.49	0.69	−29.2

.

Figure 3 shows the comparisons of runoffs between simulation results and observed data for the calibration period (2000–2010) and validation period (2011–2015) at all three stations: (a) Monschau, (b) Linnich, and (c) Stah. The model’s performance during calibration and validation was the best at the Monschau station, followed by Linnich and Stah stations. This is evident as the statistical performances are satisfied from the performance statistics (R², p, and r value) in Figure 4 and Table 1. The key reason for this can be a relative increment in the anthropogenic interventions as we move from the Monschau station in sub-basin 21 to the Stah station in sub-basin 1. From the land-use map, it is evident that the northern part of the basin has much more agricultural and urban land than the southern part of the basin. Water use by the urban population, industries, channel diversions, and dams are not considered in this model but can significantly affect the stream flow in the northern parts of the basin. This can explain why the Monschau station has better performance statistics than the other two stations in regard to both calibrated and validated values. Further, one land-use map was used to calibrate the model for all the years in which the scenario was not a real case. Our simulation was performed at a daily time step, and all absolute values of PBIAS in calibration periods were less than 25%, thus being satisfactory (Table 1). Moriasi et al. [70,74] suggested that an absolute value of PBIAS of up to 25% was satisfactory at a monthly time step. At the basin scale, the result was consistent with other studies (runoff was simulated and compared with the base-level scenario simulation for the individual gauging stations also).

3.3. Effect of LULC Change Scenarios on Stream Flows and Runoffs

Three hypothetical scenarios were constructed to assess the impact of the LULC on stream runoff due to forest area change. Most of these forests were present in parts of the Rur catchment lying in the regions of the western highlands of Germany. Deciduous forests comprise 57.39% of the total forests (coniferous and deciduous combined). For the sake of simplification, the entire deciduous lands in the three scenarios were assumed to be converted to the use of urban settlements, agriculture, or grasslands, respectively (Figure 5). Table S5 shows the description of the three scenarios. The three scenarios were compared with the base-level scenario that was based on the 2006 land-use map. Figure S5 shows that ET presents no changes across any of the scenarios, but there are clear differences regarding percolation, return flow, and runoff between the scenarios. The frsd_urbn scenario has the biggest runoffs, followed by the frsd_agrl scenario and the frsd_rnge scenario.

3.3.1. Forest to Urban Residential

In the frsd_urbn scenario, the simulation results showed that the long-term average runoff increased by 43.05% from 94.48 mm in the base-level scenario to 135.15 mm in the frsd_urbn scenario in the entire basin (Figure 6). This was mainly due to the decrease in percolation from 249.83 mm to 210.61 mm after scenario application. The results also showed that the average runoff in the entire basin increased by 13.47% for the frsd_agrl scenario and 3.45% for the frsd_rnge scenario (Figure 6). This was primarily attributed to decreased percolation in all the scenarios. The decrease in percolation was highest for the frsd_urbn scenario because the impervious surface increased after scenario application. However, the average yearly runoff values decreased for the frsd_urbn scenario at Monschau and Stah while they increased at Linnich station.

Figure 7 shows the time series plot of the relative change (with respect to the base level) in the cumulative stream flow for the period from 2000 to 2015. It was found that there was pronounced seasonality in this scenario. The relative change in the cumulative flow showed a decreasing trend (Figure 7a). This signifies that less water for this scenario will be available from stream flow than that for the base-level scenario at Monschau station. However, at Linnich, the change increased first and then decreased (Figure 7b). There was a peak when the linear trend line fitted in the time series dataset. However, the change showed lower seasonality compared to that at Monschau station. This means that the upstream was more sensitive to seasonality than the downstream. This was consistent with the results in

Table 2. Station-specific effects of various scenarios. (% changes used are with respect to the base-level scenario).

Station	Scenario	Average Daily Runoff Change (%)	Average Daily Runoff Change (Stdev) (%)	Average Long-Term (16 Years) Runoff Change (%)
Frsd_urbn	Monschau (sub-basin 21)	7.42	56.20	−0.36
	Linnich (sub-basin 6)	−6.47	29.87	2.17
	Stah (sub-basin 1)	−9.42	27.30	−8.97
Frsd_agrl	Monschau (sub-basin 21)	1.74	7.08	0.81
	Linnich (sub-basin 6)	0.55	1.95	0.54
	Stah (sub-basin 1)	1.18	7.64	0.46
Frsd_rnge	Monschau (sub-basin 21)	−1.93	4.84	−0.72
	Linnich (sub-basin 6)	−0.65	1.39	−0.60
	Stah (sub-basin 1)	−0.10	2.86	−0.45

. In Table 2, the average percentage change in the daily runoff was a positive value of 7.42% at Monschau station while it was a negative value of −6.47% at Linnich station and −9.42% at Stah station. This means that the average daily runoff increased at Monschau station and decreased at the other two stations compared to their baseline level. To obtain a sense of variability in the average daily runoff values, the standard deviation for the daily relative change was calculated to be 56% at Monschau station (%) but 29.87% at Linnich station. This meant that there was a lot of variability observed in the daily runoff values at Monschau station. This can be explained by the fact that when the forests were removed, there was a relatively quick runoff response to rainfalls at a daily time step. The high value of variability can be explained due to the location of the Monschau station in the forested area, with steeper slopes and shallow soil layers.

On the other hand, the annual surface runoff values showed a mixed trend. The average long-term (16 years) runoff change decreased by −0.36% at Monschau station and increased by 2.17% at Linnich station. This is because Linnich station is located in the plane region of the Rur basin, with a slope mostly within the 0–5 % range (Figure S4). The long-term average of the runoff at Stah station became more negative than at Linnich with a value of ~8.97%. The percentage changes in the average monthly and average yearly values were also computed (Figure 8). Percentage change in the monthly values averaged over the period of 2000–2015 showed a positive change in the period April- September compared to the baseline level, while a negative change was observed for the period October–March. The average yearly values showed a mixed trend but a consistent negative relative change from 2011 onwards, although most of these changes were only marginal except for in a few years. Several reasons are attributed to this observation. Firstly, the basin is a complicated, integrated system and the runoff has a different sensitivity to LULC change for different sub-basins. The runoff trend at a certain sub-basin can decrease even if urbanization occurs at the cost of deforestation. There can be a masking effect of some sub-basins [11].

To delve a bit into this issue, the Mann–Kendall test was applied to these three stations. The Mann–Kendall test, as a rank-based non-parametric test, is widely used to test for statistically significant increasing or decreasing trends in the data series [75,76]. The n time series values (Y₁, Y₂, Y₃, …, Y_n) are replaced by their relative ranks (S₁, S₂, S₃, …, S_n); the statistical value P of Mann–Kendall test is

$$P={{\displaystyle \sum }}_{1-1}^{n-1}\left[{{\displaystyle \sum }}_{j=i+1}^{n}sgn \left(S_j-S_i\right)\right]$$

(15)

where sgn(y) = 1 for y > 0;

sgn(y) = 0 for y = 0;

sgn(y) = −1 for y < 0.

If the null hypothesis H_o is true, then Z is approximately normally distributed with

μ = 0;

σ = n (n − 1) (2n + 5)/18.

The z-statistic is therefore

Z = P/${\sigma }^{0.5}$.

A positive value of p signifies that there is an increasing trend and vice versa [77]. At a significance level of 0.05, if

Z < − 2, then the time series is said to have a statistically significant decreasing trend.

If Z > 2, then the time series is said to have a statistically significant increasing trend.

The z-statistics were calculated (Table S6). The observed values of runoff for the Stah station showed a clear decreasing trend from the years 2007 to 2015 (Figure 9). This was because there was a significant increase in the urban fabric and a decrease in the deciduous forest when comparing the LULC map with the CLC data of 2012 and 2000.

3.3.2. Forest to Agriculture Land Conversion

In the frsd-agrl scenario, the deciduous forest was converted to agricultural land. Thus, agricultural land comprises ~34% of the total basin, with winter wheat, sugar beet, potato, corn, and barley being the dominant crops. The results showed that the long-term average of runoff increased by 13.47% from 94.48 mm in the base-level scenario to 107.21 mm in the frsd_agrl scenario for the entire basin (Figure 6). This was mainly due to the decrease in percolation after scenario application.

Figure 10 shows that the simulated cumulative stream flows increased linearly with little seasonality change at Monschau station and Linnich station. The percentage change of the average daily runoff value was positive for all three stations (Table 2). However, the variability in the relative change of daily runoff was less than that of the frsd_urbn scenario. This implied that the daily runoff was more sensitive to the frsd_urbn scenario than the frsd_agrl scenario in the basin. The long-term average runoff showed an increment at all the three stations (Table 2).

Figure 11 shows the percentage changes in the average monthly and average yearly values that confirm the results of Table 2. The percentage change in the monthly values averaged over the period of from 2000 to 2015 indicated that most of the average monthly runoffs showed either a positive increase (June to September) or a marginal decrease (December to May) in spring and fall. Compared to the frsd_urbn scenario, the long-term averages in the runoff increased for the three scenarios for the entire basin, with exceptional margins across several years. This signifies that more water will be available from stream flow for the frsd_agrl scenario than the base-level scenario. The increase was primarily due to the greater runoff generation of the agricultural land, attributed to the relatively high SCS curve number of the agricultural land as compared to the deciduous forest.

3.3.3. Forest to Perennial Grassland Conversion

In the frsd-rnge scenario, the deciduous forest was converted to perennial grasslands. Grasslands comprise ~16% of the total basin. The results showed that the percentage change in the runoff slightly increased by 3.45% from 94.48 mm in the base-level scenario to 97.74 mm for the entire basin (Figure 6). In Table 2, the percentage change in the average daily runoff was a negative value at all three stations, while the standard deviation in the daily runoff value was positive. The simulated cumulative stream flows decreased over time at the Monschau and Linnich stations compared to the base-level scenario (Figure 12). The long-term average values of the runoff showed a marginal decrease at all three stations. The percentage changes in the average monthly runoff were mostly negative, except for several months in winter at Stah station (Figure 13). The average yearly values showed a consistent decrease for almost all of the years. This was mainly due to the increase in percolation after scenario application compared to the other two scenarios (Figure S5). This implied that the frsd_rnge scenario could reduce the runoff compared to the base level.

3.4. Discussion

The hydrological response to the LULC change is a complex function of the slope, soil, vegetation, sensitivity of the landscape to runoff, basin management operations, etc. The influence of LULC change on runoff can be masked by the complexities found in a large basin [78,79]. Normally, LULC change should not be specific land use or cover. It may involve different land-use changes occurring in mutual ways. For example, for urban growth, the land-use changes may include changes from forest to urban, crop area to urban, and grassland to urban. On the other hand, the greenbelt increase is a land change from buildings to plants. Therefore, the specific changes have several reasons: (1). It is difficult to reflect the reality in the mathematical representation because this needs complex processing in the case of spatial and temporal changes. Thus, LULC changes have to be simplified more or less in any modelling studies. (2). If different LULC changes are considered at the same time, as a result, it will be difficult to explain which land-use change is more important than the other land-use changes. Thus, we think it is better to define a single LULC change to be able to compare the impacts of the different LULC scenarios. (3). The scenarios could be designed for some extreme conditions, such as climate change RCP8.5. Although they are not present in reality, they can somewhat represent the complexity.

The frsd_urbn scenario showed a strong increase in average daily discharge changes compared to the baseline level at the station located in the mountain range, but a decrease at the stations in the lowland agricultural region. Zhang et al. [17] also indicated that upstream forests play an important role in regulating the stream flow through canopy interception, evapotranspiration, and soil interception. In contrast, at the upstream station (Monschau), there were only marginal increases for the frsd_agrl scenario and decreases for the frsd_rnge scenario compared to the base level. The marginal changes in the average daily runoff for the frsd_agrl and the frsd_rnge scenarios can be due to effects of the percolation and return flows where the runoff response of deforestation was not decreased in a sub-basin. For the frsd_urbn scenario, the average long-term (16 years) runoff decreases at the upstream and downstream stations, but marginally increases compared to the base level. In contrast, the average long-term runoff increases for the frsd_agrl scenario and decreases for the frsd_rnge scenarios compared to the base level. These differences in the short- and long-term responses of runoff at the different stations to LULC changes highlight the importance of this study. For example, for the frsd_urbn scenario, a fast increased response of daily and monthly average runoffs implies that urbanization may cause regional flooding and seasonal drying. On the other hand, the grasslands are better able than the forests to prevent runoff and soil erosion. Therefore, the quantitative analysis for the three LULC change scenarios is helpful in developing the best land management practices. In addition, it is found that the runoff response to the LULC change can exhibit a threshold where the change in the surface runoff is not significant unless the threshold in the LULC change is reached. This is consistent with the results by Ghaffari et al. [80]. Further studies remain to be conducted to quantify the threshold for different LULC changes.

SWAT might show a better performance in some Literature [70,74]. It is still not easy for the model to fully capture high peaks and valleys. R² and r are oversensitive to high extreme values, and insensitive to additive and proportional differences between model prediction and observed data. Therefore, this depends on the dataset quality available rather than monthly or daily data. Generally, the calibration and validation of stream flow are better than runoff, base flow, or water quality. This is mainly because the runoff data, base flow data, and water quality are much less than the steam flows in most basins. In Table 8 of the recent reference [74], for the stream flow, the model performance at the basin scale can be evaluated as “good” if daily values are 0.70 ≤ R² ≤ 0.85 or “satisfactory” if daily values are 0.50 < R² < 0.70. For sediment, model performance at the basin scale can be evaluated as “good” if the daily values are 0.65 ≤ R² ≤ 0.80 or “satisfactory” if the daily values are 0.40 < R² < 0.65. Therefore, our results should be “good” since we used runoffs for the calibration and validation. In addition, it can be found that the performances in Figure 4 were much better than those in Figure 3. All R² in Figure 4 are greater than 0.79.

This study is limited in some respects. For instance, the climate change impact and water-use patterns were not accounted for. In a recent study, a climate change simulation and trend analysis of extreme precipitation and floods demonstrated a significant impact on the hydrology of the Rur basin in the mesoscale Rur catchment in western Germany until 2099, using the Statistical Downscaling Model (SDSM) and the SWAT model [42]. Furthermore, the SWAT model does not account for possible changes in soil properties and slopes related to LULC changes and agricultural management, such as grazing, water-use, tillage, and compaction. Particularly, these operations could significantly influence the potential maximum soil retention. In this regard, the current algorithm for the potential maximum soil retention has not represented the importance of vegetation growth and LULC change, and needs to be improved. In addition, the model could be enhanced by including dynamic land-use change, and more data from the ongoing TERENO project (Bogena et al., 2018) [23] could be used for more detailed model evaluation. Finally, the robustness and predictive power of the model can be improved by incorporating the reservoirs, dams, industrial and agricultural water-use data, and climate change scenarios in future studies.

4. Conclusions

Four LULC scenarios were used for assessing the impacts of the LULC changes on stream flow and runoff in Rur basin using the hydrological model SWAT. Three of them comprised deforestation, wherein the deforested land is converted to urban settlements, agriculture, and grasslands, respectively, and one was the baseline scenario. The SWAT model was calibrated for the period 2000–2010 and validated for the period 2011–2015 for all three runoff gauging stations (e.g., Monschau, Linnich, and Stah). The simulation model results showed that the hydrological responses to land use/land cover (LULC) changes in the Rur basin were very different, both for the individual stations (Monschau, Linnich, and Stah) and for different times of the year. The land-use change in the upstream areas can affect more stream flow and runoff than that in the downstream. The relative average daily runoffs were the most sensitive to the partial forest to urban (frsd_urbn) conversion scenario, particularly at the upstream station (Monschau). Therefore, rainfall events lead to higher runoffs or potential local flooding. The long-term average runoff for the entire basin increased by 43% from 95 mm in the base-level scenario to 135 mm in the frsd_urbn scenario due to the percolation decrease. The model also predicted an increase of 12 million cubic meters (m.c.m.) at Monschau and 50 m.c.m. at Linnich in the relative cumulative stream flows and runoffs for the scenario from partial forest converted to agricultural land (frsd_agrl), which may cause soil erosion. In contrast, for the partial forest conversion to rangeland (frsd_rnge), the average daily runoff and annual runoff decreased from −0.1% to −1.93% at all three stations compared to the base-level scenario. This suggested that the presence of more grasslands could reduce runoff and soil erosion. Critically, there are thresholds within the various type and the level LULC conversion wherein a particular sub-basin is not affected until some configuration of land-use change reaches the onset of the threshold. Finally, our study opens up new perspectives for qualifying the impacts of various LULC changes on runoffs in the Rur, and ultimately the Meuse basin, and it provides a powerful tool for runoff, sediment, and implicit nutrient management at the basin scale. Future work needs to account for the impacts of climate change and water-use patterns. The SWAT model should be improved to consider changes in the soil properties and slope when LULC changes in the SWAT.