**The Long-Term E**ff**ects of Land Use and Climate Changes on the Hydro-Morphology of the Reno River Catchment (Northern Italy)**

#### **Donatella Pavanelli 1,\*, Claudio Cavazza 2, Stevo Lavrni´c <sup>1</sup> and Attilio Toscano <sup>1</sup>**


Received: 3 July 2019; Accepted: 28 August 2019; Published: 3 September 2019

**Abstract:** Anthropogenic activities, and in particular land use/land cover (LULC) changes, have a considerable effect on rivers' flow rates and their morphologies. A representative example of those changes and resulting impacts on the fluvial environment is the Reno Mountain Basin (RMB), located in Northern Italy. Characterized by forest exploitation and agricultural production until World War II, today the RMB consists predominantly of meadows, forests and uncultivated land, as a result of agricultural land abandonment. This study focuses on the changes of the Reno river's morphology since the 1950s, with an objective of analyzing the factors that caused and influenced those changes. The factors considered were LULC changes, the Reno river flow rate and suspended sediment yield, and local climate data (precipitation and temperature). It was concluded that LUCL changes caused some important modifications in the riparian corridor, riverbed size, and river flow rate. A 40–80% reduction in the river bed area was observed, vegetation developed in the riparian buffer strips, and the river channel changed from braided to a single channel. The main causes identified are reductions in the river flow rate and suspended sediment yield (−36% and −38%, respectively), while climate change did not have a significant effect.

**Keywords:** farmland abandonment; LULC changes; climate change; runoff/suspended sediment changes; river morphology dynamics; Italian Apennines

#### **1. Introduction**

The development and specific characteristics of rivers and streams are influenced by surrounding landscapes [1,2]. Our current understanding of rivers' dynamics incorporates a conceptual framework of spatial nested controlling factors in which climate, geology, and topography at large scales influence geomorphic processes that shape channels at intermediate scales [3]. However, direct human impact on the environment cannot be neglected at a local scale, especially in the last century. In particular, land use/land cover (LULC) changes have a significant impact on basin water cycles and soil erosion dynamics. Human factors also include water abstraction for irrigation, flow regulation, the construction of reservoirs, and mining. An extensive bibliography analyzing the effects of dams and reservoirs on the geomorphic responses of rivers has been produced [4–7]. However, the influence of other drivers, such as climate and LULC changes has been much less documented over time, although recent studies highlight their importance in inducing river changes [8–11].

It can be said that agriculture and land abandonment are two complementary aspects of human impacts on the landscape. Land abandonment can affect net soil losses [12], while recolonization of natural vegetation can lead to a reduction in soil loss and a progressive improvement of soil characteristics [13]. Moreover, land abandonment and agriculture can also lead to changes of river stream morphologies, in particular narrowing and incision [13]. Liebault and Piegay [14] observed on the Roubion River (France), that colonization of unstable gravel bars tends to channel minor floods which, in turn, form a new and narrower channel in the existing river bed. A number of studies [3,15] have documented statistical links between LULC and stream conditions, using multisite comparisons and empirical models.

Hydrological alteration is one of the principal environmental factors by which LUCL influences stream ecosystems [3]. It alters the runoff-evapotranspiration balance, can cause an increase or decrease in flow rate's magnitude and frequency, and often lowers river's base flow. In addition, hydrological alteration contributes to a change of channel dynamics, including increased erosion of the channel and its surroundings, and a less frequent overbank flooding [3]. Zhang and Schilling [16] noted that increasing streamflow in the Mississippi River was mainly due to an increase in base flow, which in turn was a consequence of LULC (conversion of perennial vegetation to seasonal row crops). Many researchers have studied the effects of LULC changes on river flow and most of them have indicated that intensified afforestation will reduce both runoff peak and total runoff volume [17]. On the other hand, even a modest riparian deforestation in highly forested catchments can result in the degradation of a stream habitat, owing to sedimentary input. A comparison of different catchments showed that an increased forest area results in lower concentrations of suspended sediments, inferior turbidity at base flow, lower bed-load transport, and less embeddedness [3].

Numerous studies have demonstrated that LULC, the abandonment of rural activities, and consequently, a decrease of human pressure on mountain areas, has contributed to increase of vegetation cover [18–20]. In the case of Reno River mountain basin, Pavanelli et al. [21] documented that recolonization of natural vegetation and a consequent increase of actual evapotranspiration, was the key hydrological variable that caused the decrease of the river flow rate. However, not many studies have addressed the effect of the redevelopment of natural vegetation at the river-basin scale on river flow, sediment yield, and riverbed morphology. Picco et al. [22] noted a consistent increase of riparian vegetation within the corridor of the Piave River (Northern Italy) during the last five decades, concluding that it depended on human activities, both in the main channel and at basin scale. LULC and local climate change (e.g., precipitation, temperature and evapotranspiration), may induce notable alterations of watershed hydrology [13,23]. Several studies showed that precipitation increase alone is insufficient to explain increasing flow rate trends in agricultural watersheds [24,25], since changes in agricultural land use can also result in increased flow rate.

Collectively, these studies provide strong evidence for the importance of the surrounding landscape, and human activities for the hydrological and morphological characteristics of rivers [3,15]. However, it is often difficult to separate human from naturally driven activities [26]. In addition, most of these articles were mainly conducted on small spatial scales, within a few hundred meters of a stream, without considering larger spatial units. Finally, only few studies have addressed the effect of redevelopment of natural vegetation at the river basin scale on river discharge, sediment yield, and bed river morphology.

Therefore, the aim of this paper was to explore relationships between local climate change that occurred in the last century and agricultural land abandonment (and consequent LULC changes), which culminated in the 1950s, on the one hand; and modifications in morphology and hydrology of the Reno River (Northern Italy), on the other hand.

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

#### *2.1. Study Area*

The Reno River, located in the Northern Italy, flows into the Adriatic Sea. It is the sixth biggest Italian river, with a catchment area of 5965 km2 and a length of 211.8 km. The mountain hydrographic network of the Reno is rather ramified and dense, and it is composed of eight major rivers, 12 secondary rivers and 600 torrents. This study concentrates on one part of the Reno River; namely, the Reno River

Mountain Basin (RMB), located in the Northern Italian Apennines (Emilia Romagna and Tuscany Regions), with a catchment area of 1061 km<sup>2</sup> and length of 80 km. The RMB's average altitude is of 639 m a.s.l.; it ranges from a maximum elevation of 1945 to a minimum one of 60.35 m a.s.l. at the dam of Chiusa of Casalecchio (44◦47' N, 11◦28' E), which is the RMB outlet (Figure 1).

**Figure 1.** Orography map of the Reno River Mountain Basin with the main river channel in blue (**left**) and the simplified land use/land cover (LULC) of the area in 2003 (**right**).

The RMB can be considered representative of the environmental and anthropogenic changes that have occurred in the Italian central and northern Apennines in the last century. The Apennine agricultural area of the Emilia Romagna Region (RER) decreased by almost 50% from 1960 to 2000 [27]. After 1950, the population moved to the cities and valley, a phenomenon that affected the whole Italian Apennines. Until World War II, this area had an average population density of 85 inhabitants km<sup>−</sup>2, with agro-forestry and pastoral farming as the main activities. Currently the density is reduced to less than 70% of the previous one [27]. After World War II, due to industrialization and the development of agricultural mechanization, the landscape rapidly transformed. Agriculture remained where it was cost-effective, while the rest gave way to permanent meadows, scrub, and woodland [28]. Currently, land cover is characterized by oak woods, beeches, shrubs, and pastures at higher altitudes. Chestnut woods are present at medium altitudes and on coppices, pastures, and crops on hillsides. Crops, vineyards, orchards and urban areas cover the catchment valley (Figure 1).

The RMB consists mainly of erodible sedimentary rocks. Land cover, runoff, soil erosion, and suspended sediment in the river are closely related to each other. The bedrock consists of resistant limestone, sandstone, and meta-sandstone in the upper part of the watershed and of weakly cemented marl, mudstone, sandstone, and conglomerate in the middle and lower part of the watershed [29]. The fluvial terraces are preserved from the outlet (Casalecchio) to about 20 km upstream (~150 m a.s.l.). Upstream of that point, landslides and earthflows preclude significant preservation of terraces.

The RMB average precipitation is 1305 mm year−<sup>1</sup> (Table 1), with the following distribution: Winter 337 mm, spring 307 mm, summer 182 mm, and autumn 411 mm. The average temperature is 10.7 ◦C; July is the hottest month, with peaks up to 29.5 ◦C (1998). Winters are generally very cold, with the average minimum monthly temperature dropping to −8.9 ◦C (January 1942). The fluvial regime of the Reno River is linked to rainfall, with floods occurring in autumn and spring. Seasonal floods are characterized by short times of concentration (time it takes to reach the basin outlet), owing to the long and narrow shape of the catchment. The evaluations of the 30 and 200-year recurrence interval floods at Casalecchio (Chiusa) are 1541 and 2280 m<sup>3</sup> s<sup>−</sup>1, respectively.


**Table 1.** Available data set and its properties.

\* The two periods correspond to the two meteorological stations—one in the valley and another in the mountains, that were analyzed separately in order to highlight differences between them.

The dam "Chiusa of Casalecchio" is the RMB outlet. It is the oldest hydraulic building in Europe, and has allowed the social and economic development of the city of Bologna through hydraulic energy since the XII century. It is also included in the UNESCO program's list of Patrimony Messengers of a Culture of Peace. The Chiusa dam controls the downstream base level, making the reach from the source to the Chiusa geomorphologically independent. Census of the hydraulic works detected 51 weirs on that reach, and most of them were built in the first decades of the 20th century. Similarly, in Rio Maggiore, a tributary of Reno, there is a total of 162 dams, which equals a 10 dam km−<sup>2</sup> density [28]. Moreover, in the early 1900s, five hydroelectric dams were built in the tributaries of the Reno. Even though dams and reservoirs do not influence river water budget on the longer time scale, they do act as river sediment traps and can affect the river flood regulation [21]. During the last few decades, part of water has been diverted for domestic and irrigation purposes; however, these withdrawals represent a very limited percentage (3%) of the Reno river's flow rate [21].

#### *2.2. Methodology*

The parameters considered in this study were LULC (since the 1950s), river morphology (in 1954 and 2003), and hydro-climate changes (since the 1920s) of the RMB. They were used to assess different environmental changes that occurred in the RMB and whether if they were due to agricultural land abandonment and a consequent renaturalization after the 1950s. The main hypothesis of this research is that the year 1960 is the date from which effects of these changes can actually be seen.

#### 2.2.1. Land Use Changes

As already said, the year 1960 was taken as a point when effects of big LULC changes started to be visible, and the year 1954 was taken as representative of the period before agricultural land abandonment. The two LULC maps are both to the scale 1:25000 (Soil Use Maps RER 1954 and 2003) were derived from the black and white aerial photographs flight G.A.I. 1954 (Military Geographic Institute-Italy) and the satellite images from 2003 (Quickbird). They were analyzed and compared with Geographic Information System software ARCVIEW 3.2 (Environmental Systems Research Institute, Redlands, CA, USA) to evaluate the RMB LULC changes. The relatively heterogeneous classes of the two different maps were reclassified according to the CORINE land cover classes [30] in order to be compared. Then, 20 and 51 land cover classes of the 1954 and 2003 Soil Use Maps respectively, were grouped into urban areas, water bodies, fallows, forests, and crops. This approach reduced the error due to the different sources of images of the maps. In Figure 1 fallows and forests are grouped together to highlight natural vegetation and renaturalization effects.

#### 2.2.2. Morphological River Changes

To assess changes to the river belt and river morphology, aerial photos from the year 1954 (GAI-IGMI) and satellite images from Quickbird, 2003 (Figure 2), were used. The effects of LULC changes and human impact on the Reno river terraces and banks were analyzed on a reach from Porretta Terme (349 m a.s.l.) to Casalecchio di Reno (60 m a.s.l.), with a length of 54.5 km and a mean bed river slope of 0.53%. The area was within 250 meters from each side of the river bed, done using GIS software. The transects (Figure 2 left) were drawn using the same reference points for both 1954 and 2003, to estimate variations in the river bank width, vegetation, and stream bed. The river's morphological changes, the width, the river's channel area, and the riparian LULC of each transect, were evaluated. The river channel was identified in the images as a non-vegetated part of river bed corresponding to the physical confine of the normal water flow. Banks, on the other hand, are subject to water flow only during high water stages and a riparian buffer strip is the vegetated area near a stream. All of them constitute the river corridor.

**Figure 2.** The Reno river transects (in yellow) between Casalecchio di Reno and Porretta Terme on an IGMI-GAI photo from 1954 (**left**), and the reaches R1, R2, and R3 on a satellite image from Quickbird, 2003 (**right**).

The studied 54.5 km portion of the Reno river was split into three reaches (Figure 2, right), on the basis of similar morphological characteristics (mean width and slope of the river and its banks): R1 from Porretta Terme to Vergato (21 km long), R2 from Vergato to Sasso Marconi (23 km long) and R3 from Sasso Marconi to Casalecchio Chiusa (10.5 km long). The properties of the three reaches were obtained from field survey of the transects and from the 1954 and 2003 images. The cartographic and photo-interpretation data were managed with GIS ARCVIEW 3.2 software used for map drawing and calculating the extent of the areas covered by vegetation. On the other hand, data on the river morphology (e.g., type of riverbed material) and vegetation type in the transects were obtained during field surveys.

#### 2.2.3. Climate and Hydrological Data

The Reno River hydrological and RMB climate data (Table 1), were processed on a monthly, yearly and seasonal basis, and they were divided into two periods—before and after 1960, the year that was taken as a date from which effects of LULC changes and renaturalization can be seen. The monthly data of the river flow rate (Q) and the suspended sediment yields (SSY) are came from samples collected at the outlet of the RMB (Casalecchio di Reno gauge, 60 m a.s.l.) by the Italian Hydrographical Service (SIMI) and by the Regional Agency for Environmental Protection of the Emilia Romagna Region (ARPAE).

In addition, precipitation and temperature data collected by the SIMI and ARPAE were analyzed to evaluate impact of climate change on the RMB. For the minimum and maximum monthly temperatures of the two RMB stations: One in the mountains (Monteombraro, at 704 m. a.s.l.) and one in the valley (Anzola at 42 m. a.s.l.) were examined (Table 1).

#### 2.2.4. Data Analysis

Statistical analyses were performed on the data records with the STATGRAPHICS®Centurion XVI software (StatPoint Technologies, Inc., The Plains, VA, USA). Statistical tests were used to verify whether the river flow rate data (Qmean and Qmax) detected before 1960 (1921–1959) were statistically different from the data collected after 1960 (1960–2013). Two samples were compared using F test, discriminant analysis, box and whisker plots, *t*-tests, the Kolmogorov–Smirnov Test, and analysis of variance tests. Discriminant analysis was run on flow rate samples to verify if each value was correctly classified in the two periods. The *t*-tests were run to compare means of the two samples and to test the null hypothesis that the two means were equal. The results were verified with a Kolmogorov–Smirnov test and the distributions of the two samples were compared, for both the Qmean and the Qmax. Finally, a box and whisker plot was run to demonstrate a significant difference (p < 0.05) between the flow rate data coming from the two periods considered.

Moreover, the river flow rate, SSY and climate data were also analyzed with linear trend analysis in order to show linear regression in the examined periods with 95% confidence limits and the prediction limits of the least squares fit model. Trend-line slopes (b) were calculated by the least-square linear fitting method.

#### **3. Results and Discussion**

#### *3.1. LULC Changes*

At present, more than 60% of the entire mountain Reno catchment is covered with forest [31], while in the past, forest was scarce due to exploitation (Figure 3). Cultivated land in the RMB catchment decreased from about 37% in 1954 to 5% in 2003, partially losing space to forest, which is, in fact, mainly 50–60 years old (Figure 4). In detail, forests and fallows increased from 39.5% to 57% and from 19% to 28%, respectively. Urban areas increased from 0.45% (1954) to 6.5% (2003), with the majority of population concentrated in the Reno valley. Lastly, water surfaces were found to be 40% smaller—they reduced from 9.3 km<sup>2</sup> in 1954 to 5.7 km2 in 2003 (Figure 4a). The disappearance of rocky outcrops that were predominantly clayey badlands, the so called "Calanchi," is also interesting: They decreased from 1.67% to 0.46% of the RMB area, a phenomenon that was repeated throughout the Apennines.

**Figure 3.** Photos of Tresca Mountain (1473 m. a.s.l.) near Porretta Terme (RMB): In the past, say 1914 (**left**), forests were sparse due to exploitation, while currently the area is forested (**right**).

**Figure 4.** LULC changes in the RMB between 2003 and 1954 (**a**), and the LULC of the riparian buffer strips (**b**).

LULC changes noted in the riparian buffer strips between 1954 and 2003 were an increase of forests, urban areas, and uncultivated land, at the expense of cultivated land (that decreased from 60.25% to 11%) (Figure 4b). That decrease was the most prominent near the city of Bologna (reach R3, Figure 2), where cultivated land decreased from 85% (1954) to 20% (2003), while urban areas and woods showed an increase, from 7% to 42%, and from 8% to 38%, respectively.

The most widespread species (40–50%) in the riparian buffer strips are *Populus nigra* and *Salix alba*. Besides them, other species present are *Quercus pubescens* (about 20%) and *Salix alba* (10–15%), that form mixed populations, and different arboreal and shrub forms. *Alnus glutinosa* (about 10%) is present in the upper river banks where the anthropic impact is less. The shrub layer is composed of *Sambucus nigra*, *Corylus avellana*, *Cornus sanguinea*, and *Prunus spinosa*. Exotic species that have found, in this fluvial habitat, the ideal conditions of development, are *Robinia pseudoacacia* and *Acer negundo*, while among the shrubby, an infesting species is *Amorpha fruticose*.

#### *3.2. Hydrology Changes*

The monthly mean flow rate of the Reno River is 23.4 m3 s<sup>−</sup>1, while the maximum monthly value recorded was 143 m<sup>3</sup> s−<sup>1</sup> in December 1959 (Table 2). Between 1921 and 2013, the mean yearly flow rate was reduced by 11 m<sup>3</sup> s−<sup>1</sup> (b = <sup>−</sup>0.12) or 36% (Figure 5), while the Qmax reduction was about 30.3%. The correlation coefficient (R2) values of linear regression for Qmean and Qmax were 0.22 and 0.08, respectively. Since the *p*-values were both less than 0.05, there was a statistically significant relationship at 95% confidence level. The monthly flow rate values are given in the box and whisker plot (Figure 6a). Starting from the hypothesis that the renaturalization and consequent hydrological changes caused by agricultural land abandonment in the 1950s were detectable after 1960, river flow rate data was divided into two sub-periods: Before and after this date. The Qmean values were 26.03 m3 s−<sup>1</sup> and 21.08 m<sup>3</sup> s<sup>−</sup>1, for the 1921–1959 and 1960–2013 periods, respectively (Figure 6b). The F-test and *p*-value were, respectively, equal to 9.79 and 0.0026, according to the ANOVA test. Since the *p*-value was less than 0.05, there was a statistically significant difference between the two flow rate means.

**Table 2.** Summary statistics of mean and maximum annual flow rate (Q) and suspended sediment yield (SSY)—average yearly data.


**Figure 5.** Time series of the Reno River mean yearly flow rate and linear trends with 95% confidence limits (green) and the prediction limits (grey) of the least squares fit model.

**Figure 6.** Box and whisker plots of monthly flow rates between 1921 and 2013 (**a**), and of the subpopulations: 1921–1959 and 1960–2013 (**b**).

Discriminant analysis was run on the two groups of Qmean: Before and after 1960. From the 77 values used to fit the model, 72.7% were correctly classified in the groups, out of which 66.7% and 78% were for 1921–1959 and 1960–2013, respectively. A statistically significant difference (p < 0.05) was found between the two groups, for both Qmean and Qmax (Figure 6b). It was evident that the lowest flow rate values were mostly concentrated during the 1960—2013 sub-period. In fact, there was a marked change in the trend of Qmax and Qmean around 1960. Additionally, significant differences between the two periods were shown by discriminant analysis, as dispersion and as data trend. Dispersion of the mean flow rate data for the first period (1921–1959) was higher in respect to the second one (1960–2013), as evidenced by correlation coefficients R2 (0.04 and 0.20 respectively), and therefore the two groups are statistically different. This difference, particularly the higher value of dispersion data of the first period, indicated presence of floods events, including disastrous ones [30].

#### *3.3. Suspended Sediment Yield (SSY)*

LULC change is an important factor that affects soil erosion-runoff and SSY. The relationships between climate and LULC changes on one side and SSY on the other, was investigated by PJ Ward et al. [32], with use of geo-referenced model WATEM/sedem. The authors found that the sediment increase in the Meuse, a northern European river, was almost entirely due to LULC change (conversion of forests to agricultural land). They concluded that increase of riparian buffer strip and development of riparian vegetation can result in the reduction in SSY, as they can be used as barrier to soil runoff [33].

Similarly, in the Apennines, SSY can be used to assess a real loss of the catchment soil due to runoff, rills erosion, gullies, and badlands, but also losses due to agricultural land abandonment and LULC change [33,34]. Currently, soil erosion prevails in the badlands and agricultural areas that are concentrated on the slopes, and that are easily accessible to mechanization. The average yearly SSY in the period 1942–1978 was 934.3 Mg km<sup>−</sup>2. The maximum monthly value was 960 Mg km−<sup>2</sup> in December 1976 and the yearly maximum value was 2225 Mg km−<sup>2</sup> in 1951 (Table 2). Yearly and seasonal SSYs are given in the Figure 7a. It can be seen that the analysis of the SSY linear trend indicated a 17.5% reduction during the year 31 of data, or a 38% reduction in the period 1921–2013.

**Figure 7.** Seasonal SSY and the linear trend of yearly SSY (**a**). Seasonal SSY variations of 1942-59 and 1960-78 years respect to seasonal average of 1942 to 1978 years (**b**).

For the period 1942-1978, the seasonal SSY averages were: 426 Mg km−<sup>2</sup> (winter), 231 Mg km−<sup>2</sup> (spring), 32 Mg km−<sup>2</sup> (summer) and 244 Mg km−<sup>2</sup> (autumn) (Figure 7a). The average SSYs for the two sub-periods 1942–1959 and 1960–1978 (Figure 7b) were 981 and 902 Mg km−<sup>2</sup> respectively. Moreover, some interesting seasonal variations occurred. For example, compared to the 1942–1978 seasonal average, SSY in the second period was lower in winter and autumn (−4.5% and −9.6%, respectively). On the other hand, it increased in spring and summer by 2.6% and 13.2% respectively (Figure 7b).

In general, the highest average monthly values of SSY were in February and November (213 and 188 Mg km−2, respectively), but they decreased by 37% and 32% in the two following decades. If seasons of the two periods are compared, it can be noted that SSY reduced in summer and spring, while it increased in winter and autumn (Figure 7b). Finally, linear relations between the yearly SSY and Qmean and Qmax, showed that SSY is better correlated to the maximum flow rate (Figure 8).

**Figure 8.** Linear relation between SSY and mean (Qmean) and maximum (Qmax) flow rate.

The average annual erosion in the RMB, calculated on the basis of 31 year of SSY data, was about 9.3 Mg ha−<sup>1</sup> year−<sup>1</sup> or 0.6 mm year−1. Various bibliography estimates available for the Region of Emilia Romagna give different results. The European Agency for the Environment, using the model Pesera [35], estimated a soil loss of 2.42 Mg ha−<sup>1</sup> year−1, slightly below the average Italian value (3.11 Mg ha−<sup>1</sup> year<sup>−</sup>1). In addition, the Emilia Romagna Region (including the RMB) was defined as a soil erosion risk area [36]. It was found that around 21% of the region has a medium to high risk of soil loss. The average loss in higher grounds of the region was estimated to be around 6 Mg ha−<sup>1</sup> year−<sup>1</sup> [37]. The difference between the erosion value reported in this study and other researchers' ones is mainly due to the heterogeneity of estimation models and basic data being estimated. Considering the current conditions of vegetation cover, they are all lower than those calculated on the basis of historical data from the 1942–1978 period.

#### *3.4. Climate Change*

An important aspect of a basin water budget is climate: Temperature and precipitation trends. In fact, temperature variations influence hydrology and river flow rate, and since they have an impact on evapotranspiration from vegetation, water surfaces and soil. Figure 9a shows a linear increase for both minimum and maximum mean yearly temperature in the RMB. That is visible for the mountain and valley gauge. The minimum temperature (Tmin) showed a similar rising trend for the two stations: +4 ◦C/100 years (R2 = 0.49) for the Mountain gauge and +5 ◦C/100 years (R2 = 0.59) for the Valley gauge. On the other hand, Tmax trends were more complex—decreasing in the mountain areas (−1.9 ◦C/100 years; R<sup>2</sup> = 0.15) and showing an increase in the valley (+2.8 ◦C/100 years; R<sup>2</sup> = 0.21). Since only the minimum, hence night, temperatures increased in the mountain areas where vegetation is also more developed than in the valley, and since night-time evapotranspiration is much smaller than in the day-time, the temperature increase most probably did not have a big influence on the observed flow rate reduction.

Precipitation in the RMB was lowered by 10.67% between 1921 and 2013 (Figure 9b), corresponding to 145 mm reduction in 92 years. However, that value is not statistically significant (R<sup>2</sup> = 0.034). A strong reduction (from about 0.6 to 0.4) in the catchment runoff coefficient (flow rate/precipitation), that was observed in the last 90 years (R<sup>2</sup> = 0.43) (Figure 10), and is mainly due to reduction of the flow rate.

**Figure 9.** Local climate change: Linear trends of maximum and minimum temperature (Tmax and Tmin) for mountain (Mg) and valley (Vg) gauges (**a**); yearly RMB precipitation with linear trend with 95% confidence limits (green) and the prediction limits (grey) of the least squares fit model (**b**).

**Figure 10.** Runoff coefficient observed in the last 90 years (R<sup>2</sup> = 0.43) with linear trends, prediction and confidence limits.

#### *3.5. Morphological Stream Changes: Riparian Bu*ff*er Strips*

The main proprieties and changes of the three reaches considered, based on aerial and satellite images, and field surveys, are reported in Table 3. Figure 11 gives the relationship (exponential equation, R2 = 0.75) between the average width and mean altitude of the Reno valley. The two values that are above the curve are of two villages (Porretta and Pioppe di Salvaro) that are in a larger area due to tributary torrents. The changes of the Reno River reaches (R1–R3) concern the banks and the river morphology (Figure 12), from upstream to downstream:



**Table 3.** Main features of the Reno reaches.

**Figure 11.** RenoRiver transect between Casalecchio Chiusa and Porretta Terme, with correspondinglocations.

**Figure 12.** Reach R1 near Porretta that is more torrential (**left**) and R3 near Chiusa Casalecchio with fine sediments (**right**): Riparian forest shows a strong development.

**Figure 13.** Reach 2 near Marzabotto village. Wide stream and area covered with crops in 1954 (**left**), but urban areas and a narrow stream in 2003 (**right**).

Based on the RMB's soil use maps, it was estimated that the riverbed area decreased by about 40% from 1954 to 2003. However, based on the river transects, the reduction was as high as 80%. Fluvial banks with woods along the Reno were mostly absent in 1954, and the area was used for farming. In 2003 riparian forests appear to be well developed along the entire stream in the RMB. This is especially visible on the right-hand side (looking downstream) of the river, where fluvial terraces are narrower or absent, and are therefore under a lower human impact. In general, it is observed that the river bed gives way to riparian forests. For example, a reduction in the active riverbed corresponds to formation and/or expansion of the riparian buffer strips, and the stream reaches are colonized by riparian vegetation.

#### **4. Conclusions**

This study has examined relationships between two major geomorphological changes (channel narrowing and formation of wide vegetated banks) that took place in the Reno River mountain basin in the last century, and the hydrological, climatic, and basin re-naturalization factors that have contributed to these changes. The two phenomena are strongly and positively covariant, indicative of cause and effect, and in fact, wide vegetated banks are formed at the expense of the riverbed. While riparian buffer strips were mostly absent in 1954, currently they are well developed along the entire stream. In addition, the shape of the river channel changed from braided to a single one and the width of the river bed was reduced by around 80%. The riverbed in the past, mostly consisted of gravel bars and sand deposits, and currently clay and silt prevail at the basin outlet.

Based on this study, the steering factors for those changes and significant aspects were:


The effect of agricultural land abandonment that occurred in the 1950s can be recognized after 1960, confirming the initial hypothesis that this year can be taken as a starting point for the basin change. After that date a decrease in the Reno flow rate was observed and dispersion of data is significantly reduced. All statistical analyses confirm that the hydrological flow data measured after 1960 (period 1960 to 2013) are significantly different from those measured when the basin was still heavily agricultural (period 1921 to 1958). However, although this study has identified the human factor as one of the main causes of the above-mentioned changes, it can often be challenging to separate human from naturally driven activities, and future research is needed in order to do it. The geomorphological evolution of the Reno River shows how these changes are mainly related to the hydrological dynamics and catchment re-naturalization. Although climate changed in the period studied (precipitation reduction of 10.67% and 4–5 ◦C increase of Tmin) it had little bearing on the observed environmental changes.

This study, applied to a typical North Apennine river, illustrates the effectiveness of combining historical data (hydrological and climate data, so as aerial and satellite images) on the one hand, and the use of modern technology (geographic information systems) and direct surveys on the other. Combining these techniques can certainly contribute to a sustainable management of river systems. Although further research is needed, this study gives an insight to the past and present factors that regulate water course, its hydrology, and morphology. An in-depth knowledge of this factors can certainly make it possible to predict the evolution and dynamics of the Reno river flow and its morphology.

**Author Contributions:** Conceptualization, D.P.; methodology, D.P. and C.C.; validation, D.P.; formal analysis, D.P.; data curation, D.P. and C.C.; writing—original draft preparation, D.P.; writing—reviewing and editing, D.P., S.L., and A.T.; visualization, D.P. and S.L.; supervision, D.P. and A.T.

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

**Acknowledgments:** This research was undertaken as part of the collaboration agreement for teaching and study with the Emilia Romagna Region, Civil Protection, Basin Authority of the Reno River and ARPAE-RER, (Italy).

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

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Hydrologic Impacts of Land Use Changes in the Sabor River Basin: A Historical View and Future Perspectives**

**Regina Maria Bessa Santos 1, Luís Filipe Sanches Fernandes 1, Rui Manuel Vitor Cortes <sup>1</sup> and Fernando António Leal Pacheco 2,\***


Received: 3 July 2019; Accepted: 10 July 2019; Published: 15 July 2019

**Abstract:** The study area used for this study was the Sabor river basin (located in the Northeast of Portugal), which is composed mostly for agroforestry. The objectives were to analyze the spatiotemporal dynamics of hydrological services that occurred due to land use changes between 1990 and 2008 and to consider two scenarios for the year 2045. The scenarios were, firstly, afforestation projection, proposed by the Regional Plan for Forest Management, and secondly, wildfires that will affect 32% of the basin area. In this work, SWAT (Soil and Water Assessment Tool) was used to simulate the provision of hydrological services, namely water quantity, being calibrated for daily discharge. The calibration and validation showed a good agreement for discharge with coefficients of determination of 0.63 and 0.8 respectively. The land use changes and the afforestation scenario showed decreases in water yield, surface flow, and groundwater flow and increases in evapotranspiration and lateral flow. The wildfire scenario, contrary to the afforestation scenario, showed an increase in surface flow and a decrease in lateral flow. The Land Use and Land Cover (LULC) changes in 2000 and 2006 showed average decreases in the water yield of 91 and 52 mm·year<sup>−</sup>1, respectively. The decrease in water yield was greater in the afforestation scenario than in the wildfires scenario mainly in winter months. In the afforestation scenario, the large decrease varied between 28 hm3·year−<sup>1</sup> in October and 62 hm3·year−<sup>1</sup> in January, while in the wildfires scenario, the decrease was somewhat smaller, varying between 15 hm3·year−<sup>1</sup> in October and 49 hm3·year−<sup>1</sup> in January.

**Keywords:** SWAT; water balance components; Land Use and Land Cover changes; wildfires; afforestation

#### **1. Introduction**

Land Use and Land Cover (LULC) are considered the most critical factors affecting the intensity and frequency of surface flow, as well as soil erosion and the loss of nutrients [1,2]. Watershed-level studies have indicated that rapid LULC changes could have significant impacts on water resources [3] and were reported to be an essential factor for controlling water resources on local and global scales during the last century [4,5]. The factors responsible for land use changes are population density fluctuations, changes in the agricultural policy, and the conditions of national and international markets. The main changes occurred due to the abandonment of farmland in less productive mountain areas, the expansion of some subsidized crops to marginal lands, and intense soil erosion during extreme rainstorm events.

Several studies developed in the Mediterranean region showed that inadequate land use and soil cover accelerate water erosion processes and, consequently, lead to land degradation [1,2,6]. Soil erosion is one of the leading causes of the reduction of water quality due to the amount of sediment that arrives at the watercourses and reservoirs [1]. However, many authors have demonstrated that both runoff and sediment loss decrease exponentially as the percentage of vegetation cover increases [1,6]. Thus, forests are often used as a management strategy to improve the provision of ecosystem services in watersheds, because they have a substantial widespread positive influence on climate, hydrology, soils, and biodiversity [7–10].

The studies developed in the 19th century, based on the too-hasty generalization of single point observation, believed that natural and planted forests increased total flow and base flow [3]. Nowadays, studies prove that forests have an impact on the water balance at the basin scale, as forest water consumption is generally higher than that of other vegetation types [3,10,11]. The rise in shrub and forest cover produces declines in water yield, surface, and groundwater flows, and an increase in transpiration. The decrease in the groundwater flow in forests can be justified by the fact that trees with deep roots and high transpiration rates may act as pumps that remove water from the soil and return it to the atmosphere [1,11,12].

In the Mediterranean, erosion is the consequence of complex interactions between environmental and human-related factors. The erosion processes are products of the occurrence of intense rainstorms and long-lasting droughts, high evapotranspiration, the presence of steep slopes, topographic diversity, and recent tectonic activity as well as the recurrent use of fire, deforestation, overgrazing, farming, and construction activities [1,2]. In agreement with the CORINE program, Spain and Portugal are the Mediterranean countries in the European Union facing the highest risk of erosion [13,14]. CORINE means "coordination of information on the environment", and one of the priorities of the program is the evaluation of natural resources and environmental problems in the southern part of the European Community (e.g., soil erosion, water resources, land cover, and coastal problems) [14]. In Portugal, areas of high erosion risk cover almost one-third of the country [15]. For this reason, soil erosion is one of the most intensively studied issues in the Mediterranean region, Portugal included [1,2,6,15]. Erosion and land degradation became a problem in Portugal when arable farming expanded into marginal areas, namely cereal cultivation until the middle of the twentieth century [1]. The introduction of modern agriculture led to the abandonment of traditional or semi-traditional agriculture in mountainous areas as well as in areas with difficult access, resulting in fundamental transformations to the landscape, characterized by the spread of natural vegetation, including both shrubland and forestland [1]. The natural afforestation has resulted in a decline in water resources and surface flow and has decreased soil loss and sediment delivery, as well as caused a progressive improvement in soil characteristics [1]. It has been extensively documented that the rise in shrub and forest cover has promoted lower soil losses and sediment yields than in arable land [2,6].

However, wildfires have been responsible for sudden increases in erosion rates, and they can also result in land degradation and sometimes desertification over the long-term [16]. The occurrence of forest fires affecting thousands of hectares each year is a significant problem in the Mediterranean basin [2,16]. In Portugal, fires are essentially human-caused, and their extension and severity are dependent on extreme weather [17]. In central Portugal, at Pedrógão Grande-Góis, a tragic wildfire occurred on June 17, 2017, with an official death toll of 64 people, almost 500 buildings destroyed, and a continuous patch of more than 42 thousand hectares burned in one week [18]. The areas burned correlated well with socioeconomic and environmental characteristics (e.g., population density and land use, climate, weather, topography, and vegetation cover) [16,17]. However, in this tragic fire, the climate and meteorology played key roles in the initiation and spreading of the wildfires [18]. The warmer and drier than average spring made the landscape prone to the occurrence of large fires and extreme weather conditions favored the ignition and spread of wildfires [18].

It is generally accepted that fires increase runoff, soil erosion, and nutrient and pesticide transport to the river [16,19]. In the Mediterranean region, the most substantial erosion rates and nutrient losses tend to occur during the first rainstorms after wildfire occurrence during the dry season. The increases in runoff are due to the reduced infiltration capacity of the soil caused by the removal

of vegetation and soil organic matter by fire [16]. Wildfires lead to the release of nutrients through the combustion of vegetation and exposure of the soil to erosive factors. Consequently, wildfires are responsible for hydrologic problems and the degradation of water quality caused by excessive input of nutrients, such as phosphorus and nitrogen [16,20]. The excess nutrients affect primary production and, consequently, the eutrophication of aquatic systems. The relationship between the risk of eutrophication and nutrient exports from burned areas has been widely documented [16,19,20].

Hydrological models are computational tools that perform a mathematical representation of hydrological processes, such as infiltration of water into the soil, recharge of aquifers, runoff and drainage network flow [21,22], as well as hydrochemical processes such as weathering or contaminant transport [23–31]. They are also the basis of decision support systems, which can help watershed managers in the control of extreme events such as floods [32,33], or in the assessment of water resource availability [34–40] and threats to water quality [41–48]. The SWAT (Soil and Water Assessment Tool) was the hydrological model used to model the physical processes that occurred in the Sabor river basin. The SWAT model is one of the most widely used water quality watershed and river basin scale models worldwide [49]. It has been applied from small hydrographic basins to the continental scale. Some examples of applications are in North America [50], Europe [51], and Australia [52,53]. These studies were done in four main categories: hydrologic modelling, sediment transport, nutrient and pesticide transport, and scenario analyses.

Within this framework, the following three demands were addressed: (i) establishment of the hydrologic baseline with calibration over a 39-year period (1960–1999), (ii) assessment of the land cover and land use changes between 1960 and 2008 and their effects on water balance components, and (iii) assessment of the forecast of the afforestation and wildfire scenarios for 2045, also on water balance components. A majority of articles discussing environmental effects of land use change and dealing with scenario creation only assess some water balance components, namely water yield or surface flow. This work intends to study the changes occurring not only in water yield or surface flow but also in evapotranspiration, lateral flow, and in groundwater flow.

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

#### *2.1. Study Area*

The Sabor river basin is located in the northeast of Portugal and drains into the International Douro basin (Figure 1). The Sabor river has an extension of 212.6 km. Its source is located in Spain at an altitude of about 1600 m, and its mouth is located in the Douro River at about 88 m [54]. The Sabor river basin covers an area of about 3834.5 km2, of which 3170.7 km2 is in Portugal. The average slope of the basin is 16.2% according to the digital elevation model [54]. The slope below 10% is located in the upper zones and plateaus of the basin, and the most rugged slope is located on the scarps and banks of the Sabor river and its tributaries. The main tributaries are the Maças river (drainage area: 720 km2) and the Angueira river (drainage area: 540 km2). In the Sabor river, the Sabor hydroelectric system was built. It comprises two retention dams, known as upstream and downstream walls, which are located at 12.6 and 3 km from the mouth of the Sabor river, respectively [55].

The climate is close to the Mediterranean type. It is characterized by warm-dry summers and precipitation concentrated in the winter and spring seasons. The average annual basin values of rainfall and evapotranspiration are approximately 730 and 540 mm·year<sup>−</sup>1, respectively [56].

The soil characterization in the Sabor river basin is mainly lithosols (87% of the catchment area), but also with the presence of cambisols (7.4%), alisols (3%), anthrosols (1.4%), and fluvisols (0.7%) [57]. The 1990 CORINE Land Cover data published by the European Environmental Agency [58] identified the area as having 59% agricultural areas, 32% semi-natural areas, 9% forest, and less than 1% artificial areas and water bodies (Figure 1b). Agriculture and livestock farms are the main economic activities in the region, with olive and almond being the main crops [59]. According to the 2011 demographic census, in the Sabor river basin, the population density was 20 inhabitants per km2 [60].

**Figure 1.** The spatial distribution of (**a**) the drainage network, the topography, and the sub-basins of the Sabor river basin; (**b**) the temperature and the CORINE Land Cover in 1990; and (**c**) the area burned between 1990 and 2017.

#### *2.2. CORINE Land Cover Changes and Wildfires*

In the study of LULC, changes were used to form the maps of the CORINE Land Cover for 1990, 2000, and 2006 [58]. These maps were reclassified to SWAT land cover classes, which is the format read from ArcSWAT. The reclassification is found in the Supplementary Material (worksheet 1). Figures 2 and 3 illustrate the changes between the SWAT land cover classes of 1990 and 2000 and between the SWAT land cover classes of 1990 and 2006. These changes represented 6% of the basin area for the SWAT land cover classes of 1990 and 2000 (Figure 2a), and 10% of the basin area between 1990 and 2006 (Figure 2b). For better visualization and analysis, in Figure 3, only areas greater than 500 hectares are shown. The major changes in the area of SWAT land cover classes from 1990 to 2000 comprise the replacement of range—brush by forest (7080 ha), essentially coniferous (FRSE, 4009 ha) and deciduous forests (FRSD, 2212 ha), as well as burned areas (BARR, 1643 ha). The major change in the area of SWAT land cover classes from 1990 to 2006 was the increase in range—brush (RNGB, 14,384 ha) mainly from range—grasses (RNGE, 3635 ha), agricultural land—generic (AGRL, 2584 ha), forest—mixed (FRST, 2395 ha), forest—deciduous (FRSD, 1915 ha), forest—evergreen (FRSE, 1725 ha), and burned area (BARR, 1421 ha). As well as the replacement of agricultural land—row crops (AGRR) in agricultural land—generic (AGRL, 4155 ha).

The cartography of the burned areas covering the period 1990–2017 is available at the Conservation of Nature and Forests [61] (Table 1). In this period in the Sabor river basin, 103,033 ha was burned, which corresponds to 32% of the basin area (Figure 1c). The largest burned area (approximately 13,000 ha) occurred in 2013, followed by burned areas above 6700 ha, which occurred in 1994, 1998, 2000, 2012, and 2017 (Figure 4).

**Figure 2.** The spatial distribution of changes (**a**) between the CORINE Land Cover 1990 and 2000; and (**b**) between the CORINE Land Cover 1990 and 2006, in the Sabor river basin. Only the areas with changes greater than 500 hectares are presented. The SWAT codes are agricultural land—row crops (AGRR), agricultural land—generic (AGRL), range—brush (RNGB), range—grasses (RNGE), barren (BARR), forest—deciduous (FRSD), forest—evergreen (FRSE), and forest—mixed (FRST).

**Figure 3.** Graphic of changes between the CORINE Land Cover in 1990 and 2000, and between the CORINE Land Cover in 1990 and 2006 in the Sabor river Basin. Only the areas with changes greater than 500 hectares are presented. The SWAT codes are agricultural land—row crops (AGRR), agricultural land—generic (AGRL), range—brush (RNGB), range—grasses (RNGE), barren (BARR), forest—deciduous (FRSD), forest—evergreen (FRSE), and forest—mixed (FRST).

**Figure 4.** The area burned between 1990 and 2017 in the Sabor river basin.

**Table 1.** Maps and records were used in the construction of the hydrological model, in the evaluation of Land Use and Land Cover (LULC) changes between 1990 and 2008, as well as in the analysis of the afforestation and wildfire scenarios. All data types used have reference to the purpose, their owner institution, and the source of data.


#### *2.3. Conceptual Framework Model*

The conceptual framework model was developed to assess the effects of LULC changes and the forecast of afforestation and the occurrence of wildfires on water balance (Figure 5). SWAT was used to construct the hydrological model of the Sabor river basin with the following set of data: digital elevation model, land use, soil types, and weather station (Table 1). With these data, the drainage network, the delimitation of the basin and sub-basins, and the hydrologic response units (HRU) were defined. After that, the model was calibrated by SWAT-CUP with streamflow discharge data on a daily basis. Then, the hydrological model was simulated with the values of the parameters obtained in the calibration procedure. This model was used in the diagnostic phase and in future scenarios. In the diagnostic phase, the model was used to assess the effects of LULC changes between 1960 and 2008 on water balance components with the CORINE Land Cover in 1990, 2000, and 2006. In future scenarios, the model was used to evaluate the effects of the forecast of afforestation and the occurrence of wildfires on water balance components.

**Figure 5.** The conceptual framework model for constructing and modelling the water resources in the SWAT (Soil and Water Assessment Tool).

#### *2.4. SWAT Model*

The SWAT is based on a physical process to be simulated in a watershed [49]. In a watershed, SWAT continuously simulates the hydrological pattern processes, water balance, water quality, nutrients, and sediment exportation [49,50]. For simulations, information collected from various sources (e.g., topography, land cover) is necessary, and this is then assimilated into the database (e.g., soil characteristics, meteorological data). For modelling purposes, the SWAT model divided the watershed into a sub-basin linked in a cascade by the stream network [62]. In turn, the land area in a sub-basin is divided into HRUs. Each HRU is a portion of a sub-basin that comprises unique land cover, soil, and slope attributes [49]. The SWAT treats the HRU as a homogeneous unit of land use, management techniques, and soil properties and then quantifies the relative impacts of vegetation, management, soil, and climate change within each HRU [63]. The output of the hydrological model (e.g., runoff, sediments, nutrients) is calculated separately from each HRU and then added together to determine the total loading from a sub-basin. The advantage of HRUs is the increase in accuracy which adds to the prediction of loading from a sub-basin [49].

The SWAT uses a modified Soil Conservation Service Curve Number (SCS CN) methodology, and the Penman–Monteith equation to calculate the runoff and evapotranspiration, respectively [50]. Curve numbers have been developed and published for a wide range of land cover types and uses and can be found in [64]. The Penman–Monteith equation requires daily values of precipitation

(mm), maximum and minimum temperature (◦C), solar radiation (MJ/m2 day), relative humidity (%), and wind speed (ms<sup>−</sup>1), but this observed data is usually available with gaps and thus may limit the performance and results of the model [49,62]. To fill the gaps, SWAT includes the WXGEN stochastic weather generator model to generate climatic data [63,65]. The WXGEN generates daily weather information that is missing from the monthly average data summarized over a number of years [49].

#### *2.5. SWAT Input Data*

The 2012 version of ArcSWAT was used to build a hydrological model of the Sabor river basin. ArcSWAT is an ArcGIS extension and graphic user input interface for SWAT, which is used worldwide and is continuously under development [66]. We selected ArcSWAT because it has been used in numerous hydrologic, decision-making, and environmental applications [62,65,67], and the authors have experience with using ArcGIS for the processing, overlay, and combination of multiscale and multi-type spatial data in thematic surveys or projects focused on the collection and interpretation of spatial data [34,68–73].

The input data used to construct the hydrological model were (i) the topography of Trás-os-Montes and Alto Douro to generate a digital elevation model and slope, (ii) the drainage network of the Sabor river basin to define the stream and delineate the basins and sub-basins, (iii) the CORINE Land Cover to create a land use grid, (iv) the soil types of Trás-os-Montes and Alto Douro to create a soil grid, and (v) weather data of stations located inside or near the basin to simulate climatic data (Figure 6, Table 1).

A precipitation dataset was compiled between 1957 and 2008 by the various climatologic stations located inside and near the basin (Figure 6b). The gaps in the daily precipitation time series were calculated by inverse distance weighted interpolation. This method is frequently used in climatic predictions and has already been proven to provide good results [17,62]. The time series of weather information insert in the WXGEN weather generator model was provided by the SNIRH meteorological station of Folgares (06N/01C) (Figure 6b). The WXGEN filled the daily missing data based on an average of 46 years of weather data.

A total of 37 sub-basins were defined in the basin area (Figure 1a), with an average area and standard deviation of 86 and 69 km2, respectively. A total of about 523 HRUs were defined within sub-basins, with a threshold value of 10% for the land use, soil, and slope classes. The SWAT model was executed on a daily basis from 1957 to 1999 with a warm-up period of 3 years. A warm-up period is recommended to initialize the simulation process with the objective of ensuring the establishment of basic flow conditions and hydrologic processes equilibrium as well as to help to minimize the model values for the initial hydrological conditions [62,67].

The land use grid used in the construction of the hydrological model was the CORINE Land Cover 1990 (CLC) (Table 1), but this land cover does not provide information that corresponds to the SWAT land cover classes. Therefore, the CLC classes were reclassified, which led to several CLC classes having the same SWAT land cover classes (Figure 6c) [62,65]. For example, all classes of the heterogeneous agricultural areas of CLC classes were generalized into the SWAT land cover class agricultural land—generic. Among the different SWAT land cover classes, the most appropriate reclassification for the burned areas of CLC classes was barren (BARR). The reclassification is presented in the Supplementary Material (worksheet 1).

The State Soil and Geographic (STATSGO) is the soil database integrated into the SWAT [74]. These soil categories are unavailable in Portugal, and a match could not be found between those categories and the available ones for the research area. The soil map available in the Sabor river basin was extracted from the digital soil map of the Trás-os-Montes and Alto Douro region (Table 1). The soil categories were lithosols, cambisols, alisols, anthrosols, fluvisols, and urban land (Figure 6d). Therefore, data, including the soil component parameters and soil layer parameters, were inserted into the SWAT database for each soil type, except for urban land, which already existed.

**Figure 6.** The SWAT (Soil and Water Assessment Tool) input data used for constructing the hydrological model were (**a**) the elevation; (**b**) the slope, meteorological, and hydrometrical data; (**c**) the SWAT codes of land cover classes; and (**d**) the soil types. The SWAT codes are agricultural land—close-grown (AGRC), agricultural land—generic (AGRL), agricultural land—row crops (AGRR), barren (BARR), forest—deciduous (FRSD), forest—evergreen (FRSE), forest—mixed (FRST), vineyard (GRAP), olives (OLIV), orchard (ORCD), pasture (PAST), range—brush (RNGB), range—grasses (RNGE), commercial (UCOM), residential—high density (URHD), residential—low density (URLD), transportation (UTRN), and water (WATR).

The hydrological model was calibrated from 1960 to 1999 and validated from 2000 to 2008. A hydrometric station called Quinta das Laranjeiras (Figure 6b) was used to calibrate and validate the streamflow on a daily basis.

#### *2.6. SWAT-CUP Model Calibration*

The computer program used for calibrating the SWAT 2012 models was the SWAT-CUP 2012 (Calibration and Uncertainty Procedures) [75]. The SWAT-CUP 2012 consists of five different calibration procedures, including functionalities for validation and sensitivity analysis. The calibration procedure SUFI-2 (Sequential Uncertainty Fitting) was used in this work. The SUFI-2 algorithm is quite efficient for large-scale, time-consuming models [75,76]. In SUFI-2, the uncertainty analysis is based on the discrepancy assessment between observed and simulated values, taking into account potential sources of uncertainty, like observed data, the conceptual model, and parameters. The degree of uncertainty is quantified by the 95% prediction uncertainty (95PPU), measured by the P-factor and R-factor. The P-factor is the percentage of measured data bracketed by the 95PPU and varying from 0 to 1, wherein 1 indicates a 100% bracketing fit of the observed values. The R-factor measures the calibration quality and indicates the thickness of the 95PPU. A P-factor of 1 and an R-factor of 0 mean a perfect fit for the observed and calibrated values [75].

SUFI-2 has several criteria and quantitative statistical methods to evaluate the outcome simulation values from SWAT, as compared with the observed data [75]. The objective function selected as the calibrated parameter set was the coefficient of determination (R2). The statistical methods used to assess the model performance were the following: the Nash–Sutcliffe efficiency (NS), the percent bias (PBIAS), the ratio of the root mean square error to the standard deviation of measured data (RSR), and the coefficient of determination (R2). The model performance is considered satisfactory whenever <sup>R</sup><sup>2</sup> and NS are greater than 0.5, RSR is less than 0.7, and PBIAS is less than <sup>±</sup>25% for streamflow [77,78].

#### *2.7. Diagnostic Phase of Land Use and Land Cover Changes*

In the diagnostic phase, CLC 1990, CLC 2000, and CLC 2006 were used to determine the effects of LULC changes on the water balance components of the Sabor river basin. The SWAT hydrological model was constructed and calibrated with CLC 1990. In forthcoming sections, this hydrological model will be referred to as the reference model. After that, two more hydrological models were constructed, one with the CLC 2000 and another one with the CLC 2006. The values of the calibration parameters of the reference model were inserted into these two hydrological models. Using the same parameter values in all hydrological models ensures that changes in the streamflow are exclusively due to LULC changes. Thus, the study of the LULC changes consisted of comparing the values of the water balance components between the reference model and both the CLC 2000 hydrological model and the CLC 2006 hydrological model. The reference model was simulated between 2000 and 2008, the CLC 2000 hydrological model was simulated between 2000 and 2005, and the CLC 2006 hydrological model was simulated between 2006 and 2008.

#### *2.8. Future Scenarios of A*ff*orestation and Wildfires*

In the Sabor river basin, the afforestation scenario was based on the Regional Plan for Forest Management (RPFM) of Douro [79] and Northeast [80], and the wildfires scenario was based on the burned area between 1990 and 2007 [61]. The Regional Plan for Forest Management (Portuguese Regulate Decree no. 3/2007, published on 17 January 2007) includes a plan to be accomplished until 2045, and the main objectives are to reduce the risk of wildfire occurrence, and to adjust the proportions of resinous and deciduous species based on the application of correct forest management models. In the RPFM, the Sabor river basin occupies eleven homogeneous sub-regions in the Douro and Northeast regions, which are referenced in the Supplementary Material (worksheet 2).

In order to create a model of the afforestation scenario, a map was created with the coniferous and broad-leaved areas proposed by RPFM until 2045. Afforestation consisted of counting the areas of coniferous, broad-leaved, and mixed forest in CLC 1990 and calculating the missing areas until reaching the percentage of afforestation proposed by the RPFM for each sub-region (Figure 7a). According to this technical report, the area of afforestation of agricultural land will include 40% coniferous forest and 60% broad-leaved forest, and the afforestation of semi-natural areas will include 70% coniferous forest and 30% broad-leaved forest [79,80]. To make the map, the afforestation of agricultural land was carried out in the following order of CLC classes: agro-forestry areas, land mainly occupied by agriculture, with significant areas of natural vegetation and complex cultivation patterns. The afforestation of semi-natural areas was done in the following order of CLC classes: transitional woodland-shrub, sclerophyllous vegetation, moors and heathland, and burned areas. A map of the afforested area is presented in Figure 7b, and calculation of the respective areas is found in the Supplementary Material (worksheet 2). The map of the afforestation scenario inserted in the SWAT was the CLC 1990 updated with afforestation areas. Also, the map of the wildfires scenario inserted in the SWAT was the CLC 1990 updated with all burned areas in the Sabor river basin between 1990 and 2017 (Figure 7d).

**Figure 7.** The maps of afforestation and wildfires scenarios inserted in the SWAT for constructing the hydrological models were (**a**) broad-leaved, coniferous, and mixed forest from the CORINE Land Cover 1990 and the projected percentages of afforestation per sub-region for 2045, proposed by the Regional Plan for Forest Management; (**b**) the afforestation projection for 2045; (**c**) areas burned between 1990 and 2017; and (**d**) the CORINE Land Cover 1990 updated with areas burned between 1990 and 2017.

After creating the maps, two hydrological models were constructed, one with the map of the afforestation scenario and another with the map of the wildfires scenario. The values of the calibration parameters of the reference model were inserted into these two hydrological models. Once again, the use of the same parameters in both hydrological models allowed the scenarios to be compared with the reference data. In other words, the factors used to compare the afforestation and wildfires scenarios with the reference model included the water yield surface flow, evapotranspiration, lateral flow, and groundwater flow.

#### **3. Results**

#### *3.1. Calibration and Validation of the Streamflow*

The streamflow was calibrated for a 39-year period (1960–1999) and validated for a 9-year period (2000–2008) both on a daily basis. The parameters and their respective values resulting from the calibration are shown in Table 2. The graphics in Figure 8a,b show good agreement between the observed and simulated streamflow for both calibration and validation. The calibration and validation, illustrated in Figure 8a,b, respectively, were performed on a daily basis, but for best visualization, they are represented on a monthly basis. The goodness-of-fit indicators for the streamflow calibration (Table 3), based on the R2, RSR, and NS show satisfactory performances (with values of 0.63 and 0.62) and PBIAS shows a very good performance (2.7%) [50,51]. The same goodness-of-fit indicators were obtained for the validation with a very good performance for R2 (with 0.8), and satisfactory performances for RSR, NS, and PBIAS with 0.63%, 0.61%, and −24%, respectively (Table 3). The negative value of PBIAS indicates a model overestimation bias [77].


**Table 2.** The parameters used in the calibration procedure of streamflow. In the legend of methods, R is relative and V is the replacement value.


**Table 3.** Goodness-of-fit indicators for daily calibration between 1960 and 1999 and validation of streamflow between 2000 and 2008 in the Sabor river basin.

**Figure 8.** The comparison of observed and simulated streamflow during (**a**) the calibration (between 1960 and 1999); and (**b**) validation (between 2000 and 2008) in the Sabor river basin. The simulation of the streamflow was executed on a daily basis, but for best visualization was present on a monthly basis.

#### *3.2. Land Cover and Land Use Changes*

Figure 9 shows the graphics of the water balance components of the reference model (with the CLC 1990, simulated between 2000 and 2008) the hydrological model of the CLC 2000 (simulated between 2000 and 2005), and the hydrological model of the CLC 2006 (simulated between 2006 and 2008). The results show that the LULC changes that occurred in 2000 and 2006 led to a decrease in the water yield and an increase in evapotranspiration (Figure 9a,b). The water yield decreased by an average of 91 and 52 mm·year−<sup>1</sup> for the LULC changes in 2000 and 2006, respectively. The evapotranspiration increased by an average of 90 and 55 mm·year−<sup>1</sup> for the LULC changes in 2000 and 2006, respectively. The values of surface flow and groundwater flow decreased, while the lateral flow increased (Figure 9c,e). On average, for the LULC changes in 2000 and 2006, the surface flow decreased by 28 and 23 mm·year<sup>−</sup>1, and the groundwater flow decreased by 91 and 50 mm·year<sup>−</sup>1, respectively. The lateral flow increased by 10 and 5 mm·year−<sup>1</sup> for the LULC changes in 2000 and 2006, respectively.

**Figure 9.** Quantity of water (mm) (**a**) in the water yield; (**b**) evapotranspiration; (**c**) surface flow; (**d**) lateral flow; and (**e**) groundwater flow. The blue bars are the reference model (simulated between 2000 and 2008 with the CLC 1990), the orange bars are the CLC 2000 model (simulated between 2000 and 2005), and the grey bars are the CLC 2006 model (simulated between 2006 and 2008).

#### *3.3. A*ff*orestation and Wildfires Scenarios*

The values of the water balance components were calculated on a monthly basis between the reference model (CLC90) and both the afforestation and wildfires scenarios. The results are illustrated in Figures 10 and 11, which represent the monthly values and spatial distribution, respectively. The graphics of Figure 10a,b show that, in both scenarios and in every month of the year, the water yield decreased and the evapotranspiration increased. The decrease in the water yield was more pronounced in the rainy season (autumn and winter) than in the dry season (spring and summer) and in the afforestation scenario compared with the wildfires scenario. In the afforestation scenario, these large decreases varied from 28 hm3·year−<sup>1</sup> in October to 62 hm3·year−<sup>1</sup> in January, while in the wildfires scenario they were somewhat smaller, varying from 15 hm3·year−<sup>1</sup> in October to 49 hm3·year−<sup>1</sup> in January. The month of August registered much lower water yield decreases: 3 hm3·year−<sup>1</sup> in the <sup>a</sup>fforestation scenario and 2 hm3·year−<sup>1</sup> in the wildfires scenario. The evapotranspiration increase was more pronounced in January, February, and March, with values ranging from 46 to 82 hm3·year−<sup>1</sup> for the afforestation scenario and from 33 to 61 hm3·year−<sup>1</sup> for the wildfires scenario. In the afforestation scenario, the lowest increase in evapotranspiration was registered in August (3 hm3·year<sup>−</sup>1), while in the wildfires scenario, a small decrease in evapotranspiration was registered in September (0.4 hm3·year<sup>−</sup>1).

**Figure 10.** Quantity of water (hm3·year<sup>−</sup>1) (**a**) in the water yield; (**b**) in evapotranspiration; (**c**) in surface flow; (**d**) in lateral flow; and (**e**) in groundwater flow. The blue bars show the difference between the reference model and the afforestation scenario model, and the orange bars show the difference between the reference model and the wildfires scenario model.

The water available for surface flow and lateral flow showed different behaviors in both scenarios. In the afforestation scenario, the surface flow decreased and the lateral flow increased, and the opposite occurred for the wildfires scenario (Figure 10c,d). In the afforestation scenario, the major decrease in the surface flow occurred in January, November, and December, with values ranging from 31 to 37 hm3·year<sup>−</sup>1, while the smallest decrease occurred in July with 0.8 hm3·year<sup>−</sup>1. The major increase in the lateral flow occurred between October and January with values ranging from 10 hm3·year−<sup>1</sup> in October to 13 hm3·year−<sup>1</sup> in December, while the smallest increase occurred in July (0.2 hm3·year<sup>−</sup>1). In the wildfires scenario, the surface flow increased between January and October with values ranging between 2 hm3·year−<sup>1</sup> in July and 8 hm3·year−<sup>1</sup> in February, and it decreased in November and December with values of 1 and 6 hm3·year<sup>−</sup>1, respectively. The lateral flow decreased in every month of the year, with values ranging from 1 hm3·year−<sup>1</sup> in August to 6 hm3·year−<sup>1</sup> in January. The major decrease in the groundwater flow occurred in the rainy season in both scenarios (Figure 10e). In the afforestation scenario, these large decreases varied from 31 hm3·year−<sup>1</sup> in October to 69 hm3·year−<sup>1</sup> in January, while in the wildfires scenario, they were somewhat smaller, varying from 25 hm3·year−<sup>1</sup> in October to 62 hm3·year−<sup>1</sup> in January. The lowest groundwater flow was registered in August, with 4 hm3·year−<sup>1</sup> in the afforestation scenario and 3 hm3·year−<sup>1</sup> in the wildfires scenario.

**Figure 11.** SWAT output maps of (**a**–**c**) the surface flow and (**d**–**f**) the groundwater flow. The surface flow matches (**a**–**c**) the reference model, the afforestation scenario, and the wildfires scenario, respectively. The groundwater flow matches (**d**–**f**) the reference model, the afforestation scenario, and the wildfires scenario, respectively.

SWAT output maps of the surface flow and the groundwater flow for the reference model and scenarios are illustrated in Figure 11. Figure 11a–c shows the surface flow of the reference model, for the afforestation scenario, and for the wildfires scenarios, respectively. Figure 11d–f shows the groundwater flow of the reference model, for the afforestation scenario, and for the wildfires scenario, respectively. The maps show that the highest surface flow and groundwater flow values in the

reference model (Figure 11a,d) and scenarios (Figure 11b,c,e,f) are located in the sub-basins of the Northwest. When comparing the maps of the reference model with the scenarios, it can be seen that the surface flow and groundwater flow values decrease. However, these decreases are major in the afforestation scenario.

#### **4. Discussion**

#### *4.1. Performance of Streamflow and TP Calibration*

The performance statistics of the model for the daily streamflow during calibration and validation periods were satisfactory and good, respectively. Similar results were achieved by other authors [81,82] using the SWAT for streamflow simulation in both the calibration and validation periods in the Alto Sabor river basin (located in the upstream part of the Sabor river basin). The results of the performance statistics indicated that SWAT was able to capture and reproduce the average flows and seasonal variations of the Sabor river basin (Figure 8a). However, the SWAT also shown had difficulty reproducing the observed extreme values of the streamflow, as shown by the regular underestimation of the most important peak events for the calibration period. Many studies have also detected the inability of the SWAT models to reproduce the highly complex hydrological processes during the high flow periods [52,53,81,82]. Major high flows or extreme conditions occurred in the winter of the years 1960, 1966, 1969, 1985, and between 1976 and 1979 (Figure 8a). Several authors analyzed how the rainfall-runoff translation occurs in reality compared with that simulated by the model [52]. These authors verified that the same amount of annual rainfall did not correspond to the same amount of observed streamflow. In contrast, the model presented a similar streamflow for the same amount of rainfall. The authors stated that the curve number method in SWAT uses the average daily rainfall and does not take into account the intensity and duration of rainfall. This is a factor that can explain why peak flows were underestimated by the model. One example of the high intensity of rainfall was registered in four basins of the North of Portugal, between 2000 and 2006, which had a considerable number of days that exceeded 20 mm of precipitation [68].

Despite the goodness-of-fit indicators for the streamflow calibration being considered satisfactory, there are some potential factors which compromise the performance of the model. We identified the following potential factors: (i) the climate stations cannot represent the orographic effect on the precipitation rates, (ii) the meteorological data are considerably affected by orographic variations, (iii) the human activities across the catchment may not have been accounted for in the model (e.g., artesian wells and boreholes, irrigation, and domestic consumption), (iv) the observed flow values are not 100% accurate. The precipitation gauges are located between 441 and 905 m of elevation with an average altitude of 659 m, and they cannot represent the orographic effect which varies between 88 and 1464 m of elevation. The orographic variations have a considerable effect on meteorological data; for example, the vertical temperature gradient, the depressions, and the shade in the valley result in lower minimum temperatures. These local variations are not entirely covered by weather stations. The last factor that influences the effectiveness of calibration is the problem that the hydrometric gauges faced during data collection. The discharge flow was estimated from hydrometric levels which were transformed into streamflow by a state discharge curve, but it was only valid within a range of hydrometric heights. Thus, the calculation used was not able to capture the actual hydrological pattern within extreme events.

#### *4.2. E*ff*ects of LULC Changes and Future Scenarios for Water Resources*

A major change in the LULC was the replacement of range—brush (scrub and herbaceous vegetation associations) with forest (essentially coniferous and deciduous) from CLC 1990 to CLC 2000 and the replacement of range—grasses (natural grasslands), agricultural land—generic (heterogeneous agricultural areas), forest (coniferous, deciduous, and mixed) and burned area in range—brush from CLC 1990 to CLC 2006 (Figure 3). These LULC changes are a result of a complex process of plant recolonization of shrubs and forests as a consequence of farmland abandonment. This natural reforestation has the advantage of limiting soil erosion. However, the authors of [56,83] showed that many Mediterranean areas evolved into simplified and continuous landscapes with plenty of flammable wood, which led to fire occurrences. The wildfires increased soil erosion and progressively reduced the potential for recolonization of plants, leading to stony soils. This explains the many degraded panoramic landscapes in Mediterranean areas [16].

The results showed that the LULC changed from CLC 1990 to 2000, and from CLC 1990 to CLC 2006, causing decreases in the water yield, surface flow, and groundwater flow and increases in evapotranspiration and lateral flow. The major changes were observed in groundwater flow and evapotranspiration. The decrease in the groundwater flow in 200 and 2006 was, on average, 91 and 50 mm·year−1, respectively (Figure 9a). The increase in evapotranspiration in 2000 and 2006 was, on average, 90 and 55 mm·year−1, respectively (Figure 9b). The same results were obtained with the afforestation projection scenario proposed by the Regional Plan for Forest Management to 2045. In both the afforestation scenario and LULC changes, the decreases in the water yield, surface flow, and groundwater flow were due to a large amount of water being lost by evapotranspiration.

The increase in the lateral flow suggests that the forest provides an increase of water infiltration into the soil, but the high transpiration rate of the trees removes the water from the groundwater into the atmosphere. The increase in soil water infiltration by afforestation was confirmed by Ilstedt et al. [12]. The authors applied a meta-analysis based on several articles and concluded that the infiltration capacity increased, on average, by approximately three-fold after afforestation in agricultural fields. Another study [11] selected 20 studies to compare water flows in tropical watersheds under natural or planted forests and non-forest lands to provide useful results for valuing the watershed ecosystem services. The planted forests were Pines and Eucalyptus and lowland natural forests. The main results showed significantly lower total flows or water yields and groundwater flows under planted forests than under non-forest land uses. These results were explained by the high transpiration rates of Eucalyptus. According to these results, the authors argue that the effect of natural forests on base flow results from two competing processes: the high infiltration under the forest contributes to soil water recharge and the base flow increase, while the high transpiration of the trees contributes to the base flow decrease. Furthermore, the infiltration of soil water may be higher under planted forests than under non-forest land uses, but may not be sufficient to offset the loss of water by transpiration [12]. The compilation of paired-watershed results, based on a total of 137 basins, showed that reforestation results in a decrease in water yield and deforestation results in an increase [3]. Several authors have observed that reforestation could reduce the amount of sediment entering streams [84], as vegetation reduces the streamflow and increases infiltration [8,12]. The results of Serrano-Muela et al. [85] confirmed that forest conservation in the Central Spanish Pyrenees reduces floods and soil erosion, particularly on steep slopes. In this study, the major decreases in water yield, surface flow and groundwater flow, as well as the increases in lateral flow and evapotranspiration, were obtained in the wet season. For example, the decrease in water yield varied between 28 hm3·year−<sup>1</sup> in October and 62 hm3·year−<sup>1</sup> in January. Moreover, an increase in the lateral flow occurred between October and January with values ranging from 10 hm3·year−<sup>1</sup> in October to 13 hm3·year−<sup>1</sup> in December. The evapotranspiration was also higher in winter, because there was more water available to evaporate, but the highest values were obtained in February and March (with values ranging between 46 and 82 hm3·year<sup>−</sup>1). The reason for this is the combination of the rainfall and the increase of temperature that support evaporation and the growth of the canopy trees, which increases transpiration.

In contrast, in the wildfires scenario, increases in the surface flow and evapotranspiration and decreases in the lateral flow and groundwater flow were observed (Figure 10b–e). Similarly to the afforestation scenario, the major increases or decreases of water in the wildfires scenario mainly occurred in winter. For example, the water yield varied between 15 hm3·year−<sup>1</sup> in October and 49 hm3·year−<sup>1</sup> in January, and the increase in evapotranspiration was more pronounced in January, February, and March, with values ranging between 33 and 61 hm3·year<sup>−</sup>1.

Wildfires are agents of change that can dramatically alter evaporation and transpiration rates, and their short-term effects on surface evaporation rates may be complex. The increase in runoff after the fire has been attributed to destroyed vegetation as well as to reduced evapotranspiration [16,86–88]. The study by Miranda et al. [89] showed that immediately after a fire, evapotranspiration rates decrease. However, this situation is reversed with the first rains, because the burned area rapidly becomes a stronger sink for CO2 and has higher evapotranspiration rates than nearby unburned areas. This difference persists throughout the wet season and is attributable to the greater physiological activity of the growing vegetation in the burned area. Similar rates of evaporation immediately after rainfall for both burned and unburned plots were observed in Campo Sujo, near Brazil [90]. However, within a few days of the soil drying, lower evaporation rates from the burned plot were observed. About 1 month after the fire, evaporation rates from the burned plot were typically as much as 80% that of the unburned control early in the wet season [90]. This was attributed to greatly enhanced rates of soil evaporation in the burned plot, especially immediately after the rain. Nagra et al. [91] also obtained similar results and justified the increase in the evaporation rates post-fire as a result of low albedo and reduced vegetation cover (lack of shading).

The increase in surface flow and the decreases in lateral flow and groundwater flow were due to the reduced infiltration capacity of the soil caused by the removal of vegetation by fire. Previous studies [19,92,93] showed that changes in vegetation cover and topsoil (composed of organic matter) have important impacts on the hydrological regime. A review of the literature developed in the European Mediterranean showed that wildfires make the soil more susceptible to removal by water erosion, less likely to allow infiltration, and more likely to promote surface flow [16]. The increase in surface flow, as well as soil erosion caused by fires, were proven in several studies in Portugal [19,92,93] and in the Mediterranean region [16,94,95]. For example, in burnt Eucalyptus and Pinus forest plantations in the Águeda Basin, North-Central Portugal, the surface flow on burned plots was 5%–25% higher than on unburnt ones [96]. For burned pine plantations in a small catchment located in Central Portugal, a runoff of up to 48.5% was found (1.1 km2) [19]. In the Arbúcies basin, North-East Spain, there was an increase of 30% in flood runoff [94]. In the Rimbaud catchment in South-East France, an increase in the annual runoff of 30% was also verified during the first year after the fire [86].

A reduction in the infiltration capacity and the absence of vegetation improved the flood peak and soil erosion [16,93]. Nunes et al. [1] reported that during the study period (2005 and 2006), more than 80% of the total rainfall fell in autumn and winter. Consequently, the monthly rainfall erosivity, based on the Modified Fournier Index, increased considerably during these months. This result suggests that the energy available for erosion and transport was highly concentrated in this season. They also verified that the significantly larger amounts of runoff caused high soil loss. The overland flow was three times higher in 2006 than in 2005 which resulted in a sediment yield twice as high. However, the soil losses fundamentally depended on the amount and intensity of rainfall. The stability of the soil structure and aggregates is usually thought to be reduced by fire, producing more easily eroded soil [16].

Even so, modest post-fire soil losses could be important for soil longevity in some areas, because of soil organic matter and nutrient losses in solution or adsorbed onto eroded sediment particles [16,97]. Much of the post-fire nutrient content is in the form of ash, which is prone to removal by wind and water [19]. However, it can also benefit the soil quality through nutrient-rich ash becoming incorporated into the soil [16]. The high fire frequency and thin soils in a nutrient-poor ecosystem, typical of many fire-prone Mediterranean areas causes the risk of soil fertility depletion to be high [16,19].

#### *4.3. Overview of the E*ff*ects of the LULC Changes on Water Resources*

Changes in the LULC of the Sabor river basin have had impacts on the water balance. The main changes were the increases in forest and shrubs between 2000 and 2008. These changes caused decreases in the water yield and groundwater flow and an increase in evapotranspiration as well as the afforestation scenario. The decrease in water yield was due to the reductions of surface flow and groundwater flow. The decrease in the surface flow was due to the increase in soil water infiltration by forests. Additionally, the decrease in groundwater flow can be partly attributed to the high transpiration rate of the trees. Despite forests causing substantial reductions in water yield, they also improve the water quality, reduce erosion, and benefit biodiversity. According to Francis et al. [10], a low to intermediate tree cover can improve soil hydraulic properties by up to 25 m from its canopy edge, which means that the hydrologic gains can be proportionally higher than the additional losses from the increased transpiration.

In the afforestation and wildfires scenarios, a decrease in the water yield was observed as well as an increase in evapotranspiration. The increase in evapotranspiration in the wildfires scenario was due to the evaporation of the water as a result of low albedo and a reduction in the vegetation cover. The comparison of the two scenarios showed that the trees were responsible for the reduction of groundwater flow. However, they contributed to the decrease in sediment yield and nutrients in surface flow. On the other hand, the wildfires scenario increased the surface flow as well as the sediment yield and the concentration of nutrients in the surface flow. The recurrence of fires has long-term implications because of the increase in soil erosion and the loss of nutrients, making it difficult for plants to recolonize. The water quality was affected by the increase in nutrient concentration due to eutrophication of the aquatic environment when associated with the increase in temperatures in the summer. The construction of two dams in the Sabor river for the production of electricity has increased the problem of eutrophication. According to the technical report [98] and the article [99], the changes in water quality caused by dam construction and the consequent stream water impoundment were significant. The increases in temperature and electric conductivity and accumulation of phosphorus and nitrogen in the reservoirs triggered the growth of algae, the increase of chlorophyll a, and a drop in the transparency of the water. The consequences of water deterioration for the aquatic fauna were severe, marked by abrupt declines of native fish species and the invasion of exotic species.

#### **5. Conclusions**

In this study, the SWAT2012 model was applied to assess the spatiotemporal dynamics of the water balance in the Sabor river basin. First of all, the LULC changes between 1960 and 2008 were analyzed with the CLC 1990, CLC 2000, and CLC 2006. Secondly, two scenarios were created, afforestation and wildfires. The afforestation scenario used was the projection proposed by the Regional Plan for Forest Management until 2045. The wildfires scenario used was based on areas burned in the Sabor river basin between 1990 and 2017.

The overall results of LULC changes between 1960 and 2008 show that the increase in forest and shrubs is a consequence of farmland abandonment in a basin where the population density is less than 20 inhabitants per km2. The changes caused decreases in the water yield, surface flow, and groundwater flow and increases in evapotranspiration and lateral flow. The major changes occurred in the water yield and evapotranspiration. The water yield decreased in 2000 and 2006, on average, by 91 and 52 mm·year−1, respectively. Evapotranspiration increased in 2000 and 2006, on average, by 90 and 55 mm·year−1, respectively. The decrease in groundwater flow and increase in evapotranspiration were attributed to the high transpiration rate of the trees. Similar results were obtained for the afforestation scenario: decreases in the water yield and groundwater flow and an increase in evapotranspiration. The major decrease occurred in winter for groundwater flow with values varying between 28 hm3·year−<sup>1</sup> in October and 62 hm3·year−<sup>1</sup> in January. The increase in evapotranspiration was more pronounced in January, February, and March, with values ranging from 46 to 82 hm3·year<sup>−</sup>1.

In contrast, in the wildfires scenario, increases in the surface flow and evapotranspiration and decreases in the lateral flow and groundwater flow were observed. Similarly to the afforestation scenario, the major increases or decreases of water in the wildfires scenario occurred in winter. The major decrease in groundwater flow occurred in winter, for which values varied between 25 hm3·year−<sup>1</sup> in October and 62 hm3·year−<sup>1</sup> in January. The increase in evapotranspiration was more pronounced in January, February, and March, with values ranging from 33 to 61 hm3·year<sup>−</sup>1.

The afforestation projection, proposed by the Regional Plan for Forest Management, will cause greater decreases in water yield, surface flow and groundwater than the wildfires scenario. A decrease in groundwater has implications for streamflow, because it is this flow that feeds the rivers in the dry season. However, it generates an increase in lateral flow, with more available water for irrigation. This is an important factor, because the basin is essentially agroforestry. With this study, we have shown that afforestation cannot simultaneously maximize all environmental benefits, but any form of afforestation should provide environmental improvements over agricultural land. Reductions in water yields could be minimized by planting trees with low water use at low densities and by avoiding landscape positions with access to the groundwater. An improvement in water quality could be achieved by using the contrasting strategies.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/7/1464/s1, The Supplementary Materials comprise the following worksheets: Worksheet 1—Reclassification of CORINE Land Cover classes into SWAT counterparts; Worksheet 2—Forest occupation and predicted afforestation until 2045, in the eleven homogeneous sub-regions in the Douro and Northeast regions.

**Author Contributions:** Conceptualization, R.M.B.S. and F.A.L.P.; methodology, R.M.B.S. and F.A.L.P.; software, R.M.B.S.; validation, F.A.L.P.; formal analysis, L.F.S.F.; investigation, R.M.B.S.; resources, L.F.S.F. and R.M.V.C.; data curation, R.M.B.S., F.A.L.P., L.F.S.F., and R.M.V.C.; writing—original draft preparation, R.M.B.S.; writing—review and editing, F.A.L.P.; visualization, R.M.B.S.; supervision, F.A.L.P. and L.F.S.F.; project administration, R.M.V.C.; funding acquisition, R.M.V.C.

**Funding:** This research was funded by the INTERACT project "Integrated Research in Environment, Agro-Chain and Technology", no. NORTE-01-0145-FEDER-000017, in the line of research entitled BEST "Bio-economy and Sustainability", and co-financed by the European Regional Development Fund (ERDF) through NORTE 2020 (the North Regional Operational Program 2014/2020). For authors at the CITAB research centre, this research was further financed by the FEDER/COMPETE/POCI—Operational Competitiveness and Internationalization Programme under Project POCI-01-0145-FEDER-006958, and by National Funds of FCT—Portuguese Foundation for Science and Technology under the project UID/AGR/04033/2019. For the author in the CQVR, this research was additionally supported by National Funds of FCT—Portuguese Foundation for Science and Technology under the project UID/QUI/00616/2019.

**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**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **E**ff**ects of Human Activities on Hydrological Components in the Yiluo River Basin in Middle Yellow River**

#### **Xiujie Wang 1, Pengfei Zhang 1,2, Lüliu Liu 3,\*, Dandan Li <sup>1</sup> and Yanpeng Wang <sup>1</sup>**


Received: 4 March 2019; Accepted: 29 March 2019; Published: 3 April 2019

**Abstract:** Land use and land cover change (LUCC) and water resource utilization behavior and policy (WRUBAP) affect the hydrological cycle in different ways. Their effects on streamflow and hydrological balance components were analyzed in the Yiluo River Basin using the delta method and the Soil and Water Assessment Tool (SWAT). The multivariable (runoff and actual evapotranspiration) calibration and validation method was used to reduce model uncertainty. LUCC impact on hydrological balance components (1976–2015) was evaluated through comparison of simulated paired land use scenarios. WRUBAP impact on runoff was assessed by comparing natural (simulated) and observed runoff. It showed that urban area reduction led to decreased groundwater, but increased surface runoff and increased water area led to increased evaporation. LUCC impact on annual runoff was found limited; for instance, the difference under the paired scenarios was <1 mm. Observed runoff was 34.7–144.1% greater than natural runoff during November–June because of WRUBAP. The effect of WRUBAP on wet season runoff regulation was limited before the completion of the Guxian Reservoir, whereas WRUBAP caused a reduction in natural runoff of 21.6–35.0% during the wet season (July–October) after its completion. The results suggest that WRUBAP has greater influence than LUCC on runoff in the Yiluo River Basin. Based on existing drought mitigation measures, interbasin water transfer measures and deep groundwater exploitation could reduce the potential for drought attributable to predicted future climate extremes. In addition to reservoir regulation, conversion of farmland to forestry in the upstream watershed could also reduce flood risk.

**Keywords:** changes in hydrological components; effects of human activities; LUCC; WRUBAP; Yiluo River

#### **1. Introduction**

The climate and human activities are two factors that can affect the hydrological cycle in different ways. Climate change alters hydrological systems by inducing both spatiotemporal variations of regional precipitation and changes in temperature [1–3]. Compared with climate change, human activities are more controllable; thus, alteration of human activities constitutes the principal measure for dealing with the potential impacts of climate change on hydrological systems. With recent developments of society and technology, human activities have gradually increased and their consequential impacts on the hydrological cycle on different spatiotemporal scales, such as river basins, have become widely recognized. In general, human activities in river basins can be divided into land use and land cover change (LUCC), e.g., vegetation degradation, deforestation, and urbanization, and water resource

utilization behavior and policy (WRUBAP), e.g., agricultural irrigation [4], reservoir regulation [5–7], deep groundwater extraction, and interbasin water diversion.

Numerous studies have shown that LUCC can have considerable impact on watershed hydrology in terms of water quantity and water quality. For example, vegetation degradation causes decline in the water storage capacity, which ultimately causes changes in the hydrological balance components (e.g., evapotranspiration, infiltration, and baseflow), or vice versa [8,9]. Moreover, changes in surface albedo, surface aerodynamic roughness, leaf area, and rooting depth caused by changing land use can influence the hydrological cycle via different processes. Land use conflicts have impacts on concentration or yield of nitrates, phosphorus, and sulphates, etc. in both surface water and ground water [10–12]. Hence, understanding the hydrological response of a watershed to LUCC constitutes an important step toward sustainable water resources management. Reasonable WRUBAP is another focus of global attention, particularly in semiarid and subhumid agricultural regions where appropriate watershed management is extremely important in relation to the prevention of droughts/flooding. WRUBAP can affect runoff directly via water intake, water transfer, and reservoir operation and indirectly via other components of the water balance, e.g., groundwater, actual evapotranspiration (ET) and infiltration. The mechanisms via which LUCC and WRUBAP affect hydrological processes are different. Therefore, it is highly important to separate the impacts of LUCC and WRUBAP for sustainable water resources utility and water management, which could ultimately ensure manageable agricultural development and ecological environment protection [13–15].

The paired catchments approach [16–18], statistical analyses, and other modeling approaches [19,20] are methods commonly adopted to further the understanding of how human activities might influence basin hydrology. The paired catchments approach has certain limitations attributable to the sizes and characteristics of the watersheds. Analyses based on statistical methods represent a data-driven approach and thus are not based on physical processes. Distributed hydrological models are physically based models that have been used widely to simulate hydrological responses to human activities [13,21,22]. The Soil and Water Assessment Tool (SWAT) [23–25], which is one of the models used most commonly in hydrological response studies, can be used to assess the influence of human activities on the hydrological balance. The delta approach is often applied in association with the SWAT both to analyze the impacts of LUCC through comparison of hydrological responses under two different land use scenarios for the same time frame, and to evaluate the impacts of WRUBAP by comparing observed values with simulated values that are considered representative of natural runoff.

The Yiluo River Basin in China is located in a region where a semiarid area borders a subhumid area. It supplies water for an important grain-producing region in Henan Province. Regional agricultural production and food security are frequently threatened by serious flooding and droughts attributable to the monsoon climate. Increasingly severe problems regarding water resources are expected in the future because drier springs and increasingly severe flooding over long return periods (25 and 50 years) have been projected under the 1.5 and 2 ◦C global warming scenarios [26]. However, changes in regional land use have been detected since the 1990s, e.g., the area of cultivated land has decreased and the area of urbanized land has increased [27]. The effects of such changes in land use on hydrological processes have been investigated in previous studies [21,28]. For example, the change characteristics of intra-annual and interannual runoff, as well as the relative contributions of climate change and human activities on the decrease in annual runoff, have been explored in earlier research [29,30]. However, their effects on other hydrological variables have rarely been investigated.

This study had two primary objectives: (1) to assess the impacts of LUCC on hydrological balance components, and to determine the relative contributions of individual land use types and (2) to analyze the impact of WRUBAP on streamflow during different periods. The remainder of this paper is organized as follows. Section 2 provides an overview of the study river basin, data, and methods. Section 3 presents the results of both the calibration and validation of the hydrological model and the responses of the hydrological components to LUCC and WRUBAP. The approaches adopted to

reduce the uncertainty of model parameters and the effects of LUCC and WRUBAP on the hydrological processes are discussed in Section 4. Finally, the main conclusions are summarized in Section 5.

#### **2. Materials and Method**

#### *2.1. Study Area*

The Yiluo River Basin (33◦–35◦ N, 109◦–113◦ E) covers an area of around 18,563 km<sup>2</sup> (Figure 1). The river, which is 449 km long, runs from Shaanxi Province through Henan Province and into the Yellow River in the city of Luoyang. The Heishiguan hydrological station is located at the outlet of the Yiluo River in the northeast of the basin. The Yiluo River has two principal tributaries: the Yihe River and the Luohe River, on each of which is a large reservoir. The Lushun Reservoir, which is located in the middle reaches of the Yihe River, was completed in 1965. The Guxian Reservoir, located in the middle reaches of the Luohe River, was completed in 1994. There are ten types of soil within the Yiluo River watershed: Acrisols, Alisols, Andosols, Arenosols, Anthrosols, Cambisols, Fluvisols, Gleysols, Solonchaks, and Vertiaols (Figure 2). The major land use classes of the region are cultivation (AGRR), barren (BARR), forest (FRST), grassland (RNGE), urban (URBN), and water (WATR) (Figure 3).

**Figure 1.** Study area with weather stations, reservoirs, and watershed outlets.

Dominated by a typical continental monsoon climate, precipitation is concentrated in the wet season, and about 60% of the total precipitation falls during June–September [30]. The annual mean temperature is about 12 ◦C, and the annual averaged precipitation is in the range 414–1066 mm (1960–2016). The fluctuation of precipitation causes corresponding changes in the hydrological components. The main source of water supply in the basin is local surface runoff, shallow groundwater and deep groundwater. During the 1980s, local surface runoff accounted for approximately 58–69% of the total. Local surface runoff is lower in dry years than in normal years, while water demand increases in dry years because of the requirements of agriculture; consequently, such conditions lead to water shortages [31]. Generally, water consumption increases with socioeconomic development and

an increasing tendency has been observed during the past 20 years in China from Water Resources Bulletins from 1997–2017 (www.mwr.gov.cn).

**Figure 3.** Map of land use maps in 2010 in the Yiluo River Basin (maps for 1980, 1990, and 2000 not shown).

#### *2.2. Dataset*

A digital elevation model (DEM) maps of soil types and properties, LUCC, meteorological variables, and river discharge data were used in this study. The digital elevation model with 90-m resolution was obtained from the Geospatial Data Cloud of China (http://www.gscloud.cn). The LUCC maps with 1-km resolution for 1980, 1990, 2000, and 2010 were acquired from the Resource and Environment Data Cloud Platform (http://www.resdc.cn). These were used to analyze the changes in land use and to investigate the effects of such changes on the water balance. Soil maps with 1-km resolution were obtained from the Harmonized World Soil Database v1.1 (http://www.fao.org/soils-portal/soil-survey/en). Soil parameters used in the hydrological model were estimated using the Soil–Plant–Air–Water budgeting tool.

Daily minimum temperature, maximum temperature, precipitation, relative humidity, wind speed, and sunshine duration at 12 meteorological stations and daily evaporation at the Guxian hydrological station during 1969–2015 were obtained from the National Meteorological Information Center of the China Meteorological Administration, as shown in Figure 1. Evaporation was monitored using E601B pan.

Monthly data of observed streamflow through the Heishiguan hydrological station during 1971–2015 were obtained from the Hydrological Yearbook of the Yellow River. These data were used both to calibrate and validate the hydrological model and to analyze the effects of WRUBAP on runoff.

#### *2.3. Hydrological Modeling*

#### 2.3.1. Model Introduction

SWAT [23,32] is a semidistributed hydrological model that is used widely to simulate the effects of varying soils, land use, and management practices on hydrological conditions, sediment loading, and contamination on the basin scale [32,33]. This model was applied in this study to the Yiluo River Basin to assess the impacts of LUCC and WRUBAP on hydrological balance components: surface runoff, evapotranspiration, return flow, and lateral flow. The daily water-balance equation can be expressed as follows:

$$SW\_t = SW\_0 + \sum\_{i=1}^{t} \left( R\_{day} - Q\_{surf} - ET - W\_{sep} - Q\_{\mathcal{G}\mathcal{W}} \right) \tag{1}$$

where *SWt* is the final soil moisture content on day *t* (mm), *SW*<sup>0</sup> is the initial soil moisture content (mm), *Rday* is the amount of precipitation on day *i* (mm), *Qsurf* is the amount of surface runoff on day *i* (mm), *ET* is the amount of actual evapotranspiration on day *i* (mm), *Wseep* is the amount of water transferred from the soil profile into the gas zone on day *i* (mm), and *Qgw* is the return flow of day *i* (mm).

#### 2.3.2. Model Setup

The Yiluo River Basin was divided into 29 sub-basins (Figure 4) and the hydrological response units were based on the land uses, soil types, and slope classes. The modified Soil Conservation Service Curve Number method [34,35], Penman–Monteith method [36,37], and Muskingum routing method were used to simulate the hydrological components.

**Figure 4.** Sub-basins within the study area.

#### 2.3.3. SWAT Calibration and Validation

The SWAT model was calibrated using observed monthly runoff and evaporation data series from 1971–1982, and it was validated using monthly historical data from 1983–1985. The 1985 cutoff date was based on the recognition that natural runoff has been changed by human activities since 1986, because a change point in the annual runoff series was detected around 1986, whereas no change point was observed in the annual precipitation series [30].

To reduce model parameter uncertainty and equifinality, the multivariable (runoff and actual evapotranspiration) calibration and validation method (MCVM) was used in this study. Parameters were auto-calibrated using Particle Swarm Optimization (PSO) [38] in the SWAT-CUP (SWAT Calibration and Uncertainty Programs) package [39], based on monthly observed discharge through the Heishiguan hydrological station and the set of 20 initial sensitive parameters listed in Table 1. The initial two years of the simulated outputs were used as a warm-up period [40]. After autocalibration, the performance of ET was improved further by manual adjustment of EPCO (Plant uptake compensation factor) and ESCO (Soil evaporation compensation factor), which are parameters related to soil evaporation and plant transpiration [41]. The formula for calculation of ET is as follows:

$$ET = K \times ET\_{\text{pm}} \tag{2}$$

where *ET* is the actual evapotranspiration (mm), *ETpan* is the observed evaporation monitored using an E601B pan and *K* is a conversion coefficient of evaporation, the value of which is 0.81 according to [42].

Three indices including the Nash–Sutcliffe efficiency (*NSE*) [43], coefficient of determination (*R2*), and percent bias (*PBIAS*) were used to evaluate the model's performance. The performance of the SWAT model can be judged satisfactory if *R2* and *NSE* are both >0.5 and *PBIAS* is <±25% [44].

$$NSE = 1 - \frac{\sum\_{i=1}^{n} \left( Q\_i^{\text{sim}} - Q\_i^{\text{obs}} \right)^2}{\sum\_{i=1}^{n} \left( Q\_i^{\text{obs}} - \overline{Q^{\text{obs}}} \right)^2} \tag{3}$$

$$\mathcal{R}^2 = \frac{\left[\sum\_{i=1}^n \left(Q\_i^{obs} - \overline{Q^{obs}}\right) \left(Q\_i^{sim} - \overline{Q^{sim}}\right)\right]^2}{\sum\_{i=1}^n \left(Q\_i^{obs} - \overline{Q^{obs}}\right)^2 \sum\_{i=1}^n \left(Q\_i^{sim} - \overline{Q^{sim}}\right)^2} \tag{4}$$

$$PBIAS = \frac{\sum\_{i=1}^{n} \left(Q\_l^{\text{obs}} - Q\_l^{\text{sim}}\right) \times 100}{\sum\_{i=1}^{n} \left(Q\_l^{\text{obs}}\right)} \tag{5}$$

where *n* is the total number of sample pairs, *Qi* obs is the observed value, *Q*obs is the mean of the observed values, *Qi* sim is the simulated value and *Q*sim is the mean of the simulated values.


**Table 1.** Optimal parameter values calibrated for the Yiluo River Basin.

v\_\_ means the current parameter value is to be replaced by a given value, a\_\_ means a given value is added to the current parameter value, and r\_\_ means the current parameter value is multiplied by (1 + a given value). \*means the optimal parameter value was finally determined by manual calibration based on actual evaporation (*ET*).

#### *2.4. Evaluation E*ff*ects of LUCC and WRUBAP*

The delta approach was applied to assess the impact of LUCC. First, four scenarios (S1, S2, S3, and S4) were defined based on available LUCC data for 1980, 1990, 2000, and 2010 and the contemporaneous climate, which was taken as the ten-year mean during 1976–1985, 1986–1995, 1996–2005, and 2006–2015, respectively. Among them, S1 was set as the baseline period. The runoff simulated under the four actual LUCC scenarios approximately emulated the natural runoff on the monthly scale. Then, three LUCC scenarios (S2\*, S3\*, and S4\*) were defined using the same climate series as in S2, S3, and S4 and LUCC data for 1980 (Table 2). The effects of LUCC were evaluated by comparing S2 against S2\*, S3 against S3\*, and S4 against S4\* using the following Equation:

$$\text{Deviation}\_{\text{lucc}} = \mathcal{W}\_{\text{sim}}^{\text{S}} - \mathcal{W}\_{\text{sim}}^{\text{S\*}} \tag{6}$$

where Deviationlucc is the change of the water balance term under the actual scenario relative to the baseline scenario, *WS*<sup>∗</sup> *sim* is the water balance term under the baseline scenario and *<sup>W</sup><sup>S</sup> sim* is the water balance term under the actual scenario.

The effects of WRUBAP on river runoff were evaluated by comparing the simultaneously simulated and observed runoff under each of the S1, S2, S3, and S4 scenarios using the following Equation:

$$\text{Deviation}\_{\text{WRUBAP}} = Q\_{\text{sim}} - Q\_{\text{abs}} \tag{7}$$

where DeviationWRUBAP is the difference between the simulated and observed runoff, *Qsim* is the simulated runoff under the actual LUCC and *Qobs* is the observed runoff.


**Table 2.** Scenario setting based on climate and land use and land cover change (LUCC).

#### **3. Result**

#### *3.1. Model Calibration and Validation*

The *NSE* value for monthly runoff during the calibration periods was 0.85, the *R<sup>2</sup>* value was 0.88 and the *PBIAS* value was 13.37. The *NSE* for yearly *ET* during calibration was 0.84, the *R<sup>2</sup>* value was 0.87 and the *PBIAS* value was −1.36, indicating good agreement between the monthly simulated values and the observation values during the calibration period. Moreover, satisfactory results were also obtained during the validation, with *NSE*, *R2*, and *PBIAS* values of 0.72, 0.75, and 12.01, respectively, in the runoff validation and 0.80, 0.81, and −6.06, respectively, in the *ET* validation. The results of the peak value were reasonably accurate, but the base flow was underestimated (Figure 5a). Similarly, the peaks of *ET* > 150 mm/month were also underestimated (Figure 5b). Overall, the simulations for water discharge and *ET* indicated satisfactory agreement between the observed and simulated values; thus, the parameters were considered credible and the model deemed suitable for simulating the natural water balance of the Yiluo River Basin.

**Figure 5.** Simulated and observed monthly (**a**) discharge and (**b**) actual evapotranspiration (*ET*) during calibration period (1971–1982) and validation period (1983–1985).

#### *3.2. Analysis of LUCC*

The areas of AGRR, FRST, RNGE, URBN, WATR, and BARR in 1980, 1990, 2000, and 2010, as well as their changes relative to 1980, are listed in Table 3. In 1980, the main land use types were AGRR, RNGE, and FRST, which accounted for 94.8% of the total watershed area. Through urban centralization and the "Returning farmland to grass" policy, sporadic areas of URBN were converted to AGRR, while some AGRR land was converted to RNGE during 1980–1990. The RNGE area increased by 115 km2 from 3038 km2, but the area of AGRR, URBN, and FRST decreased by 59, 48, and 7 km2, respectively, in 1990 relative to 1980. After 1990, because of reservoir construction and enhanced urbanization, the areas of WATR and URBN both increased, while the areas of the other four main land use types decreased. The WATR area increased by 35 km<sup>2</sup> in 2000 and by 38 km<sup>2</sup> in 2010, URBN increased by 121 km<sup>2</sup> in 2000 and by 187 km<sup>2</sup> in 2010, but AGRR decreased by 72 and 148 km2 in 2000 and 2010, respectively, relative to 1980.

**Table 3.** LUCC in 1980, 1990, 2000, and 2010 relative to 1980.


#### *3.3. Responses of Hydrological Components and Runo*ff *to LUCC*

Table 4 summarizes the annual rainfall, *ET*, surface runoff (SURQ), return flow (GWQ), and runoff depth (RD) under each scenario and their difference between paired scenarios. Comparatively, the values of ET and GWQ are higher by 0.2 and 0.1 mm but SURQ and RD are lower by 0.3 and 0.4 mm, respectively, under scenario S2 relative to S2\*. The values of ET, SURQ, and RD are higher by 1.0, 0.5, and 0.2 mm, respectively, and GWQ is lower by 0.5 mm under scenario S3 relative to S3\*. The values of *ET* and SURQ are higher by 0.9 and 0.7 mm, respectively, and GWQ is lower by 0.3 mm under scenario S4 relative to S4\*.


**Table 4.** Water balance components and the difference between paired scenarios.

Generally, ET occurs in three forms: vegetation transpiration (VT), soil evaporation (SE), and water evaporation (WE) [45]. In the efficient consumption of ET, vegetation land uses (AGRR, FRST, and RNGE) possess greater VT and SE than URBN land use. The total amount of WE from WATR is much higher than from the other five LUCC types. Under scenario S2, the increase in ET was caused mainly by the reduction of the URBN area and the expansion of RNGE. It also caused a slight decrease of RD. The URBN and RNGE under scenarios S3 and S4 presented an opposite trend to scenario S2 However, the trend of increase in ET was enhanced, which could be explained by the obvious expansion of the WATR area. In addition, plant root systems change the physical properties of soil, improve infiltration and GWQ and ultimately decrease SURQ. This suggests that the expansion of vegetation land uses leads to increase of both GWQ and ET but to decrease of SURQ, which consequently leads to higher water-holding capacity. These effects can be seen under scenario S2 relative to S2\*. The effects of URBN expansion are opposite to those of vegetation expansion because URBN does not allow much soil evaporation or water infiltration; thus, most precipitation within an URBN area flows into rivers in the form of SURQ. The continuous urbanization after 1990 in the study area has played a positive role in the increase of SURQ.

#### *3.4. Response of Runo*ff *to WRUBAP*

This study analyzed the effect of WRUBAP on RD only because of data availability. In this part, the simulated RD under the baseline scenario (S1) and the actual LUCC scenarios (S2, S3, and S4) are assumed to represent natural RD.

Scenario S1 (1976–1985): At Hesihiguan station, the observed monthly runoff under scenario S1 was obviously higher than the natural runoff in the first seven months of the year but lower in the final three months (Figure 6a). The main cause was reservoir regulation performed to meet the needs of agricultural irrigation during April–June. Consistency between the dry season and the irrigation period in the study area requires reservoir regulation. Water storage in the reservoir during October–December in most cases caused the observed runoff to be lower than the natural runoff during this period. To supply crop growth with sufficient water during April–June, the watershed management commission would release reservoir storage. That action resulted in observed runoff (78.1 mm) being higher than the natural RD (56.7 mm) during the dry season, and WRUBAP accounted for about 37.7% of the

natural RD (Table 5). In addition, because of flood prevention, reservoir storage capacity was reserved by discharging water in the early wet season; therefore, observed runoff was slightly higher than natural runoff during July. Ultimately, the effect of WRUBAP on runoff was reasonably small during the wet season under scenario S1.

**Figure 6.** Comparison between the simulated and observed annual runoff values for scenarios (**a**) S1, (**b**) S2, (**c**) S3, and (**d**) S4.


**Table 5.** Comparison between the simulated (natural) and observed mean discharge during the dry season and the wet season under scenarios S1, S2, S3, and S4.

Notes: Deviation = observation − simulation; Deviation rate = Deviation/simulation×100%. Dry season: January–June, November and December, Wet season: July–October.

Scenario S2 (1986–1995): Observed annual runoff under scenario S2 was greater than natural runoff in most years except 1992 and 1995, as shown in Figure 6b. Moreover, the observed monthly runoff was greater than the natural runoff in both the dry season and the early wet season (July), as shown in Figure 7b. This might be attributable to reservoir regulation and to irrigation through the extraction of deep groundwater. Storing water in the reservoir upstream of the Heishiguan station during August and September led to the observed runoff being smaller than the natural runoff. Conversely, discharging water during the dry season led to the observed runoff being greater than the natural runoff. The extraction of deep groundwater to supply irrigation increased the base flow, which ultimately led to the observed RD (72.4 mm) being higher by 144.1% than the natural RD (29.7 mm) during the dry season (Table 5). Consequently, the impact of WRUBAP on runoff in the wet season under scenario S2 was reasonably small.

**Figure 7.** Simulated and observed monthly runoff under scenarios (**a**) S1, (**b**) S2, (**c**) S3, and (**d**) S4.

Scenario S3 (1996–2005): As shown in Figure 6c, annual observed runoff was smaller in the relatively wet years (1996, 1998, 2000, 2003, and 2005) but greater in the relatively dry years (1997, 1999, and 2002) than natural runoff. This might be due to stronger interannual water resource management than under scenarios S1 and S2 to prevent floods during the wet years and to prevent droughts during the dry years. The effects of regulation can also be seen in Figure 7c and Table 5. The observed RD (165.2 mm) was smaller by 35.0% than the natural RD (254.2 mm) during the wet season, but the observed RD (63.3 mm) was higher by 34.7% than the natural RD (47.0 mm) during the dry season. It indicates that the effect of WRUBAP on runoff under scenario S3 was enhanced during the wet season.

Scenario S4 (2006–2015): Similar to scenario S3, observed annual runoff was smaller in relatively wet years but greater in relatively dry years (Figure 7d) than natural runoff. Natural annual RD under scenario S3 (116.7 mm) was higher by 27.9 mm than under scenario S4 (88.8 mm). However, annual observed RD under scenario S3 (97.3 mm) was almost equal to that under scenario S4 (95.4 mm) (Table 5). Observed monthly RD (145.5 mm) was lower by 21.6% than natural RD (185.7 mm) during the wet season but higher by 34.7% than natural RD (40.8 mm) during the wet season. It suggests the effect of WRUBAP under scenario S4 alleviates the unbalanced temporal distribution of the water resource.

#### **4. Discussion**

Parameter uncertainty and equifinality can influence numerical model performance. At the basin scale, multiobjective optimization algorithms [46–48], multiperiod parameterization [49,50], and the MCVM are methods commonly adopted to reduce uncertainty. Generally, the MCVM is based on observed discharge at multiple hydrological stations. In this study, the two variables of discharge and *ET* were used to calibrate/validate a numerical model to reduce parameter uncertainty and to improve *ET* performance. Although human activities, especially WRUBAP, were not accounted for in the calibration and validation of SWAT model, the model performance is satisfactory. Simulation, not observation, during the baseline period was compared with other scenarios to explore the effects of human activities.

In this study, the CN value of FRST was lower than other five land use types, which indicates that FRST has a stronger ability of water storage and is likely to reduce more SURQ. It has previously been found that Grassland and cultivation have similar influences on hydrological balance components in this area [14], and the mutual transformation between them has little effect on water yield within the Yiluo River basin. Compared with herbaceous plants (grassland and cultivation), woody plants have more developed root systems and canopies; thus, conversion of herbaceous plants into forests would enhance the capacity of catchment to conserves soil and water, and reduce water yield [51,52], which consequently affects the yield of phosphorus and nitrates [12,53]. Considering that the upstream regions of the Yiluo River Basin are mountainous areas and are not the main crop region, returning farmland to forests is conducive to preventing land degradation, increasing infiltration, and thus

reducing flood risk [54,55]. However, quantitative analysis has not been performed in this study. This will be carried out in further study.

Hydrological changes have important impacts on the local use of natural water [56]. However, compared with LUCC, WRUBAP has more dynamic and greater impact on hydrological processes, particularly in monsoon climate zones. It is noticed that the Guxian Reservoir greatly improved the capacity of runoff regulation in the wet season. Therefore, reservoirs could be used to respond to climate variety, which increases runoff in the dry season while decreasing runoff in the wet season. A recent study made similar conclusions, indicating that large and small reservoirs can have equally large effects on runoff [57].

#### **5. Conclusions**

In this study, the SWAT model and the MCVM were used to assess the impact of human activities (LUCC and WRUBAP) on the catchment hydrology of the Yiluo River Basin during 1976–2015. Based on the analysis, the following conclusions were drawn.


**Author Contributions:** Conceptualization, X.W. and L.L; data curation, P.Z.; Methodology L.L. and P.Z.; Resources, X.W. and L.L.; Software, P.Z. and L.L.; Supervision, Z.W. and L.L; Writing—Original Draft, L.L., X.W. and P.Z.; Manuscript revision/review D.L. and Y.W.

**Funding:** This study was supported by the National Key R&D Program of China (Grant Nos. 2017YFA0605004, 2016YFE0102400, and 2018YFC1508403), and National Natural Science Foundation of China (51579173).

**Acknowledgments:** The authors wish to acknowledge the assistance of all the editors and reviewers. We thank James Buxton MSc from Liwen Bianji, Edanz Group China (www.liwenbianji.cn./ac), for editing the English text of this manuscript.

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

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Flood Vulnerability, Environmental Land Use Conflicts, and Conservation of Soil and Water: A Study in the Batatais SP Municipality, Brazil**

**Anildo Monteiro Caldas 1, Teresa Cristina Tarlé Pissarra 2, Renata Cristina Araújo Costa 2, Fernando Cartaxo Rolim Neto 1, Marcelo Zanata 3, Roberto da Boa Viagem Parahyba 4, Luis Filipe Sanches Fernandes <sup>5</sup> and Fernando António Leal Pacheco 6,\***


Received: 12 September 2018; Accepted: 28 September 2018; Published: 29 September 2018

**Abstract:** In many regions across the planet, flood events are now more frequent and intense because of climate change and improper land use, resulting in risks to the population. However, the procedures to accurately determine the areas at risk in regions influenced by inadequate land uses are still inefficient. In rural watersheds, inadequate uses occur when actual uses deviate from land capability, and are termed environmental land use conflicts. To overcome the difficulty to evaluate flood vulnerability under these settings, in this study a method was developed to delineate flood vulnerability areas in a land use conflict landscape: the Batatais municipality, located in the State of São Paulo, Brazil. The method and its implementation resorted to remote sensed data, geographic information systems and geo-processing. Satellite images and their processing provided data for environmental factors such as altitude, land use, slope, and soil class in the study area. The importance of each factor for flood vulnerability was evaluated through the analytical hierarchy process (AHP). According to the results, vast areas of medium to high flood vulnerability are located in agricultural lands affected by environmental land use conflicts. In these areas, amplified flood intensities, soil erosion, crop productivity loss and stream water deterioration are expected. The coverage of Batatais SP municipality by these vulnerable areas is so extensive (60%) that preventive and recovery measures were proposed in the context of a land consolidation–water management plan aiming flood control and soil and water conservation.

**Keywords:** flood vulnerability; land use conflicts; water management; soil conservation; spatial multi-criteria analysis; geographic information system

#### **1. Introduction**

Natural floods occur when stream flow exceeds the river channel's capacity, and consequently spills over the natural or artificial embankments. Floods are particularly related with extreme precipitation events and specific geomorphologic settings. When floods reach the marginal urban and rural areas they can cause significant economic damage and even human deaths [1,2].

Land uses, which are forms of land occupation and refer to the management and modification of a natural environment, influence the intensity of processes such as infiltration and runoff. The areas of greater impermeability, such as settlements, favor the runoff process and, hence, cause a substantial accumulation of rainwater in drainage channels, especially during periods of intense rainfall. The same holds for intensively tilled grounds, such as arable fields or pastures, as well as for managed woods. Therefore, the spatial distributions of urban and agricultural areas at the neighborhood of a river channel, and especially their dynamic changes over time can be the cause of frequent and often unexpected overflows with crop and property losses [3–12].

The frequency and magnitude of floods can increase where environmental land use conflicts have developed. In general, a land use conflict occurs when there are conflicting views on land use policies, such as when an increasing population creates competitive demands for the use of the land, causing a negative impact on other land uses nearby [13]. However, a different definition has been presented to describe land use conflicts in rural watersheds mostly occupied by agriculture, pastures, and forests [14,15]. In that context, an environmental land use conflict develops where the current land use differs from a natural use set up on the basis of specific morphometric parameters, namely drainage density and hill slope. In the same watersheds, conflicts were also associated with changes in permanent preservation areas, such as destruction of riparian vegetation. The consequences of conflicts for floods are apparent. For example, in the long term the conversion of forest into agricultural land reduces soil permeability, because the use of heavy machines over many years increases compaction of topsoil layers and deteriorates the soil structure [4,16,17]. The stability of the soil influences its infiltration capacity [3,4] and its ability to preserve structure under the impact rain drops [18,19]. The overall result is increased runoff and soil erosion. Thus, flood prevention as a priority measure should be emphasized to avoid environmental land use conflicts. To accomplish this goal, it is important to improve design solutions for natural drainage by enhancing and promoting innovative forms of land use structures, as well as the correction of riverbeds and specific sectors along the river, especially riparian and floodplain areas [20].

A pre-requisite for the proposition and implementation of flood control solutions and conservation of soil and water measures, is the evaluation and mapping of flood vulnerability and land uses [21,22], besides environmental land use conflicts [14,15] if the focus problem is to assess human-induced impacts on flood prone agro-systems triggered by land use changes. However, if the analysis of flood vulnerability and environmental land use conflicts, taken separately, is relatively abundant, the potential amplifying factor of environmental land use conflicts on flood frequency and intensity, as well as the feedback of amplified floods on agro-systems has barely been discussed. This study aims contributing to fill in this gap. To accomplish this general goal the study established the following specific objectives, using a geographic information system (GIS) to implement them: (1) prepare a flood vulnerability map for the studied area, using multi-criteria analysis based on relevant flood-control attributes (e.g., terrain slope, land use); (2) prepare a map of environmental land use conflicts based on comparisons between land capabilities (natural uses) and actual uses [23]; (3) cross-tabulate the data on flood vulnerability and environmental land use conflicts to verify if flood prone regions extend beyond the areas in environmental land use conflict and estimate the overlapping areas; and (4) Using Pareto diagrams, identify the regions requiring priority action as regards flood control/soil and water conservation and propose preventive or recovery methods.

The studied area comprises the Batatais SP municipality, located in Brazil, which is a flood prone region with severe environmental land use conflicts caused by replacement of natural vegetation with sugar cane plantations. The local population and government have the right and the duty to perceive the risks they are vulnerable to, to plan flood control and adapt land uses in the sequel. It is expected that the identification of areas susceptible to flooding and areas of conflicts, offered in this study, subsidize decision-making in a context of land consolidation–water management.

#### **2. Material and Methods**

#### *2.1. Study Area*

The study area is located in the northeast of São Paulo State, Brazil (Figure 1), in the municipality of Batatais, at the latitude 20◦53 28 S, longitude 47◦35 06 W of Greenwich, and average altitude of 862 meters. The Batatais municipality covers an area of approximately 851 km2 and has a population of 56,476 inhabitants [24]. The climate is considered Cwa in the classification of Köppen-Geiger, described as tropical (mild) with dry winter. The rainy period runs from November to March. The average precipitation (period 1970–2017) is 1972 mm/year and 100.4 mm/day. The average annual temperature is 21.9 ◦C (http://www.inmet.gov.br). In the past century precipitation has steadily increased on the annual and daily basis. Both trends are displayed in Figure 2 and point to an average increase of 10 mm/year and 0.7 mm/day. Within this period, most years after 1995 were characterized by precipitation above the average, while before that most years received rainfall below the average. Relief comprises the structural plateaus of Ribeirão Preto, characterized by 500–700 m of altitude and shallow incised stream valleys, and the residual plateaus of Franca/Batatais installed at higher altitudes (800–1100 m) and characterized by deep incised valleys (http://www.abagrp.cnpm.embrapa.br). According to the Prefeitura Municipal da Estância Turística de Batatais (https://www.batatais.sp.gov.br), in 2014 land use was largely dominated by sugar cane plantations (49,065 hectares; 57.7% of municipality area) and brachiaria pastures (10,437 hectares; 12.3%). Other agriculture activities were much less represented, such as coffee (2625 hectares; 3.1%), corn (2247 hectares; 2.6%) or soya (1349 hectares; 1.6%). The predominant type of vegetation cover is the Cerrado, composed of forest savanna, wooded savanna, park savanna and gramineous-woody savanna. The Batatais municipality hosts a Conservation Unit called the State Forest of Batatais, which occupies an area of approximately 1353.3 hectares (1.6%). Other forested areas are occupied by pine (976 hectares; 1.1%) or eucalyptus (419 hectares; 0.5%), while the urban centers occupy 751 hectares (0.9%). The dominant soil type is distroferric red latosols with sandy to sandy-loam texture.

**Figure 1.** Study area, municipality of Batatais, São Paulo State, Brazil.

**Figure 2.** Precipitation records of Batatais State Forest Station relative to period 1970–2017 (red and blue lines), and associated increase trends (dashed black lines). The labels close to blue line peaks refer to precipitations at reported major flood events. The daily precipitation of 174 mm occurred in 29 January 2008 and provoked a severe flood in Batatais.

The scientific studies on floods are not abundant for the Batatais municipality, but several videos on the consequences of urban and rural flood events can be viewed through Brazilian television websites (e.g., the Globo Tv). Based on comparisons between daily precipitation records (e.g., Figure 2) and reported flood events, a precipitation >140 mm/day can be used as reference for the occurrence of a severe flood in Batatais. Figure 2 points to the occurrence of four major events in the 1970–2017 period (labeled peaks), which represents approximately one major event per decade. Figure 3 illustrates the severe flood occurred in 29 January 2008, when the daily precipitation reached 174 mm, the largest value registered in the Batatais Forest Station in the 1970–2017 period. Besides the four events highlighted in Figure 2, some local newspapers have reported the occurrence of floods in 20 January 2014, 9 January 2016, or 29 December 2016, presumably at lower daily precipitations. This is indication that moderate floods can occur at rather short return periods (e.g., 1–2 years).

**Figure 3.** Illustration of 29 January 2008 severe flood in the Batatais town.

#### *2.2. Raw Data and Software*

The image dataset comprises radar images from the SRTM (Shuttle Radar Topography Mission), with a spatial resolution of 30 m, referring to the WGS84-Zone 23K (SF-23-VA and SF-23-VC); and a mosaic of images from the Landsat 8 optical sensor, with a spatial resolution of 30 m and temporal resolution of 16 days, referred to Datum WGS84. These datasets were used to prepare all the maps and to extract data on environmental attributes. The images were obtained using public data sources

and processed in ArcGIS software from ESRI, licensed to Universidade Estadual Paulista (UNESP), Departamento de Engenharia Rural, Campus Jaboticabal. The public data sources comprised the websites of EMBRAPA—Empresa Brasileira de Pesquisa Agropecuária (http://www.relevobr.cnpm. embrapa.br) and of INPE—Instituto Nacional de Pesquisas Espaciais (http://www.dgi.inpe.br/CDSR).

#### *2.3. Preparation of Workable Data*

The flood risk map was calculated and drawn for the municipality of Batatais SP, based on four variables, called environmental attributes, represented by altitude, slope class, land use/occupation and soil type. The maps of all environmental attributes were prepared in raster format as described in the next paragraph.

The altitude and slope classes were developed from a 30 m resolution digital elevation model (DEM), following Brazilian principles and rules on cartography. The DEM images were first corrected for relief depressions and then submitted to the slope routine of ArcGIS that calculated and projected the slopes. The map of land use was drawn with four classes: agriculture, pasture, forest and urban. It was generated from the Landsat 8 images using the supervised classification method, whereby statistical algorithms are run for the recognition of spectral patterns that subsequently are labeled as land uses. The map of soil units was obtained from [25] in paper format and then scanned, geo-referenced, vectorized, and finally transformed into the required raster format. The four maps are illustrated in Figure 4.

**Figure 4.** Maps of environmental attributes within the municipality of Batatais SP.

#### *2.4. Flood Vulnerability*

The mapping of flood vulnerability involved a weighted combination of the environmental attribute maps (Figure 4) operated in ArcGIS using map algebra tools. The selection of attributes was based on the judgment of experts who know the region and could empirically relate historical flood events to these variables. The processing of attributes to produce the map of vulnerable areas followed a conventional multi-criteria analysis (MCA), whereby the original variables are first harmonized into a common dimensionless scale (rated), then are allocated a specific importance (weight), and finally are aggregated into a global vulnerability index. The MCA is used in many environmental applications [26,27].

The first step in the MCA comprised the reclassification of environmental attributes into a dimensionless scale ranging from 1 to 10, where "1" means least susceptible to floods and "10" most susceptible to floods. This common range allows equivalence among attributes that otherwise would not be comparable given the different scales. The range 1–10 is arbitrary but is regularly used in MCA analyses. The reclassification scheme is depicted in Table 1 and resulted from bibliographical surveys and debates within a multidisciplinary team composed of agronomists, forestry engineers, biologists, lawyers, and geographers.


**Table 1.** Reclassification of environmental attributes according to their flood susceptibility. The source for the soil classes is [28].

As regards the altitude attribute, the lower elevations were given the highest susceptibility (rated as 10) because they receive runoff from the higher altitudes (rated as susceptibility 1) and therefore are more prone to flooding. The attribution of susceptibility ratings to the land use attribute was based on potential infiltration capacity of land covers, assessed through curve numbers (https://swat.tamu.edu). The larger the curve number, the larger was rated the susceptibility, because big curve numbers are in favor of big overland flows. In that context, the urban areas were ascribed the largest possible susceptibility (10) because these areas comprise large portions of impermeable surface (e.g., street pavements, commercial and industrial areas), characterized by very large curve numbers (>80 in a maximum of 100). The low susceptibility end-member (rated as susceptibility 1) was the forested areas because vegetation has the capability to retain water in the canopy promoting infiltration and reducing runoff downstream. These areas are frequently characterized by low curve numbers (e.g., between 30 and 60 for forests and bushes located on permeable soils). The occupations by agriculture or pasture were given high susceptibility ratings because characteristic curve numbers in these areas are intermediate to high (60–75). The distribution of susceptibility ratings per slope classes considered the effect of topography in the share of overland flow between exportation and accumulation. The steepest slopes (>75%) were attributed the lowest rating (1) because overland flow in these areas is mostly exported under the action o gravity towards the plains where the surface is leveled promoting deceleration of water velocity, accumulation and floods in the sequel. For the opposite reasons, the low slope areas (slope < 10%) were assigned a high susceptibility rating (8–9). The susceptibility ratings of soil classes took into account the expected permeability of dominant soil types, deduced from soil texture. The textures of red and yellow-red latosols are mostly sandy or sandy-loam, which means the soils are permeable. For that reason, the susceptibility ratings of both soil types were small (3–4).

The analytic hierarchy process AHP; 29 was used to weight the importance of environmental attributes as regards flood susceptibility. This method is based on pairwise comparisons among the attributes according to the fundamental scale of absolute numbers (Table 2). Firstly, a hierarchy of attributes is defined according to expert's judgment. In this study, ranks were assigned to the four environmental attributes as indicated in the last column of Table 1. The altitude attribute was given the largest relative importance (7), because altitude determines where the overland flow is likely to accumulate producing floods. The smallest importance (1) was attributed to soils, because soil types in the studied region are relatively homogeneous as regards permeability. A similar rationale was used to rate the slope class parameter. In this case, the assignment of a rating 3 was related to the observation that the Batatais municipality is generally a flat to undulated region, incised by deep valleys only in specific sectors (e.g., the south limit). An assignment of a large relative importance to these steady attributes would tend to smooth the final susceptibility map masking important contributions from the other attributes. Land use was given the second largest relative importance (5), because the role of vegetation cover is prominent in the control of overland flow, but also because land use/occupation is a dynamic attribute that can change overtime. A large relative importance assigned to this parameter highlights the impact of a land use change, and therefore is likely to improve the efficacy of decision support systems focused on the flood control versus land use management issue.

Fourthly, attribute weights (*w*) are obtained from the main eigenvector of matrix *B* as follows (last column of Table 4):

$$w\_{i} = \frac{\sum\_{j} b\_{ij}}{\sum\_{i} \sum\_{j} b\_{ij}} \text{ with } i \text{ and } j = \text{ 1, 2, ..., } p \tag{2}$$

with:

$$
\sum\_{i} w\_{i} = 1 \text{ with } i = 1, 2, \dots, p
\tag{3}
$$

Finally, a consistency ratio (*CR*) is computed to check whether the elements of matrix *B* have been randomly generated, in which case the attribute weights are not statistical significant and for that reason are invalid within the MCA. The consistency ratio ranges from 0 and 1 and invalidate the weights when *CR* ≥ 0.1 [29]. The details on the calculation of *CR* are beyond the scope of this paper and can be consulted elsewhere [30]. For the weights depicted in Table 4, a *CR* = 0.043 was estimated, which means the weights are valid.

Having validated the weights the vulnerability index is calculated using the following aggregation equation:

$$\text{Fload susceptibility} = 0.06 \times \text{soil} + 0.12 \times \text{slope} + 0.26 \times \text{land use} + 0.57 \times \text{altitude} \tag{4}$$

where the numeric coefficients represent the attribute weights (last column of Table 4).


**Table 2.** The fundamental scale of absolute numbers, according to [29].

Secondly, the ranks are assembled in a pairwise comparison matrix, denoted matrix *A* of *a*ij elements (Table 3). Thirdly, matrix *A* is normalized as matrix *B* with elements *b*ij (Table 4), as follows:

$$b\_{ij} = \frac{a\_{ij}}{\sum\_{i} a\_{ij}} \text{ with } i \text{ and } j = 1, 2, \dots, p \tag{1}$$

**Table 3.** Pairwise comparison matrix based on the attribute ranks listed in the last column of Table 1.


**Table 4.** AHP to obtain the attribute weights. The values under the heading of each attribute represent *b*ij scores calculated by Equation (1). In these columns, the nominator values to the left of the equal sign are *a*ij scores transported from Table 3, while the corresponding denominator values represent the sum of nominator values. The attribute weights (*w*) are depicted in the last column. They result from application of Equation (2) to the *b*ij scores calculated previously.


#### *2.5. Environmental Land Use Conflicts*

Environmental land use conflicts relate to uses of the land that ignore land capability. In this study, land capability was assessed using the ruggedness number concept (*RN = Sl* × *Dd*) based on the morphometric parameters hill slope (*Sl*) and drainage density (*D*d). The *RN* is easily estimated from a digital elevation model (e.g., Figure 4; top-left map) using the adequate terrain modeling tools within the GIS software, namely the slope and density functions of ArcGIS. The relationship between the *RN* and land capability can be described as follows [31]: sectors of a watershed with a low *RN* are assumed capable for the practicing of cropping agriculture as they correspond to undulated low dissected areas. When the *RN* is high the sectors are found proper for an occupation by forests because they are sloping and considerably dismembered. Finally, sectors of the basin with intermediate *RN* values are adequate for livestock pasturing or for a mosaic of natural pastures and forests. The boundaries of each *RN* class are set up in four consecutive steps: (a) a number of sub-basins is identified and delineated across the watershed in study; (b) the mean *Sl* and *D*<sup>d</sup> of each sub-basin are estimated, from which the corresponding *RN* is obtained; (c) the *RN* amplitude (max*RN*–min*RN*) is estimated, where max*RN* is the highest and min*RN* the lowest *RN*, and then divided by a pre-defined number of land capability classes (*n*), for example the four classes used in this study (*n* = 4): 1—Agriculture, 2—Pastures for livestock production, 3—Pastures for livestock production/Forestry and 4—Forestry; (d) the sub-basins are assigned to these classes, based on their individual *RN* scores, from which a land capability map is prepared.

The environmental land use conflict is the divergence of an actual use (e.g., Figure 4; top-right image) from land capability. This conflict can be graduated to express the departure between the natural and actual land uses. The steps required to accomplish this goal were thoroughly described in the study of [15]. Firstly, land capabilities (labeled N) and actual uses (labeled A) were ranked in ascending order of their resemblance with agriculture, and given the codes N = A = 1 (agriculture), 2 (livestock production), 3 (mixed livestock production and forestry), and 4 (forestry). Secondly, the conflict was equated to the difference between codes of capability (1≤ N ≤ 4) and actual land use (1≤ A ≤ 4). According to this method, (a) the no conflict areas are represented by regions where N − A ≤ 0, with the negative values representing areas with potential for a sustainable expansion of agriculture or livestock pasturing; (b) areas where the value of the difference between capability and actual use is equal to 1 or 2 are classified as Class 1 (N – A = 1) or Class 2 (N – A = 2) conflict areas, respectively; (c) areas with potential for forestry (N = 4) but occupied with a farm (A = 1) are referred to as Class 3 (N – A = 3). According to [23], the risks and limitations allocated to conflict classes 1 to 3 are: (Class 1) these lands represent a low risk for surface water contamination and riverine ecosystem degradation when used for annual crops or pastures. The attenuation of the risks may be accomplished by the implementation of soil conservation measures resorting to vegetative and(or) mechanical techniques; (Class 2) these lands represent a medium risk for surface water contamination and riverine ecosystem degradation and, hence, are incapable for the practice of agriculture, although being adapted for livestock pasturing, reforestation or environmental preservation; (Class 3) these lands represent a high risk for surface water contamination and riverine ecosystem degradation and therefore are not capable of being used for cropping or livestock pasturing; nonetheless, they are still capable of being reforested or used for environmental preservation.

#### *2.6. Combined Analysis of Flood Vulnerability, Land Use, and Land Use Conflicts*

The potential impact of floods on agricultural systems, with or without eventual amplification caused by environmental land use conflicts, was assessed through a combined analysis of flood vulnerability, land use and environmental land use conflicts, using Pareto diagrams. A Pareto chart, also called a Pareto distribution diagram, is a vertical bar graph in which values are plotted in decreasing order of relative frequency from left to right. A point-to-point graph, which shows the cumulative relative frequency, may be superimposed on the bar graph. Pareto charts are extremely useful for analyzing what problems need attention first because the taller bars on the

chart, which represent frequency, clearly illustrate which variables have the greatest cumulative effect on a given system. The Pareto chart provides a graphic depiction of the Pareto principle, a theory maintaining that 80% of the output in a given situation or system is produced by 20% of the input (https://whatis.techtarget.com).

#### **3. Results**

The vulnerability to floods in the Batatais SP municipality is represented in Figure 5 and the associated areas quantified in Table 5. The areas with higher vulnerability (levels high and very high) are concentrated in sectors around the municipal boundary where altitudes are lower, namely the west, south, and northeast sectors. The highest vulnerability level (very high, displayed in dark blue in Figure 5) occupies an area of 39.45 km2, which represents 4.64% of municipal territory. These areas are located along the water course margins where the relief is smooth. They represent a serious problem concerning human, structural, or material losses, should an anomalous natural phenomenon occur, namely an intense and very prolonged rainfall event.

**Figure 5.** Flood vulnerability map of Batatais SP municipality.


**Table 5.** Quantification of flood vulnerability areas within the municipality of Batatais SP, as function of vulnerability class. The spatial incidence of vulnerability classes is represented in Figure 5.

The very low (20.31 km2; 2.39%) and low (152.77 km2; 17.95%) vulnerability levels occupy the southwest and southeast sectors of Batatais SP municipality. Vulnerability in these areas is generally low (sometimes very low) because these regions occupy areas of high altitude where accumulation of runoff less likely. The southwest region comprises a particularly low vulnerability sector (the dark brown area), because vegetation cover (Cerrado), permeable forest soil, and low slope gradients all potentiate infiltration increase and runoff decrease, which combined with the high altitudes further reduce flood vulnerability. The medium vulnerability level occupies the largest area (388.92 km2) representing 45.7% of municipal area. This outcome is strongly influenced by the MCA options, whereby altitude and land use were ranked as most relevant factors in the flood vulnerability assessment while agriculture and medium altitudes were ascribed relatively large ratings (Tables 3 and 4). Therefore, a large portion of medium vulnerability areas occupies agricultural areas located at medium altitudes. It is worth mentioning that the urban area of Batatais SP municipality is also described as medium vulnerability area. This is enough reason for civil defense agencies to implement monitoring programs, as well as to adopt prevention measures to minimize the potential losses of floods and the transmission of communicable diseases, such as water-borne (e.g., typhoid fever, cholera) and/or vector-borne (e.g., malaria, dengue) diseases.

Overall, the flood vulnerability analysis reveals a flood-prone municipality where the total area susceptible to floods is much larger than the less susceptible areas, since the regions of medium and high flood vulnerability account for 75.03% of all municipal area. Additionally, the most vulnerable regions, although representing a small area, are located in the vicinity of drainage channels requiring special attention from local authorities.

The map describing environmental land use conflicts within the Batatais SP municipality is illustrated in Figure 6. The percentage of municipal area occupied by conflicts is 65%, distributed as follows per conflict class: class 1—32.3%; class 2—16.4%, and class 3—16.4% (Table 6). The cross-classification of flood vulnerability and environmental land use conflict exposes large percentages of municipality area (78.78%) occupying medium (52.01%), high (23.6%), or very high (3.17%) vulnerability regions where land use is in moderate (class 1; 40.23%), high (class 2; 16.61%) or severe (class 3; 21.94%) conflict with the soils' natural use.

**Figure 6.** Map of land use conflicts within the Batatais municipality.


**Table 6.** Cross-classification of land use conflict and flood vulnerability in the municipality of Batatais SP. The grey shaded cells represent concern situations because they combine medium to very high flood vulnerability and land use conflict.

The illustration of land use conflict-flood vulnerability cross classification in the form of a Pareto diagram is depicted in Figure 7. In this case the Pareto 80/20 rule does not provide an unique solution, because the number of classes representing 20% of input (land use conflict–flood vulnerability classes) is different from the number of classes representing 80% of output (cross-classified area). The high frequency classes representing 20% of input comprise three classes, namely the areas of conflict classes 1 and 3 located in regions of medium flood vulnerability and the areas of conflict class 1 occupying regions of high flood vulnerability. The higher frequencies representing 80% of the cross classified area comprehend the previous three classes plus the areas of conflict class 2 occupying the regions of medium flood vulnerability and the areas of conflict class 1 and 2 located in regions of low flood vulnerability. From a decision-making standpoint, it was assumed that the optimal Pareto solution gathers four classes, namely all regions of medium flood vulnerability affected by land use conflicts (classes 1–3) and the regions of high flood vulnerability affected by class 1 conflict class. Under the Pareto principle, these would be the areas of Batatais SP, municipality requiring the implementation of mitigation measures in the first place.

**Figure 7.** Pareto diagram illustrating the frequency of land use conflict versus flood vulnerability classes in the municipality of Batatais SP.

The cross classification of land use and flood vulnerability is quantified in Table 7 and illustrated in Figure 8 in the form of a Pareto diagram. In this case, the Pareto principle applies because the 80% of cross-classified area is associated with 20% of land use classes, namely agriculture located in regions of low to high flood vulnerability and forest occupation in regions of medium vulnerability.


**Table 7.** Cross-classification of land use and flood vulnerability in the municipality of Batatais SP. The grey shaded cells represent concern situations because they correspond to medium and high flood vulnerability over large agricultural areas.

**Figure 8.** Pareto diagram illustrating the frequency of land use conflict versus flood vulnerability classes in the municipality of Batatais SP.

#### **4. Discussion**

The areas occupied by agriculture and simultaneously vulnerable to floods represent 516.80 km<sup>2</sup> of Batatais SP area, which means 60.73% of this municipality. From a decision-making stand point, these areas require intervention in the first place because they are agricultural areas located in regions of medium and high flood vulnerability. Although the scientific literature is scarce on actual or potential damage of floods in this region, the news periodically communicated by the social media (e.g., local television, newspapers) is enough reason to expose concerns about the potential harmful effects of recurrent floods for local agricultural systems and stream waters. Following the Pareto principle and common sense, these vulnerable areas require priority attention of public policies regarding decision making of preventive measures, public education and occupational rearrangement of the soil at catchment scale.

Floods submerge the fields used for agriculture, frequently causing production losses. The fields stay submersed for variable periods of time depending on factors such as topography, hydrologic barriers (e.g., dams), or return periods of the flood. The actual percentage of production loss will also depend on a number of factors, including crop variety, stage of plant development, flooded area, level of inundation and length of flooding period [32]. The present study does not provide numbers for all these factors but certainly indicate that large areas of Batatais SP municipality mostly used for sugar cane can be periodically submerged, especially where vulnerability is high. It is worth mentioning that the economic effects of river floods for sugar cane producers can be noteworthy. In a study in Pakistan [33], a reduction of sugar recovery from sugarcane was attributed to excessive water absorbed

by a flooded sugarcane crop which reduced the sugar contents in the cane and increased in the cane weight. The overall consequence was more payments to the local growers with less sugar recovery. In the context of a changing climate, this economic damage can aggravate if extreme precipitation events become longer and more intense [34]. Besides the direct tangible losses, the flood consequences for farmers also comprise direct intangible (e.g., business interruption), and indirect (clean up and recovery, soil and water remediation) losses [35].

Land-use changes contribute to an increased frequency and intensity of floods by increasing surface runoff [36]. Flood hazard events often occur in areas with dense population or high levels of agricultural land reclamation [37,38]. This gives urbanization and expansion of agriculture decisive roles in increasing the intensity and frequency of floods. A common cause of flood change response in modern agriculture areas relate to soil compaction from ploughing and heavy machinery [39–41]. However, flood response to land use changes in the rural environment usually begins during land conversion. Land use changes such as shifts from forestry to agriculture, from pasture to arable land, from rain fed to irrigated agriculture or from agricultural use to urbanized areas act as drivers of flood frequency and intensity increases [42,43], because these changes reduce land capacity to attenuate runoff. Some of these changes can produce amplified flood responses in case they correspond to environmental land use conflicts. In these cases, the change is performed in disagreement with land capability (natural use), for example when agriculture expands to steep forested hillsides. Flood intensity under these settings is amplified because land reduces the runoff attenuation capacity offered by forest cover while the large slope gradients accelerate runoff increasing flood potential downstream. Given the large area occupied of class 3 conflicts (Figure 6), characterized by forest-agriculture conversions, the Batatais SP municipality is presumably very prone to amplified flood intensities.

Land use conflicts may accelerate runoff and amplify flood response, but floods can feedback and amplify soil, water and biodiversity damages in areas where the conflicts have developed. Land use conflicts trigger changes in the soil composition and structure, namely reduction of organic matter and the capacity of soil to form stable aggregates [31]. Consequently, a cascade of environmental impacts develop, such as amplified soil loss [14,15], deterioration of water quality [44,45], and biodiversity decline [46]. Under these settings, recurrent floods will tend to aggravate these damages. If no prevention or recovery measures are taken, the flood-conflict issue may turn into a negative spiral ultimately leading to land degradation [47]. The return of a degraded land to a neutral condition is a priority given the economic and social costs of land degradation [48]. The neutral condition implies no net loss of the land-based natural capital relative to a reference state. Actions to achieve land degradation neutrality include sustainable land management practices that avoid or reduce degradation, coupled with efforts to reverse degradation through restoration or rehabilitation of degraded land. This rationale is independent from the causes of land degradation and hence applies to flood-related causes. However, to achieve land degradation neutrality is important to have a strong engagement of planners, stakeholders and the society, and may need to involve law enforcement initiatives [49].

The first step towards flood prevention is to map flood vulnerability [5,6]. In the present study, a map of flood vulnerability was prepared for application in a rural environment, using flood control attributes already adopted by other work groups [50], and a similar multi-criteria analysis [51]. However, to our knowledge this was the first time a flood vulnerability map has been coupled with a land use conflict map to investigate the potential impact of these inadequate land use changes on floods, and vice versa. Following the production of vulnerability maps the next step should be focused on the implementation of measures to actively control the floods, as well as measures to improve soil conservation and rural households' flood preparedness.

In general, active flood control can be accomplished through construction of detention basins (structural measures), and recent models attempt to optimize their location in a catchment [52]. However, besides the construction costs of these dams concerns exist about installing large detention

basins in catchments, because stream water stored in the dam reservoir can deteriorate extensively and rapidly [53]. To overcome this problematic situation, other recent models propose the construction of small, low-cost, and sustainable water detention basins, termed rainwater harvesting systems [54,55]. These models comprise tools to optimize the location of rain water harvesting systems in a catchment as to maximize their flood-control capacity and the complementary use for agro-forestry applications (irrigation, wildfire combat). In addition to flood control at the catchment scale, rain water harvesting systems can also be used to attenuate urban floods, which are important in various neighborhoods of Batatais town. In this case, water reservoirs are installed in rooftops to collect rainwater during storm events, thus reducing rainwater reaching the ground to become runoff. Furthermore, flood control in the urban environment the stored rainwater can be complementary used for sanitary or other purposes [56]. Therefore, the planning and installation of urban and rural rainwater harvesting systems is recommended as active flood-control measure in the studied area, especially in the context of a changing climate [57]. Other measures include non-structural flood management initiatives such as afforestation and restriction of urban development on floodplains. This set of measures has gained recognition in many countries as viable and cost-effective approaches to flood risk management [3,4,7].

For effective soil conservation, flood prone areas require the implementation of multiple complementary measures. Eventually, the most cost-effective measure is the reinforcement of manuring and composting of soil to raise the levels of organic matter [58] and produce stable aggregates that are resistant to detachment by rain drop action. Level plantation is another measure of soil erosion control, frequently used to reduce laminar erosion and increase soil water uptake. Level plantation can be successfully coupled with strip cropping [59] that involves the alternation of forages with strips of row crops, such as sugar cane. The control of soil erosion also comprises implementation of techniques that prevent soil compaction, such as no-tillage treatments [60], maintenance of crop residue to keep organic matter and nutrients in the topsoil, or more costly measures such as terraces in level or gradient since they reduce the ramp length reducing the surface drag of particles and nutrients [14].

The aforementioned measures of soil conservation act indirectly as stream water quality protective measures, because they tend to reduce sediment and nutrient loads transported in runoff. A direct protection of stream water quality from eroded soils and nutrients can be accomplished through creation of buffer strips along stream reaches [61]. Regardless of being a reliable measure to reduce deterioration of stream water quality, there is some debate around the buffer strip width capable to render these filters notable efficiency. For example, the new Forest Code in Brazil (Brazilian Law No. 12.651/12) has defined a minimum width of 30 m for efficient buffer strips, but some recent studies state that the reference should be raised to 50 m, at least [62].

All measures described above are likely to produce the proposed results (flood control and conservation of soil and water) if they can be aggregated as land consolidation–water management (LC-WM) plan capable to realize a multifunctional climate resilient rural area [63]. A decision support system will be required to adequately put this plan into practice [64,65]. Additionally, to be properly planned in the study area local watershed, committees, on which watershed management expectations stand, need to gain legitimacy to become mediators of sectorial and land use policies. To get such legitimacy, these committees need to acquire the status of State public entities offered by Brazilian Law No. 9433/97. Finally, to ensure effective implementation, a LC-WM plan needs to engage local communities, politicians and stakeholders in participatory design approaches [66] to evaluate the rural households' flood preparedness.

#### **5. Conclusions**

The Batatais SP municipality is a flood prone region characterized by environmental land use conflicts. Moderate floods can occur every 1–2 years, while severe floods have occurred every 10 years in this region. Natural flood vulnerability results from abundant precipitation over a territory with contrasting reliefs and altitude ranges. The most sensitive areas in this regard are the lower altitude lands that occupy 289 km<sup>2</sup> (34% of the municipal area). Environmental land use conflicts

developed because a vast portion of land became occupied by sugar cane plantations, even where land capability (the natural use) advised occupation by native forest vegetation. The conflicts potentiate amplification of flood intensity and inundation levels, because they reduce runoff retention capacity of soils. The study revealed 60% occupation of Batatais SP municipality with areas in environmental land use conflicts. This very large coverage is preoccupying. The study also exposed the situation and alerted for the need to engage state water planners, municipality authorities, and local communities in a process of land consolidation–water management involving a diversity of specific flood control, as well as soil and water conservation measures. The implementation of such plan would constitute a notable contribution to the sustainable use of soil and water resources in the region, with unquestionable benefits for the environment and Batatais society.

**Author Contributions:** A.M.C. contributed data curation, investigation, conceptualization, methodology, and writing original draft; T.C.T.P. and F.C.R.N. contributed to formal analysis, supervision, project administration and funding acquisition; M.Z. and R.d.B.V.P. contributed to resources and formal analysis; R.C.A.C. contributed to software and visualization; F.A.L.P. and L.F.S.F. contributed to formal analysis, validation, writing-review, editing, and funding acquisition.

**Funding:** As regards the Portuguese authors, the research was funded by national funds (FCT-Portuguese Foundation for Science and Technology) under the strategic project of the Vila Real Chemistry Research Center (PEst-OE/QUI/UI0616) and the CITAB (Centre for the Research and Technology of Agro-Environmental and Biological Sciences) project UID/AGR/04033.

**Acknowledgments:** The Brazilian authors wish to thank the São Paulo State University (Unesp)/School of Agricultural and Veterinarian Sciences, Campus Jaboticabal, São Paulo, Brazil; the CAPES—Coordenação de Aperfeiçoamento de Pessoal de Nível Superior; and the CNPq—Conselho Nacional de Desenvolvimento Científico e Tecnológico IF—Instituto Florestal, São Paulo, Brazil. The authors also wish to thank the University of Trás-Montes and Alto Douro (UTAD), the Chemistry Research Center and the Center for the Research and Technology of Agro-Enviromental and Biological Sciences (CITAB), for technical support.

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

#### **References**


© 2018 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
