**Effects of Land Use on the Ecohydrology of River Basin in Accordance with Climate Change**

Editors

**Wiktor Halecki Dawid Bedla Marek Ryczek Artur Radecki-Pawlik**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Wiktor Halecki Polish Academy of Sciences Poland

Dawid Bedla University of Agriculture in Krakow Poland

Marek Ryczek University of Agriculture in Krakow Poland

Artur Radecki-Pawlik Cracow University of Technology Poland

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Land* (ISSN 2073-445X) (available at: https://www.mdpi.com/journal/land/special issues/effects land use ecohydrology river basin accordance climate change).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-5645-1 (Hbk) ISBN 978-3-0365-5646-8 (PDF)**

Cover image courtesy of Dawid Bedla and Wiktor Halecki.

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


Reprinted from: *Land* **2022**, *11*, 535, doi:10.3390/land11040535 .................... **141**


## **About the Editors**

#### **Wiktor Halecki**

Wiktor Halecki (PhD) As a system data analyst, I specialize in computational ecology. Currently, I work as an assistant professor at the Institute of Nature Conservation Polish Academy of Sciences. Nature-based Solutions and blue–green infrastructure in urban areas are topics of my research. To conserve biodiversity, I am implementing ecosystem services, including the control of invasive species. Additionally, in my postdoctoral work, I contributed to an international project assessing how extreme weather events impact desertification using hydrological models. Among Wiktor Halecki's accomplishments is popularizing environmental protection knowledge. Over 370 articles of international scope were reviewed by me between 2017 and 2022. As a microblogger (@WiktorHalecki)/Twitter user, I present my research results in an accessible graphic form on various topics in biology and ecohydrology. A publication on forecasting climate change in drought-prone areas is in preparation.

#### **Dawid Bedla**

Dawid Bedla (PhD) In my research, I specialize in aquatic ecology, plant ecology and surface water pollution. I am currently working as an assistant professor at the University of Agriculture in Krakow. Water status assessment, blue–green infrastructure in urban areas and ecological indicators are the main topics of my research interests. I participated in an international project devoted to the implementation of bioeconomy, which was implemented from the budget of the "Horizon 2020" program supported by the European Union. I am the author of over 50 articles of international scope and a reviewer of 10 international publications in the years 2019–2022.

#### **Marek Ryczek**

Marek Ryczek (PhD) Marek Ryczek is an Associate Professor in Environmental Engineering, Mining and Power Engineering. He works at the University of Agriculture in Krakow (Faculty of Environmental Engineering and Surveying, Department of Land Reclamation and Environmental Development). His research is on soil water erosion, land reclamation, soil-science, and soil and water protection.

#### **Artur Radecki-Pawlik**

Artur Radecki-Pawlik (Prof.) Artur Radecki-Pawlik is a Professor of Civil Engineering (PhD, PE). He works at the Cracow University of Technology (Faculty of Civil Engineering). His research is on rivers and mountain streams, as well as low-head hydraulic structures, sediment transport, hydromorphology and fluvial geomorphology for engineers. He was previously employed at the Agricultural University of Krakow and by CBS and PBW HYDROPROJEKT Designing Office, Krakow.

## *Editorial* **Editorial: Effects of Land Use on the Ecohydrology of River Basins in Accordance with Climate Change**

**Wiktor Halecki 1,\*, Dawid Bedla 2, Marek Ryczek <sup>3</sup> and Artur Radecki-Pawlik <sup>4</sup>**

	- and Surveying, University of Agriculture in Krakow, 31-120 Krakow, Poland

A total of nine original publications and one concept paper are included in this Special Issue on water management and land use (Appendix A). Research topics include desertification in the Mexico Valley, the evaluation of environmental conditions and chemicals in Lithuania's water supply, and land use changes in the Polish mountains. This Special Issue also discusses climatic and anthropogenic factors that influence the Yangtze River in China. Similarly, the article on the Yellow River catchment fits quite well within Land's scope.

Climate change has changed the precipitation phase in Europe. More rain and less snow are likely to be experienced in winter. Snow replenishes groundwater resources that feed aquatic ecosystems (rivers, reservoirs, lakes) and water-dependent ecosystems (wetlands, marshes, and peatlands). Consequently, water shortages could arise and worsen as a result of a lack of snow or increased variability and intensity of precipitation. A variety of economic sectors and agriculture systems need to be properly managed.

Waterways with altered morphology must be revitalized or rehabilitated in order to restore their natural corridors and connectivity with river valleys. Groundwater retention capacity and restoration are also affected by such activities. Biodiversity enhancement, with an ecohydrological system, warrants climate change adaptation.

As an alternative to renaturation measures, the development of small, flowing water reservoirs could be a viable solution for floodplain and river banks morphology (straightened, staggered, and connected with groundwater). It is important to focus on high-quality, reusable retention water. Ecohydrological studies improve the quality of water and the environment by integrating biological indicators and hydrological process. Nature-based Solutions to water management are environmentally friendly.

#### **1. Water Management in Rural Areas**

Water retention has been greatly reduced in recent decades as the hydrological system has been remodelled to favour runoff rather than storage [1]. In agricultural areas, overdeveloped drainage channels, as well as flood protection measures, such as building dikes and straightening rivers in response to heavy rain, increase water loss. The sustainable development of catchments depends on preserving, protecting, or restoring natural watercourses [2].

Water management must integrate surface and groundwater resources, water quality, and reuse with environmental concerns [3]. Extreme weather phenomena limit the resilience of hydrological systems. Water deficits can be reduced more frequently, and the quality of the water might be improved simultaneously [4]. By integrating preventive steps, the hydrological buffer of the landscape would be increased, and more water could be retained. A rural landscape is designed to store water and maintain groundwater resources [5]. Generally, local water authorities and decision-makers have minor influence on social capital. Rural areas should, then, include the following [6,7]:

**Citation:** Halecki, W.; Bedla, D.; Ryczek, M.; Radecki-Pawlik, A. Editorial: Effects of Land Use on the Ecohydrology of River Basins in Accordance with Climate Change. *Land* **2022**, *11*, 1662. https:// doi.org/10.3390/land11101662

Received: 20 September 2022 Accepted: 22 September 2022 Published: 26 September 2022

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

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


#### **2. Blue–Green Infrastructure in an Urban Area**

Spatial planning should integrate green areas and urban structures to adapt to climate change [8]. Water retention measures should be created as an element of city planning activities (natural, semi-natural, and artificial areas). Blue–green infrastructure is undoubtedly an integral part of the built fabric (facilities combining greenery and water). In addition to significant negative environmental impacts (such as the desertification of water-dependent ecosystems, a decrease in groundwater-fed peatlands, the disappearance of alder and riparian forests), it has major economic and social consequences from an ecohydrological perspective [9,10]. The water quality of transformed rivers with a simplified biological structure and reduced ability to self-purify will be adversely affected by pollutants and elevated temperatures in the cities.

#### **3. Conclusions**

As a result of climate change, surface water quality may deteriorate and groundwater may become scarce. There are a number of meteorological factors, such as torrential rainfall, that can increase surface runoff, causing an accumulation of pollutants. An ecohydrologist explores both ecological processes and hydrological cycles. Studies and research in this field may be helpful in assessing environmental hazards. Climate change adaptation and mitigation can also be achieved through changes in land use, especially in river valleys. Small farms with hydrological-based production may also suffer from reduced productivity due to uncontrolled groundwater demand. Hence, groundwater exploitation requires monitoring and protection, which should set new standards for water resources. Nature-based Solutions can be the subject of further research with an emphasis on the eco-hydrological approach, especially in cities.

**Author Contributions:** Conceptualization, W.H. and D.B.; validation, M.R.; formal analysis, A.R.-P.; writing—original draft preparation, W.H. and D.B. All authors have read and agreed to the published version of the manuscript.

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

**Data Availability Statement:** Not applicable.

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

#### **Appendix A**


#### **References**


## *Article* **The Effects of Land Use on Concentrations of Nutrients and Selected Metals in Bottom Sediments and the Risk Assessment for Rivers of the Warta River Catchment, Poland**

**Michał Fiedler**

Department of Land Improvement, Environmental Development and Spatial Management, Faculty of Environmental Engineering and Mechanical Engineering, Poznan University of Life Sciences, 60-628 Pozna ´n, Poland; michal.fiedler@up.poznan.pl

**Abstract:** Changes in the environment, aiming at agricultural intensification, progressive urbanisation and other forms of anthropopression, may cause an increase in soil erosion and a resulting increase in the pollution inflow to surface water. At the same time, this results in increased nutrient pollution of bottom sediments. In this study, the concentrations of total nitrogen (TN), total phosphorus (TP), total organic carbon (TOC), calcium (Ca), iron (Fe) and potassium (K) were analysed using bottom sediment samples collected at 39 sites located along the entire length of the Warta River and its tributaries. Agricultural use of land adjacent to rivers was found to significantly degrade sediment quality, while anthropogenic land use (as defined by Corine Land Cover classification—CLC), unlike previous studies, reduces the pollution loads in the bottom sediments. Forest use also contributes to the reduction of the pollution load in sediments. It was found that the significance of the relationship between pollutant concentrations and land use depends on the length of the river–land interface. According to the analyses, the level of correlation between the analysed constituents depends on the use of land adjacent to rivers. The impact of agricultural land use has the strongest effect in the 1 km zone and 5 km in the case of anthropogenic land use. The results showed that the variability of total phosphorus TP concentrations is strongly correlated with the variability of iron concentrations. SPI values indicate that the risk to sediment quality is low due to TOC and Fe concentrations. In contrast, the risk of sediment pollution by TN and TP shows greater differentiation. Although the risk is negligible for 40% of the samples, at the same time, for 33% of the samples, a very high risk of pollution with both TN and TP was found.

**Keywords:** sediment; nutrient element; risk; land cover; Warta River

#### **1. Introduction**

River sediments are an important part of the cycling of materials in the aquatic ecosystem. The understanding of the pollution processes of river bottom sediments is of great importance because it can be an indicator of the ecological health of waters [1]. Sediments can accumulate pollutants and act as a buffer to absorb and release pollutants into the aquatic environment [2]. Loads of constituents such as carbon, nitrogen and phosphorus in bottom sediments and their interrelationships enable identification of the sources of organic matter and its transformations. The ratios of these constituents are affected by several environmental factors such as climate [3,4], terrestrial inflows, morphometry and use of adjacent areas [5] or the mineralogical composition of sediments [6]. Biochemical and biological transformations taking place in sediments and at the water–sediment interface [7,8], as well as the hydrodynamic processes of sediment transport in riverbeds [9], are also of great importance for concentrations and mutual proportions of constituents that can be found in bottom sediments.

The land-use pattern of the area located in the catchment is one of the more important factors affecting the concentrations of nutrients, as well as other pollutants, in river bottom

**Citation:** Fiedler, M. The Effects of Land Use on Concentrations of Nutrients and Selected Metals in Bottom Sediments and the Risk Assessment for Rivers of the Warta River Catchment, Poland. *Land* **2021**, *10*, 589. https://doi.org/10.3390/ land10060589

Academic Editors: Wiktor Halecki, Dawid Bedla, Marek Ryczek and Artur Radecki-Pawlik

Received: 16 May 2021 Accepted: 1 June 2021 Published: 2 June 2021

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

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

sediments [10]. Intense agricultural production can be a significant source of nutrient inputs to surface waters, both as diffuse and point sources of pollution [11,12]. This impact varies over time and depends, among other things, on changes in agricultural production technology [13]. Szatten and Habel [14] showed that the decrease of intensive agricultural areas results in the reduction of sediments and nutrients loads flowing into the reservoir. Additionally, urbanised areas were indicated as a source of surface water pollution [15,16] and, consequently, also bottom sediment pollution. However, as indicated by Khatri and Tyagi [15], the anthropogenic factors affecting water pollution in rural areas differ from those in urban areas. The chemical composition of bottom sediments of the water bodies located in urban and industrial areas shows anthropogenic enrichment with organic matter, phosphorus, calcium and trace elements [17]. An increase in impermeable surfaces in these areas increases surface runoff and erosion processes, causing an increased inflow of sediments to rivers [18]. In turn, afforested areas enable a reduction in the inflow of pollutants to watercourses and, consequently, reduce their deposition in bottom sediments. Riparian forest buffer systems are an effective tool for controlling the number of sediments and related pollutants carried in surface runoff [19]. Forest ecosystems also reduce the leaching of soluble nutrients into stream water [20]. In this regard, there is a frequent emphasis on the role of riparian vegetation that can act as a buffer to intercept suspended sediments and pollutants [21]. Forested riparian zones can reduce loads of suspended sediment in the runoff by 90% [22]. Grass riparian filter strips show a similar reduction of sediments in runoff as observed in forested zones, while the reduction of nitrogen and phosphorus is near 50% [23,24]. However, using only the primary type of land use may not be sufficient to explain the influx of pollutants to the aquatic environment fully. Liu et al. [6] showed that more fragmented forms of land use would give higher pollutant loads even if they are of the same primary type of land use.

The assessment of a qualitative status of bottom sediments requires the adoption of certain standards defining the concentrations of chemicals, the exceeding of which causes the sediments to be considered polluted. This is mainly due to the simplicity and ease of their application in practice. These standards, depending on a local controller, may address a variety of substances, including organic carbon, nutrients, heavy metals, organic compounds and others, and have different thresholds [25,26].

The chemicals that make up bottom sediments are interrelated. According to previous studies, the relationships between TN, TP and TOC concentrations are the result of biochemical and geochemical processes occurring in the environment and enable the identification of the sources of pollution [27,28]. Carbon–nitrogen ratio is one of the main variables to determine the source of organic matter in rivers [29]. Yun and An [30] indicated that the N:P ratio could be used in diagnosing the ecological health of a stream. Additionally, the relationships between TP concentrations and Fe, Ca and TOC concentrations are indicated, but the links between these elements can be different for each river [31].

The researches on the quality of riverine bottom sediments in lowland areas of Poland conducted in recent years has focused mainly on the analysis of concentrations and hazards caused by heavy metals [32–34], while the risks of nutrients or organic matter have not been analysed in detail.

This study aims to analyse the spatial variability of TN, TP, TOC, Ca, Fe and K concentrations and their interrelationships in bottom sediments of the Warta River and its tributaries. The analysis also considered the effects of land cover structure and P, Ca and Fe concentrations in soils, found in areas adjacent to sediment sampling sites, on sediment composition for different lengths of the contact zone between the river and surrounding areas. This enables the assessment of sediment quality and the identification of factors affecting the concentrations of nutrients in sediments.

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

#### *2.1. Site Description*

The catchment area of the Warta River is 54,519 km2. In terms of length, 808 km, it is the third-longest river in Poland (Figure 1). The tributaries of the Warta River, which were also analysed in the presented study, are the Prosna River; the Bobrowski Canal; the Mosi ´nski Canal; the Wełna, Obra and Note´c rivers as well as the additional tributaries of the Note´c River: the Gwda and Drawa rivers (Figure 1). The source of the Warta River is 380 m above sea level, and its mouth to the Oder is 12 m above sea level. The catchment area of the Warta River varies in terrain. The upper-southern part of the catchment is located in the Kraków-Cz ˛estochowa Upland, passing through the Central Polish lowlands to the northern lake districts. The tributaries of the Warta River range from 11 km (the Bobrowski Canal) to 391 km (the Note´c River). Their areas range from 28 km<sup>2</sup> in the case of the Bobrowski Canal to 17,319 km2 in the case of the Note´c River (Table 1).

**Figure 1.** Study site location.

Almost along the entire length of the Warta river as well as its tributaries drains a post-glacial landscape formed during the retreat of the Vistulian ice-sheet. The layer of glacial sediments with a thickness of more than 200 m consists of tills, glacial sands and gravels and outwash sands and gravels [35]. In the Warta River catchment area, sandy soils make up 45% of the land area, clay soils 41% and organic soils and alluvial soils 14% of the land area [36]. According to the European Soil Database v2.0 m, the main soils are Podzol—32%, Luvisols—30%, Fluvisols—12% and Arenosols—12%.

According to the CLC classification, in the whole catchment area of the Warta River, agricultural land use makes up 60.2% of the catchment area, forest and semi-natural land use make up 32.5%, anthropogenic land use makes up 5.6%, surface water makes up 1.4% and the rest is covered by wetlands and peat moors (Table 1, Figure 2).


**Table 1.** Analysed rivers' descriptions.

Crine Land Cover: 1—Artificial surfaces, 2—Agricultural areas, 3—Forest and semi-natural areas, 4—Wetlands, 5—Water bodies; \*—primary Warta tributary, \*\*—secondary Warta tributary (via Note´c).

**Figure 2.** Land cover of the River Warta basin according to Corine Land Cover.

The highest proportion of agricultural land has the catchment area of the Mosi ´nski Canal—approximately 78%. This catchment also has the lowest proportion of forest land and semi-natural land in a total catchment area—less than 16%. A slightly lower, compared to the aforementioned catchment, proportion of agricultural land in the total catchment area, 72%, is found in the catchment areas of the Prosna and Wełna rivers. On the other hand, the catchment area with the highest proportion of forest land and semi-natural forest land is the Drawa River catchment located in the northern part of the Warta catchment, where this proportion is 61.8%. This catchment also has the highest proportion of the open water area—3.9% of the total catchment area. The proportion of wetlands in all catchments is very small, reaching a maximum of 0.3% for the Mosi ´nski Canal catchment.

An important source of bottom sediment contamination is the material brought into rivers as a result of erosion processes. One of the indicators that identify areas sensitive to erosion is the Topographic Wetness Index TWI. The TWI takes into account the upslope area and its slope, allowing the calculation of the steady-state wetness and runoff across

the analysed area. Calculated TWIs show values in the range of 0.9–14.4 and are shown in Figure 3.

**Figure 3.** Distribution of Topographic Wetness Index TWI for the Warta River basin.

Geochemical conditions in the Warta River catchment, which describe the phosphorus, calcium and iron concentrations found in a 20 cm topsoil, are shown in Figure 4.

**Figure 4.** Classes of P, Ca and Fe content (GBC) in the upper layer of soils in the Warta River basin.

The Warta River catchment has the lowest annual precipitation in Poland [37]. Climatologically, according to the Köppen–Geiger's classification, the catchment area is located in the oceanic climate (Cfb) that has warm summers and milder winters [38,39]. The average annual air temperature for the WRC is approximately 8 ◦C, while the annual precipitation ranges from 520 mm in the north-eastern part of the catchment to 675 mm in the southern part [40]. The mean annual runoff is 3.86 dm3·s−1·km−<sup>2</sup> for the Warta catchment and it shows considerable spatial variability [41]. The lowest runoff values, below 2 dm3·s−1·km<sup>−</sup>2, are observed in the upper catchment of the Note´c River [42].

#### *2.2. Materials*

The analyses used data concerning the measurements of N, P, K, Ca, Fe and TOC concentrations in bottom sediments of the Warta River and its tributaries, which were obtained as part of the State Monitoring of the Natural Environment programme. Bottom sediment samples were collected from 39 sites in 2016 (Figure 1). Sampling sites were

located at the borders of the catchment area, at the mouths of tributaries and below cities where large industrial plants are located. For chemical analyses, a 5 cm top layer of sediments was collected from 4–5 sites over a 50 m distance for each location using a van Veen grab sampler. Samples for each site were mixed and rubbed through a nylon sieve with 2 mm mesh. Laboratory measurements were performed for Ca, K, Fe and P using inductively coupled plasma atomic emission spectrometry (ICP-AES) according to PN-EN 13657:2006, PN-EN ISO 11885:2009. Total organic carbon (TOC) was measured using coulometric titration according to CZ\_SOP\_D06\_07\_055 (CSN ISO 10694, CSN EN 13137). Kjeldahl nitrogen was determined using titration according to PN-EN 13342:2002.

The boundaries of the catchment area of the Warta River, its course and the course of its analysed tributaries were determined based on a Computer Map of the Hydrographic Division of Poland, which contains full hydrographic data of Poland in vector format. The land-use structure in the catchment area was made based on the CLC database provided by the Chief Inspectorate for Environmental Protection. The landform was described using a Digital Elevation Model (DEM) on a grid of at least 100 m, which was obtained from the Central Office of Geodesy and Cartography. The Ca, Fe and P concentrations found in the 20 cm topsoil were estimated using the Geochemical Atlas of Poland at a scale of 1:500,000, which was obtained from PGI—PIB (Polish Geological Institute—National Research Institute).

#### *2.3. Methods*

#### 2.3.1. Overall Statistics

To evaluate the variability of N, P, K, Ca, Fe and TOC concentrations in bottom sediments of the analysed rivers, the mean values of the concentrations, the median and their maximum and minimum values were calculated. Concentration distributions were determined using the Shapiro–Wilk test. Outliers were determined using the Interquartile Range (IQR) test. The values that fall below or above 1.5-IQR were considered as outliers of the variables:

$$\text{Llp} = \text{Q}\_3 + 1.5 \cdot \text{IQR} \tag{1}$$

$$Low = \text{Q}\_1 - 1.5 \cdot \text{IQR} \tag{2}$$

where *Up*—concentrations for values that fall above the 75th quartile (Q3) and *Low* concentrations for values that fall below the 25th quartile (Q1). Correlations were determined using Spearman's rank coefficients.

All statistical calculations and visualisation of results were performed using R version 4.0.5.

#### 2.3.2. Risk Assessment Method

An ecological risk assessment of nutrients in river bottom sediments was performed using a method based on the Single Pollution Index (SPI) [43].

The *SPIi* is calculated based on the formula:

$$SPI\_i = \mathbb{C}\_i / \mathbb{C}\_s \tag{3}$$

where *Ci* is a measured concentration of the evaluated *i* factor and *Cs* is a standard concentration of the evaluated *i* factor. Calculations were performed for TN, TP, TOC and Fe. Values of 550 mg-kg−1, 600 mg-kg−1, 1% and 2% were used as standard values for TN, TP, TOC and Fe, respectively, which are derived from safe nutrient concentration limits that can be found in the Sediment Quality Guidelines [25,43]. Based on the *Pi*, four hazard classes can be distinguished as shown in Table 2 [43].


**Table 2.** Single pollution index *SPIi* classification.

C:P:N stoichiometry enables the description of potential sources of pollution as well as geochemical and biochemical processes occurring in the environment [31,44]. TOC/TN is one of more frequently used ratios, which identifies potential sources of organic matter [43]. TOC/TP reflects to some extent the rate of conversion of organic carbon and phosphorus compounds [45]. In contrast, TN/TP reflects the dynamics of accumulation, deposition and release of nitrogen and phosphorus in water [46,47]. TOC/TN > 10 indicates primarily a terrestrial source of organic matter. In contrast, TOC/TN < 10 indicates an aquatic source of organic matter. TOC/TN~10 indicates that the organic matter found in sediments is of both terrestrial and aquatic origin [48,49].

2.3.3. Analyses of Spatial Variability of Concentrations of Selected Elements Found in Bottom Sediments

The principal component analysis (PCA) and cluster analysis (CA) were used for the differentiation of groups of sampling sites exhibiting similar variability and similar sources. Due to concentrations of analysed substances in bottom sediments, the cluster analysis of sampling sites was performed using k-medoids [50,51]. K-medoids is a technique that is less sensitive to the effect of noise and outliers in analysed data, compared to the k-means. The number of cluster groups was found using the silhouette method [52].

Land cover was defined for river-adjacent zones of a width of 500 m, in % of the total area of the zone. The adopted zone lengths were 1, 2, 5 and 10 km from the sampling sites (Figure 5). In zones defined as above, GBI's content index for P, Ca and Fe was also calculated based on the formula:

$$GBI = \Sigma(GBC\_i \cdot A\_i) / \Sigma A\_i \tag{4}$$

where *GBCi* means the topsoil content class for each element, while *Ai* is the area per class in an analysed zone.

**Figure 5.** Scheme of buffer zones for land-cover analysis.

Spatial variability data concerning the parameters were determined using QGIS 3.18 whereas statistical calculations were performed using R 4.0.

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

*3.1. General Characteristics of Concentration Variability*

The results of measurements made in 2016 for TN, TP, Ca, Fe and TOC concentrations found in the bottom sediments of the Warta River and its tributaries are shown in Table A1. Concentrations of individual elements are within N: 71.5–16,240 mg-kg−3, P: 48.8–8864 mg-kg−1, Ca: 423–54,010 mg-kg−1, Fe: 1602–68,280 mg-kg−<sup>1</sup> and TOC: 0.20–21.50 % d.m.

The average concentrations of the analysed elements ranged from 13,182 mg-kg−<sup>1</sup> for Fe through 9320 mg-kg−<sup>1</sup> for Ca, 2169 mg-kg−<sup>1</sup> for N to 1130 mg-kg−<sup>1</sup> for P and 3.76% for TOC (Table 3). The median values are 8269 mg-kg−1, 11,078 mg-kg−1, 998 mg-kg<sup>−</sup>1, 464 mg-kg−<sup>1</sup> and 2.18% for Fe, Ca, N, P and TOC, respectively. The TOC values measured in 1998–2000 are similar to those measured in bottom sediments of the Oder, which is a receiving body for the water of the Warta River, when they ranged from 0.2% d.m. to 18.2% d.m., averaging 4.9% d.m [53]. Similar values for P, Ca, Fe and TOC concentrations were observed by House and Denison [31,54] in seven rivers in southern England. Their reported data concerning concentrations were 40–47,000 mg-kg−<sup>1</sup> for Fe, 40–19,000 mg-kg−<sup>1</sup> for Ca, 10–2000 mg-kg−<sup>1</sup> for P and 0.6–19% d.m. for TOC. This may indicate similar sources of sediment pollution and similar land use of catchments.

**Table 3.** The descriptive statistics concentrations of chosen elements (mg·kg<sup>−</sup>1) and TOC (% d.m.).


The IQR values shown in Table 3 were 2685 mg-kg−1, 1084 mg-kg−1, 9934 mg-kg−1, 13,821 mg-kg−<sup>1</sup> and 3.18% for TN, TP, Ca, Fe and TOC, respectively. The IQR test showed several outliers involving N concentrations in the Warta River bottom sediments at site no. 4 and the Note´c River sediments at site no. 82. For phosphorus, concentration outliers were found in the Warta River sediments at sites no. 4, 19, 20 and in the Note´c River sediments, at site no. 82. For Ca, outliers were observed in the Mosi ´nski Canal at sites no. 51 and 53. For Fe, outliers were found in the Note´c River at site no. 82. Outlier TOC values were observed in the Warta River sediments at site no. 4 and in the Note´c River bottom sediments at sites no. 82, 83 and 84. The outliers were most frequently observed at site no. 82 (Note´c River)—four times and three times higher/lower at sites no. 4 and 23 (Warta River).

The analysis of distributions of the analysed concentrations showed that none of the concentration distributions is a normal distribution (Figure 6). In contrast, for all *p* > 0.05 parameters, all distributions are similar to a log-normal distribution.

All analysed concentrations are strongly positively correlated with one another. Spearman's correlation coefficients range from 0.50 for the correlation between TOC and Ca to 0.93 for the correlation between TP and Fe. Almost all correlation coefficients are significant at the *p* < 0.0001 level (Table 4). Only the TOC–Ca correlation is significant at the *p* < 0.01 level, whereas the TOC–K correlation is significant at the *p* < 0.001 level. The strongest correlation is between Fe concentrations and TP, Ca and K concentrations, for which the correlation coefficients are 0.93, 0.87 and 0.86, respectively. In contrast, the correlations between TOC and Ca/K are the least correlated, but still significant at the *p* < 0.01 level.

**Figure 6.** Variation of the concentrations of TN, TP, TOC, Ca, Fe and K in the bottom sediments of analysed rivers.

**Table 4.** Correlation coefficients between analysed elements.


Inorganic constituents that can combine with or adsorb phosphates are one of the factors affecting TP concentrations in bottom sediments [55,56]. Immobilising components include metals, such as iron, which combine with phosphorus, forming crystalline metal phosphates (e.g., strengite or vivianite), or they can adsorb phosphorus on the oxide/hydroxide layer. In the case of the presence of large amounts of calcium, phosphate precipitation may occur in the form of hydroxylapatite (HA) or various calcium phosphates [55,57,58]. Furthermore, iron (III) plays an important role in the phosphate sorption by humic elements, forming iron (III)–humus complexes with phosphates [55,59].

#### *3.2. C:N:P Stoichiometry*

The calculated characteristic values of C:N, C:P and N:P ratios indicate their high variability as shown in Table 5 and Figure 7. C:N values range from 4.51 to nearly 112, with a median value of 22.98. The highest C:N values of 112, 81, 68 and 61 were observed at sites no. 11, 52, 3 and 21, respectively. A total of 33 sites had C:N > 10, two sites had ratios similar to 10, and only four sites had C:N < 10 (Figure 7). This indicates that the organic matter found in bottom sediments derives mainly from runoff from surrounding areas. Similarly, high C:N values were also obtained by Lu [60], indicating that exogenous sources of pollution are the cause. On the other hand, the lowest C:N values have sediments at sites located in the vicinity of the aforementioned site no. 52, i.e., sites no. 51 and 53–4.5 mg-kg−<sup>1</sup> and 6.9 mg-kg<sup>−</sup>1, respectively. The bottom sediments at these two sites have high nitrogen concentrations as well as very high P, Ca and Fe concentrations (Table A1).


**Figure 7.** Spatial distribution of C:N, C:P and N:P ratios in the bottom sediments.

C:P ratio values in bottom sediments of rivers of the Warta River catchment ranged from 6.83 to 501.27, with a mean value of 58.3. The highest value of this parameter was obtained for site no. 83—the Note´c River, while the lowest for site no.—the Mosi ´nski Canal.

N:P ratios in the bottom sediments of rivers of the Warta River catchment ranged from 0.28 to 15.37, with a mean value of 2.32. According to Chen et al. [45], low N:P values correspond to eutrophic and mesotrophic conditions with the predominant supply of nutrients derived from external sources. The nutrients from such sources are complex and have low N:P ratios. In contrast, high TN/TP values can be observed under oligotrophic conditions when natural nutrient supply sources have high N:P ratios. Knösche [55] reports that the TN:TP ratio in river sediments decreases significantly from low-density Paleopotamal sediments (median TN:TP = 23.7) to high-density Eupotamal sediments (median 8.5). He also observed a similar correlation for the TOC:TP ratio. At the same time, Ostendorp [61] indicated that N:P:OC ratios could be explained by the variability of organic content rather than by eutrophication processes in each case.

#### *3.3. Single Pollution Index (SPI)*

The pollution status of bottom sediments of the rivers of the Warta River catchment shows great spatial variation, as shown in Figure 8. SPI values for nitrogen range from 0.13 to 29.5. The largest proportion of sediment samples, 41%, can be classified as Class I (Table 6). Class II includes 13% of the samples. Sampling sites with a low risk of nitrogen pollution are located mainly in the upper part of the Warta River (sites no. 1–14, excluding sites no. 4 and 12) and in its lower reaches (sites no. 18–24, excluding site no. 20). On the other hand, Class III includes 13% of sediment samples and Class IV, with the highest pollution, includes as many as 33% of samples. These sites are located primarily in the middle course of the Warta River and in its tributaries.

**Figure 8.** Spatial distribution of SPI for TN, TP, TOC and Fe.



SPI values for sediments polluted with phosphorus ranged from 0.08 to 14.8. Classification of sediment samples classifies 39% of samples into Class I, 15% into Class II, 10% into Class III and 36% into Class IV. These values are similar to SPI for TN. As with nitrogen, the sampling sites classified as Class I and II are located in the upper and lower reaches of the Warta River. The SPI values for TOC allow 23% of sediment samples to be classified as Class I and 54% as Class II. On the other hand, 15% and 8% of sediment samples should be classified into Class III and IV, respectively. This indicates a much lower risk of bottom sediment pollution by TOC than by TN or TP. The lowest pollution risk to bottom sediments was found for Fe. SPI values ranging from 0.08 to 3.4 enabled 56%, 23%, 18% and 3% of the samples to be classified as classes II–V, respectively.

The obtained SPI values indicate a lower pollution risk to bottom sediments of the Warta River than those obtained by He et al. [43] for Chinese rivers.

#### *3.4. Spatial Variation of Pollutant Concentrations*

The clustering analysis of TN, TP, TOC, Ca, Fe and K concentrations in bottom sediments using the Partitioning Around Medoids (PAM) identified three groups of bottom sediment sampling sites (Figure 9). A simultaneous principal component analysis (PCA) found that the first two principal components, for which the eigenvalues are greater than 1, explain approximately 87% of the total variance (Table 7).

**Figure 9.** Partitioning clustering of sampling sites according to element concentrations in sediments (**A**) and location of points (**B**).


**Table 7.** Eigenvalues, percent of variance and cumulative percent of variance for the PCA.

The land-use structure, as determined by the CLC classification, shows varying effects on the levels of constituents found in river sediments of the Warta River catchment. Anthropogenic lands (C\_1) and forests and semi-natural ecosystems (C\_3) show a negative correlation, while agricultural lands (C\_2) and lands covered by water (C\_5) show a positive correlation with concentrations of all elements; however, this relationship is not statistically significant in most cases. It should be noted that previous studies highlighted the negative impact of urbanised areas on the quality of bottom sediments in rivers. This was mainly due to the impact of direct discharge of wastewater into rivers or increasing surface runoff from areas with impervious surfaces [62,63]. However, the following phenomena that have been marked in recent years, such as the connection of most wastewater sources to wastewater treatment plants, the increase of areas that enable rainwater infiltration or, for example, measures to retain rainwater runoff from impervious surfaces, may reduce or neutralise the negative impact of urbanisation on the quality of surface water and bottom

sediments. The negative impact of agricultural areas is mainly caused by an increased input of nutrients [12], inorganic suspended solids [64] and organic matter [65] from these areas to surface waters.

Factor loadings for the first two principal components are shown in Table 8. According to the classification presented by Liu et al. [66], factor loadings were classified as "strong," "average," and "weak" for values of >0.75, 0.75–0.50, and 0.50–0.30, respectively. PC1, which explains 67.4% of the total variance, has strong positive Fe, TN, TP and TOC loadings, while also showing average K and Ca loadings.

**Table 8.** Loading for PCA matrix.


\*\*—strong effect; \*—moderate effect.

Table 9 shows the Spearman's correlation coefficient values of TN/TP/K/Ca/Fe/TOC concentrations in bottom sediments with land cover structure in the zones of a width of 1000 m, adjacent to the river on a length of 1, 2, 5 and 10 km, and the calculated Ca, Fe and P content index of the topsoil.

The significance of the land-use impact depends not only on the land-cover structure itself, but also on the length of the contact zone of land adjacent to the watercourse. A significant (*p* < 0.05) effect of agricultural land (C\_2) on the concentration levels of all analysed elements was observed for the 1 km buffer area. For longer buffer lengths, excluding K and TOC concentrations in the 2 km buffer area, that effect was no longer statistically significant. In contrast, the proportion of anthropogenic land (C\_1) is significantly linked to P, K, Fe and TOC concentrations in the 5 km buffer area. Additionally, the proportion of the area covered by water (C\_5) shows the strongest relationships with P, K, Ca and TOC concentrations in the 5 km buffer area. At the same time, the proportion of the land covered by water is the only one that is significantly correlated with phosphorus concentrations for all buffer lengths. Except for the effects of forest land proportion (C\_3) on Fe concentrations in the 10 km buffer area, there was no significant impact of that land-use category on concentrations of analysed substances found in bottom sediments.

The results indicate a significant positive relationship between the geochemical background level of phosphorus in the 2 km buffer area and N, P, Ca, Fe and TOC concentrations in the bottom sediments analysed. In contrast, the geochemical background levels of Ca and Fe do not show such a relationship.

#### *3.5. The Relationships between Water Quality Parameters and Land Use*

The analysis of land use as identified based on the CLC data in individual buffers enabled four land-cover groups to be distinguished. The clustering results and the distribution of each group are shown in Figure 10. The first group includes areas dominated by anthropogenic lands, with a low proportion of forests. The second group includes buffer areas that are predominantly forested. The third group includes buffer areas dominated by agricultural lands. The fourth group consists of buffer areas where the proportion of lands covered by water is at least 10%, regardless of other land-use categories.


**Table 9.** Correlation coefficient matrix showing element–land cover and element–geochemical background relationships in sediments.

\*—significant at *p* < 0.05; \*\*—significant at *p* < 0.01; Land Cover: C1—Artificial, C2—Agricultural, C3—Forests and seminatural, C4—Wetlands, C5—Water.

**Figure 10.** Classification of sampling points according to land use.

The Spearman's correlation coefficients between parameters describing constituent concentrations in bottom sediments for defined land-use groups are shown in Table A2. In group 1, where anthropogenic land use is the main factor, the strongest relationship is between TP and Fe, TN and TOC concentrations. However, the relationships between the elements show a higher correlation in the case of buffers that were identified based on a significant proportion of forest land. The highest correlation is between TP and TN, TN and TOC and TP and Fe. Unlike cluster 1, there is also a strong correlation between TOC and Fe in cluster 2.

The highest levels of correlation between bottom sediment constituents occur when the river is adjacent to agricultural land (cluster 3). In this case, all correlation coefficients are statistically significant, at least at the *p* < 0.05 level. The highest correlation in cluster 3 is between Ca and K, TP and Fe and Fe and K. The lands included in cluster 4 have surface water within the buffer area, the proportion of which is between 10% and 20% (highest proportion) of the total buffer area. In the case of sediment samples collected at sites assigned to this cluster, the highest correlation occurred between Ca and K concentrations, where r = 1.0. Again, very high levels of correlation occurred between Fe and TP and Ca and K.

It should be noted that regardless of the influence of the environment, phosphorus concentrations show the highest correlation with iron concentrations and slightly lower correlation with calcium concentrations, which is due to the already mentioned mutual interactions between these constituents, leading to precipitation or release of phosphorus [67–69]. The strong correlations between TP and TN concentrations for all use categories of the river-adjacent zones, except for the relatively low correlation for agricultural use, may indicate similar sources of nutrients in these areas. In agricultural lands, however, sources of pollution are of different nature—diffuse or point. Nutrients may derive from both field crops and animal husbandry, and they have different nutrient ratios, as also indicated by Lee [70] and Jones [10]. At the same time, bottom sediments collected from sites included in cluster 3 have the highest average TN, TOC, Ca and Fe concentrations and also high TP concentrations (Table 10). In contrast, the concentration values for clusters 1 and 2 (i.e., river-adjacent zones with anthropogenic and forest land uses) are significantly lower. Average sediment constituent concentrations for cluster 3 compared to concentrations for clusters 1 and 2 are 4.3–5.8 times greater for nitrogen, 2.3–4.0 times greater for phosphorus, 3.4–4.2 times greater for TOC, 4.3–5.8 times greater for calcium and 1.9–2.8 times greater for iron. The intermediate values of constituent concentrations were found in the sediments for cluster 3 determined on the basis of the proportion of surface water in the river-adjacent zone. As this proportion is between 10–20%, the remaining area of the land may have very different uses and, consequently, average out the effects of different land uses on the bottom sediment composition.


**Table 10.** Mean concentrations of elements (mg·kg<sup>−</sup>1) and TOC (%) in sediments samples for landuse clusters.

#### **4. Conclusions**

Based on the results of analyses concerning spatial variability of N, P, TOC, Ca, Fe and K concentrations in the bottom sediments of the Wart River and its tributaries, the following conclusions can be made:

• The dominance of agricultural activities in the zone of land adjacent to the river was found to be a significant factor increasing the pollutant content found in the bottom sediments of the analysed rivers. On the other hand, artificial and wooded areas reduce the concentrations of analysed substances in bottom sediments. The significance of this impact also depends on the length of the contact zone between the river and the land adjacent to it.


**Funding:** The publication was co-financed within the framework of the Ministry of Science and Higher Education programme as "Regional Initiative Excellence" in years 2019–2022, Project No. 005/RID/2018/19.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Appendix A**

**Table A1.** Concentrations of chosen elements (mg·kg<sup>−</sup>1) and TOC (% d.m.) in the bottom sediments of the Warta river.



**Table A1.** *Cont.*

<sup>1</sup> Kanał Bobrowski, \* outliers.

**Table A2.** The correlation coefficients between the concentrations of elements in the bottom sediments for buffer zones with various land covers.


\*—*p* < 0.05, \*\*—*p* < 0.01, \*\*\*—*p* < 0.001; *n*—number of samples.

#### **References**

1. Liao, J.; Chen, J.; Ru, X.; Chen, J.; Wu, H.; Wei, C. Heavy metals in river surface sediments affected with multiple pollution sources, South China: Distribution, enrichment and source apportionment. *J. Geochem. Explor.* **2017**, *176*, 9–19. [CrossRef]

2. Petticrew, E.L.; McConnachie, J.L. Sediment-water interactions. In *Sediment Dynamics and Pollutant Mobility in Rivers*; Westrich, B., Förstner, U., Eds.; Environmental Science and Engineering; Springer: Berlin/Heidelberg, Germany, 2007; pp. 217–268. ISBN 978-3-540-34782-8.


## *Article* **Hydromorphological Inventory and Evaluation of the Upland Stream: Case Study of a Small Ungauged Catchment in Western Carpathians, Poland**

**Łukasz Borek \* and Tomasz Kowalik**

Department of Land Reclamation and Environmental Development, University of Agriculture in Kraków, al. Mickiewicza 24–28, 30-059 Krakow, Poland; tomasz.kowalik@urk.edu.pl

**\*** Correspondence: lukasz.borek@urk.edu.pl

**Abstract:** The hydromorphological conditions of watercourses depend on numerous natural and anthropogenic factors such as buffer zones or human infrastructure near their banks. We hypothesised that, even in a small stream, there can be substantial differences in the hydromorphological forms associated with naturalness and human impact. The paper aims at the field inventory and evaluation of the hydromorphological conditions of a small upland stream in the conditions of contemporary human activity, against the background of meteorological and hydrological conditions. The study concerned a left-bank tributary of the Stradomka River located in the Wi´snicz Foothills (Western Carpathians). The analyses were conducted with the use of the Polish method, the Hydromorphological Index for Rivers (HIR), which conforms to the EU Water Framework Directive (WFD). The hydromorphological condition and quality of habitats were evaluated based on the Hydromorphological Diversity Score (HDS) and Habitat Modification Score (HMS). The study shows that the largest changes in stream hydomorphology and habitat conditions took place in the downstream, urbanised stream catchment area with an intensive development of construction and technical infrastructure. The hydromorphological condition of the examined stream sections was evaluated as good or poor. The best hydromorphological conditions were found in the section located in the semi-natural area, and the worst in the urbanised area. As our research shows, the strong influence of human activity, including weather extremes, and the risks and hydrological hazards of the hydromorphological conditions of the small, ungauged catchment, highlight the necessity to search for other research methods to support the decision-making cycle in the transformation of riverbeds and catchments.

**Keywords:** hydromorphological diversity; highland watercourse; human activity; catchment management; weather extremes

#### **1. Introduction**

Upland and mountainous rivers, streams and brooks have large flow variations during a year, which are caused by higher precipitation than in the lowlands, as well as the varied topography—especially the gradient (slope) of the river which, for upland and mountain streams, generally ranges from 2% to 7% [1]. The density of the river network is greater than on the lowlands; hence, the dynamics related to water flow and changes in the river catchment are significant. The evaluation of the ecology conditions/potential of a WFD body of surface water, to date, has mainly focused on larger watercourses (classified as rivers) with a catchment area greater than 10 km2, omitting their numerous tributaries (streams, brooks); the latter permanently or temporarily supply them with water, and are considered—often incorrectly—as artificial watercourses (e.g., ditches) formed as a result of land reclamation. Some of these small watercourses, often called 'wild streams', have never been examined in terms of their hydromorphology. Some of them are unnamed, not supervised and are not classified in terms of their natural or artificial origin. The beds of numerous, but mainly small mountain and upland streams are often dry, and

**Citation:** Borek, Ł.; Kowalik, T. Hydromorphological Inventory and Evaluation of the Upland Stream: Case Study of a Small Ungauged Catchment in Western Carpathians, Poland. *Land* **2022**, *11*, 141. https:// doi.org/10.3390/land11010141

Academic Editor: Ilan Stavi

Received: 21 December 2021 Accepted: 14 January 2022 Published: 17 January 2022

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

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

are classified as 'temporary sporadic', while heavy rainfall causes streams and rivers to swell, sometimes leading to landslides. The Carpathian streams also experience high water erosion which carves the flysch bedrock [2]. An important part of understanding a stream's hydromorphology is to undestand the catchment's geological structure and topography [3,4]. The Carpathian region (including the Wi´snicz Foothills) is built by flysch formations, Cretaceous and Eocene shales, marls and sandstones. They are covered with several inches of loess layers from Quaternary sediments [5]. The stream catchment area can be divided into three sections: upper, middle and lower [6]; these three sections are clearly distinguishable in mountain streams, and are hard to distinguish in upland streams. The heads of mountain and upland streams are usually located on mountain sides covered with trees or other vegetation.

The role of such watercourses in mountainous and upland areas is the natural drainage of the catchment areas [7]. Far-reaching anthropogenic changes in recent years, such as the development/sealing of river catchment areas, amongst other changes, cause various problems. Often, such small watercourses/streams/brooks are significantly changed or eliminated, e.g., for land development. The anthropogenic modifications involving the transformation of open riverbeds, streams and brooks into concrete pipelines can increase the occurrence of inundations and so-called flash floods [8–10]. Studies show that changes in climate and the use of land and water must be considered together to fully understand watershed hydromorphology [11]. A great deal of attention is given to the pollution of surface waters, and not enough attention has been paid to the effects of transformations in the beds and micro-catchments of small upland and mountain streams that form an integral part of larger river catchment areas.

The term 'hydromorphology' relates to the hydrology and geomorphology of a small watercourse, and includes the analysis of the physical attributes of the watercourse bed, such as flow type, material of the bottom and bank slopes, channel modification, and its natural elements; it also includes the characteristics of the land-use within 5 m and 50 m of the bank-top [12,13]. Hydromorphology is one of the key elements of the ecological integrity of waters, according to the EU Water Framework Directive WFD [14]. Member States shall, in accordance with this directive and European Standards EN 14614 [15] and EN 15843 [16], standardise provisions regarding the hydromorphological evaluation of a watercourse, allowing a better understanding of its functioning, and, if necessary, correct regulation of the riverbed [6]. The basic term used in hydromorphology is catchment, i.e., the area limited by the watershed in which the water flows to one place (the main watercourse/recipient). This is a dynamic system that depends on natural factors. The catchment morphodynamics vary and depend on elements of the natural environment, such as geological structure, climate and weather conditions, landform, river gradient/slope, soil and flora [17]. Under unfavourable conditions (e.g., flash floods), a seemingly small stream can cause damage to both a catchment area and human property.

Hydromorphology still occupies a small space in the water environment management of Europe and Poland. With the constant pressure of climatic change and changes in politics, society and economy, the restoration of natural habitats in our water environment is often regarded as a non-priority, despite the fact that hydromorphology supports the biodiversity of our waters. Good hydromorphological condition is a basic element of ecosystem health and a foundation of many ecosystem services and benefits for society [18]. The hydromorphological evaluation of a watercourse can be used for planning the restoration of rivers and streams by means of their regulation or biotechnical development [6]. There is evidence that even a small hydromorphological interaction can have a deep impact on the functioning of an ecosystem [19]. There are numerous methods for the hydromorphological evaluation of watercourses [20–22]. The most commonly used methods for assessing the quality of physical river habitats are: the German LAWA-FS [23]; the British River Habitat Survey, RHS [6]; the Spanish assessment of bank habitats, QBR; and the Czech comprehensive morphological assessment, HEM [24]. These methods are based on field surveys and the characterization of physical attributes of riverbeds and flow regimes, and can be

classified as river habitat surveys or physical habitat assessments [20]; despite the different assumptions, the methods lead to similar results and can be used in various countries, especially in Europe [24]. From the groups of methods used for field inventory of rivers for their restoration purposes, those recommended are: the Australian River Styles Framework [25]; the IHG method in Spain [26]; the MQI method in Italy [27]; and the Polish method of river hydromorphological quality assessment RHQ [28]. The index-based methods are being developed to allow the standardisation of methods and tools used for evaluation of hydromorphology and impact on the the ecological condition of watercourses/catchment areas [29]. In Poland, the Hydromorphological Index for Rivers (HIR) method, which conforms to the EU Water Framework Directive (WFD), has been used since 2017. The method provides the assessment of lowland, upland and mountain rivers and streams, and it can be used to evaluate natural and heavily modified watercourses, as well as artificial channels [30]. The HIR method is based on the British River Habitat Survey method, which is widely used for the hydromorphological assessment of waters around the world, as demonstrated by a number of citations in other research papers more than 2,000,000 times (on the Google Scholar website). The majority of the research on the hydromorphological assessment of watercourses has focused on large river systems. This work concentrates on the hydromorphological evaluation of a small upland stream, whose catchment area is being transformed as a result of high human activity. It is considered that, for the small ungauged aquatic ecosystem, it is advisable to know the hydromorphological conditions, which affect the flow of water and biodiversity.

We hypothesised that, even in a small stream, there is a significant variety of hydromorphological forms associated with naturalness and anthropopressure. The paper aims to provide a field inventory and evaluation of the hydromorphological conditions of a small ungauged stream catchment using the HIR method, including the meteorological and hydrological conditions.

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

The stream catchment area is located in southern Poland (Figure 1), in the Wi´snicz Foothills, part of the Western Carpathians. Empirical studies were conducted in the bed and valley of the 'unnamed stream', a left-bank tributary of the Stradomka river of the left-bank tributary of the Raba river (in the Upper Vistula Basin). The unnamed stream is a natural, periodic watercourse. In the currently valid typology of surface waters [31], it is classified as an upland carbonate stream with fine-grained substrate on loess and loess-like rocks (type code: 6). The studies were divided into 3 stages:


In terms of climate, according to the Köppen-Geiger classification [32], the examined catchment is located in a warm, humid continental climate (Dfb). The datasets from 2 meteorological stations from a multi-year period (2001–2020) were used in the analysis. The monthly precipitation totals (from the Łapanów station; GPS coordinates: 49◦51 44" N, 20◦16 32" E) and monthly average temperatures (from the Łazy station; GPS coordinates: 49◦57 54" N, 20◦29 43" E), provided by the monitoring network system and carried out by the Institute of Meteorology and Water Management–National Research Institute (IMGW-BIP), was used in the study. Both stations are located in the basin of the Upper Vistula within a distance of 20 km. The precipitation characteristics for monthly and annual periods were made on the basis of the relative precipitation index (RPI). The RPI classifies precipitation in terms of its excess or shortage [33]. The precipitation and thermal characteristics of the year 2020 were presented against the background of the long-term period 2001–2020. The RPI coefficient was calculated using the following formula:

$$RPI = \frac{P}{\overline{P}} \cdot 100 \text{ (\%)}\tag{1}$$

where *P*—precipitation sum in the studied period (mm), and *P*—average precipitation value in the studied long-term period (mm).

**Figure 1.** Location of the studies and sections with control profiles on the stream.

Based on the value of the RPI, the months and years were classified as follows: *RPI* < 25—extremely dry; 25 ≤ *RPI* < 50—very dry; 50 ≤ *RPI* < 75—dry; 75 ≤ *RPI* ≤ 125 normal; 125 < *RPI* ≤ 150—wet; 150 < *RPI* ≤ 200—very wet; *RPI* > 300.0—extremely wet [34].

The thermal classification of months, seasons or years was prepared according to Lorenc [35,36] in relation to the average monthly temperature (T) and the average monthly temperature calculated over the multi-year period (Tav), increased or reduced by the standard deviation value (δ) (Table 1).

**Table 1.** Thermal classification of months, seasons and years.


Explanations: T—air temperature in a given period; Tav—long-term average temperature in a given multi-year period; δ—standard deviation of air temperature in a given period. Source: elaborated according to Lorenc (2000).

A direct measurement of the water flow rate (*Q*) was performed using a bucket of known volume (*V*) and a stopwatch (*t*). The *Q* was calculated using the following formula:

$$Q = \frac{V}{t} \tag{2}$$

where: *<sup>Q</sup>*—water flow rate (m3·s−<sup>1</sup> or dm3·s<sup>−</sup>1); *<sup>V</sup>*—volume flow rate of water (m<sup>3</sup> or dm3); and *t*—time (s).

The average annual water flow (*SSQ*) for small ungauged catchments located in the Carpathian region was calculated by Punzet formula [37–39]:

$$SSQ = 10^{-3} \cdot SSq \cdot A \tag{3}$$

$$SSq = 0.00001151 \cdot P^{2.05576} \cdot I^{0.0647} \cdot N^{-0.04435} \tag{4}$$

where: *SSQ*—average annual water flow (m3·s−1); *SSq*—average annual surface runoff (dm3·s−1·km<sup>−</sup>2); *<sup>A</sup>*—catchment area (km2); *<sup>P</sup>*—average annual precipitation (mm); *<sup>I</sup>*—river slope indicator (‰); and *N*—soil imperviousness index (%).

The soils in the catchment area were comprised of regoliths and loess-like dusty formations formed as a result of flysch weathering and simultaneous eolian sedimentation.

The physiographic characteristics of the catchment area included the 15 parameters [40,41] listed in Table 2, which presents: catchment geometry and shape, catchment morphometry, and hydrological conditions. On a 1:10,000 vector topographic map were areas with the following forms of land use: arable land, grassland and wasteland, forests (tree-covered areas) and built-up areas. The share of individual land use forms was calculated based on the total areas.


**Table 2.** Characteristics describing the stream catchment physiographic conditions.


**Table 2.** *Cont.*

Explanations: \* [–]—dimensionless indicator/parameter.

The stream, at a length of about 250 m, carries water in the pipeline before spilling into the Stradomka river, a consequence of the progressive residential and service developments from the mid-1980s. In the middle part of the stream, a car park was built, which was also associated with the transformation of the water flow from an open stream bed into a closed pipeline. The increasing encroachment of urban areas in the stream catchment will follow the adopted land-use plan.

The hydromorphological evaluation of the 'unnamed stream' was performed in August 2020 with the use of the Hydromorphological Index for Rivers (HIR), which comprised of a field evaluation on two representative 300 m sections—semi-natural and urbanised in which 10 control profiles (marked as P0–P10) were delimited about 30 m apart from each other, as well as a synthetic evaluation of the entire stream section (Figure 1). The study involved the determination of the presence and diversity of the natural elements of the stream and of the valley, and characterisation of the range of modifications in the watercourse's morphology [30,42].

The land use in the stream valley was determined based on the orthophotomap from 2020, available at geoportal.gov.pl (accessed on 20 December 2021) (supported by a site visit), within a 100 m wide buffer on each bank; as a percentage of the area and then used to determine the weight coefficient (w) to calculate the mean HIR. Out of three basic land use forms, i.e., urbanised (U), agricultural (R) and semi-natural (S), the urbanised areas dominated in the buffer—sealed surfaces with dense or dispersed development, roads and other anthropogenic non-built-up areas, which made up about 58% of the buffer area (wU = 69.9). Next, there were semi-natural areas—forests, trees and shrubs, marshy areas, rush and herb vegetation, covering about 25% of the buffer area (wS = 30.1). The remaining 17% were agricultural areas—arable land, grassland and allotment gardens (wA = 0.0). Only the land use forms with a share ≥25% in the buffer were used in the calculation [42]. Two study sections were determined based on the land use analysis in the stream valley: the first in the semi-natural part, and the second in the urbanised part. Due to the stream features (its small length and partial pipelining) and the property conditions in the catchment area (the stream flows through private properties, often fenced), the study sections lengths were

300 m long instead of the standard 500 m, which is acceptable according to the authors of the method [30,42].

For each stretch, the HIR value was calculated. The multimetric HIR index combined two indices: The Hydromorphological Diversity Score (HDS) and the Habitat Modification Score (HMS). The HDS informed about the presence of natural attributes in the channel, coastal zone and the river valley. Each of the HDS attributes delivered a range of points, enabling the calculation of the HDS of the river stretch. HMS provided information on the hydromorphological modifications; it included various forms of fluvial ecosystem transformations, such as profile modifications and reinforcements, and the presence and abundance of engineering facilities. The detailed procedure of scoring is presented in the HIR method manual. In practice, HDS and HMS values usually do not exceed 100 [30,42].

$$HIR = \frac{\left(\frac{HDS - HMS}{100}\right) + 0.85}{1.8} \tag{5}$$

The obtained values were compared with the classes for a given type of watercourse: class I *HIR* ≥ 0.824; class II *HIR* ≥ 0.715; class III *HIR* ≥ 0.600; class IV *HIR* ≥ 0.485; and class V *HIR* < 0.485.

The assessment of the hydromorphological state of the stream catchment consisted of calculating the weighted average of the HIR values, taking into account the weighting factor for three forms of land development (calculated on the basis of their percentage share in the buffer), according to the following formula:

$$HR\_{mean} = \frac{(HIR\_U \cdot wLI) + (HIR\_A \cdot wA) + (HIR\_S \cdot wS)}{100} \tag{6}$$

where: *HIRmean*—average HIR value for the whole water body; *HIRU*—HIR calculated for the survey sites located in an urbanised area; *wU*—weighting factor for urbanised areas calculated on the basis of their percentage share in the buffer; *HIRA*—HIR calculated for the survey sites located in an agricultural area; *wA*—weighting factor for agricultural areas calculated on the basis of their percentage share in the buffer; *HIRS*—HIR calculated for the survey sites located in a semi-natural area; and *wS*—weighting factor for semi-natural areas calculated on the basis of their percentage share in the buffer.

#### **3. Results**

#### *3.1. Meteorological and Hydrological Characteristics*

The characteristics of the meteorological conditions are presented in Figures 2 and 3. The climate is fairly dry in the winter half-year, and is characterised by fairly hot summers, where short rain showers are quite common and often come in intense bursts.

**Figure 2.** The characteristics of the precipitation conditions: (**a**) average annual precipitation in 2001–2020; (**b**) average monthly precipitation from the multi-year period 2001–2020 and in 2020; (**c**) classification of precipitation conditions on the basis of Relative Precipitation Index (RPI) for the months and years.

Typically, 30% of the precipitation falls during the winter half-year and 70% during the summer half-year. In the multi annual period from 2001–2020, the average annual air temperature was 9.1 ◦C, and total precipitation was 792 mm. Increased average temperatures and rainfall are observed over the analysed period, which is in line with a report by the Institute of Meteorology and Water Management—National Research Institute, which has been monitoring Poland's climate for over 100 years on an ongoing basis [43]. To assess the meteorological conditions of the region where the study was conducted, the relative precipitation index (RPI) values were calculated, as a measure of the precipitation efficiency in a given month and year. Based on the RPI values calculated for the months from 2001–2020, as shown in Figure 2c, it was determined that 112 of the months were dry, 63 of the months were optimal and 65 of the months were wet. In terms of precipitation, 2020 belonged to the wet year. The average multi annual temperature varies between 7.8 ◦C and 10.3 ◦C (Figure 3c). The water temperature depends on the air temperature which, in turn, translates into the amount of oxygen in the water [44]. The annual average temperature in 2020 was 0.7 ◦C higher than the multi annual average. In the analysed multi-year period, two years were very cold, six years were slightly cold, six years were normal, two were slightly warm, three years were warm and one year was very warm. In general, the variable meteorological conditions in the study area, characterised by periods of excesses and shortages of water, indicate the need for good maintenance of the hydromorphological conditions of watercourses, mainly drainage functions.

**Figure 3.** The characteristics of the thermal conditions: (**a**) average annual temperature in 2001–2020; (**b**) average monthly temperature from the multi-year period 2001–2020 and in 2020; (**c**) monthly thermal classification proposed by Lorenc (2000).

The average annual water flow of the unnamed stream is placed at *SSQ* = 0.004 m3·s−<sup>1</sup> and corresponds to the flow measured directly (Q = 0.005 m3·s−1). The catchment basin of the stream receives water both from rainfall and from snow melt in the slope. The higher precipitation levels in 2020 lead to a higher supply of water to the stream. Based on conversations with residents, there have been dry years or months in which the stream did not carry water.

#### *3.2. Catchment Characteristics*

The basic physiographic parameters describing the catchment of the examined stream are presented in Table 3. The stream catchment area is very small (<10 km2). The catchment length is just over 1 km, its mean width is 0.23 km, and the perimeter is about 3 km. The form coefficient (*Cf*), elongation coefficient (*Cw*), circularity coefficient (*Ck*) and compactness coefficient (Cz) indicate that the stream catchment is narrow and elongated. The denivelation of the catchment area is 59 m, with a watershed slope of about 23.9‰. The mean altitude in the catchment is about 258 m a.s.l.


**Table 3.** Physiographic parameters of the studies stream catchment.

The dominant land use is arable land with slightly over 40% of the catchment area. The catchment afforestation rate is 6%. The trees and shrubs occur mainly in the stream valley, thus forming the natural biological development of the watershed bed. The urbanised area (dispersed and dense development) is 25%. Grassland (meadows/pastures) and wasteland with grass vegetation form about 28% of the catchment area (Table 3, Figure 4).

**Figure 4.** Land use structure in the stream catchment.

#### 3.2.1. Characteristics of the Semi-Natural Section

Marshy areas without a distinctively formed bed can be found near the control point P01—the spring area. Valuable environmental elements include marshy meadows, as a spring area in the form of a weak outflow of groundwaters to the surface (seeps, bogs). The dominant land use forms include forest, trees/shrubs and grassland (marshy meadows) (Figure 5A). From cross-section P02, the bed-bottom material is classified as clay and silt, sometimes covered with a thin layer of mud (Figure 5B). The average water depth is about 3.0 cm. No modifications in the bed and natural morphologic elements were found. Soil (loess) is deposited on the slopes. The banks' cross-sections in the majority of the control profiles are eroded (as a result of water erosion), with vertical or underwashed banks of the stream. No modifications were found on the slopes, and the natural morphological elements include eroding bank undercuts with visible plant roots and numerous bank outwashes not fixed with vegetation (Figure 5C). There is no vegetation in the bed, and the vegetation appears on the slope as a uniform structure (Figure 5C) with ferns and grass vegetation, and as a complex structure on the slope tops with continuous tree stands (mainly maple, alder, oak, and spruce and grass vegetation). The stream valley in the semi-natural section is very shaded (50–75%). The land use in the areas near the banks is forest, trees and shrubs. The synthetic evaluation of the entire examined section (300 m) indicates an unobstructed flow in the bed, despite an occasional weak natural swelling by minute wood bed load (Figure 5D) and the vegetation in the bed and on the slope being typical of marshy areas (i.a. sedge). Detritus in the form of small-size dead organic matter from falling leaves and weak lateral tributaries can be found in the upper course of the stream. The laminar water flow dominates in the bed (this was the case when the studies were conducted, but rapid flow and overflow was also observed, mainly upstream in places where erosion faults/incisions formed small waterfalls). Florae were found in the stream valley, in the form of ivy (*Hedera helix*) and small balsam (*Impatiens parviflora* DC.) Faunae were also found, in the form of Roman snail (*Helix pomatia*) and slug. Debris classified as refuse/waste was also found the stream valley.

**Figure 5.** Selected elements of the hydromorphological diversity index in the semi-natural section: (**A**)—marshy areas and spring zone; (**B**)—laminar water flow; uniform vegetation on the slopes; (**C**)—eroding bank undercuts; bank outwash not fixed with vegetation; (**D**)—natural swelling by small wood bed load.

#### 3.2.2. Hydromorphological Characteristics of the Urbanised Section

In the initial control profiles of the urbanised section (P01–P03), the dominant bed material is fine sand and dust covered by a layer of muddy sediments. The water flow in the bed is laminar (Figure 6A). The water level in the bed is about 6 cm. The left-hand and right-hand banks, in the 1 m wide strip, are made of loose material (soil). In the 5 m strip from the bank, there are semi-continuous trees and shrubs, and a simple vegetation structure on the top and slopes of the banks; these are mainly alder (*Alnus*), maple (*Acer*), oak (*Quercus*) and spruce (*Picea abies*); there is also moss (*Bryophyta*), particularly on the tree branches and ferns (*Polypodiopsida*). No water plants were found in the bed. The bank profile ranges from gentle, which is classed as <45◦ (left-hand bank), to steep, which is >45◦ (right-hand bank). Typical land use in the strip 50 metres away from the bank top was trees/shrubs and buildings (school, private houses). Shading of the bed is visible, as well as exposed roots and outgrowth on the slopes of both banks. The natural erosion-related morphological elements include an undercut on the right-hand bank.

**Figure 6.** Selected elements of the hydromorphological diversity index in the urbanised section: (**A**) overgrown bed with laminar flow; (**B**)—outlet of Ø 1.00 m concrete culvert; (**C**)—grass vegetation (wasteland) and built-up areas; (**D**)—trees and shrubs. Grass vegetation (wasteland).

Farther in the urbanised section, the dominant bed-bottom materials are still fine sand and silt. The water flow is laminar. For about 30 metres, the stream flows in the pipe (the closed bed) with a paved surface over it (the car park). Then, the water flows out of the Ø = 1.00 m concrete pipe (the culvert) (Figure 6B), and the flow is <sup>Q</sup> ≈ 0.005 m3·s<sup>−</sup>1. There are no cross-section modifications in the open bed, and no natural morphological elements were found in the bed. The material of the right-hand and left-hand slopes is soil (sand and dust are dominating). No modifications on the righthand and left-hand slopes, and no natural morphological elements were found. The bed is heavily overgrown with vegetation; these are mainly plants above the water surface or recumbent amphiphites. In section P04–P07, the watercourse runs through uncultivated land. There is simple vegetation structure on the right-hand bank in the form of single trees/shrubs—mainly willow (*Salix*), hazel (*Corylus*) and oak (*Quercus*)—and green vegetation—mainly reed (*Phragmites australis*) and sedge (*Carex*) (Figure 6C,D). Uniform vegetation—mainly grass vegetation (*Phragmites australis)*—and tall herb vegetation mainly nettle (*Urtica dioica*)—dominate on the left-hand slope. Vegetation typical of marshy habitats and areas has also been found, i.e., purple loosestrife (*Lythrum salicaria*) and marsh horsetail (*Equisetum palustre*). The plot through which the stream flows is built up/developed on both sides. There is a road to the north with housing development along it, and single-family houses to the east (Figure 6C). Refuse/waste was found near the bank on the south side, behind the fence.

From profile P07 to P10, the stream bed is heavily overgrown with grass vegetation, with clear signs of anthropopression in the form of bottom- and slope- grading. The dominant bed-bottom material is fine sand. There is soil on the slopes. The water flow

in the bed is laminar. Farther along, the stream is pipelined (concreted) and flows under the urban development; before flowing into the Stradomka river, the stream bed is open again with visible signs of grading and vegetation mowing (inter-embankment zone of the Stradomka river). In addition, from among the elements recorded in the synthetic evaluation, gullies and culverts were found. The bed dimension for the semi-natural and urbanised section are presented in Figure 7A,B, respectively.

**Figure 7.** Bed dimensions: (**A**)—semi-natural section; (**B**)—urbanised section.

The HDS index indicates the hydromorphological diversity of the stream. The HDS values calculated for the individual study sections are varied from 27 for the urbanised section to 56 for the semi-natural section (Table 4). The HDS values for the semi-natural section are mainly influenced by: the presence of the natural morphological elements of the slopes in the control profiles (17.8%); the diversity of elements accompanying the trees (14.3%); and the naturalness and heterogeneity of the stream valley use (14.3%). For the urbanised section, the highest percentage of HDS value is found in the bed material heterogeneity (15.0%), the natural morphological elements of the slopes (15.0%), the vegetation diversity in the channel bed (15.0%), and the structure of the bank vegetation (15.0%). The HMS values reflecting the degree of anthropogenic changes in the hydromorphology of streams and rivers indicate that the lower part of the examined stream is the most modified (Table 4). The score for the lower-urbanised stream section is 15.5, indicating a significant modification of the stream related mostly to the modifications in control profiles (38.7%), especially profiling the bottom and slopes of the stream and the disturbance of the connectivity with the river valley (25.8%); the latter is especially influenced by urban areas. The upper examined section, for which the HMS is 0, is a hydromorphological habitat which has not been significantly modified. It is located within the tree stand boundaries and remains almost in semi-natural conditions.

It is unequivocally found that the value of the HIR index varies significantly for the tested sections. The final HIR value for the semi-natural section is 0.78, and is 44% higher in relation to the urbanised section (0.54); the changes in its numerical value mainly depend on the heterogeneity of the water flow, the natural morphological elements of the slopes, the diversity of elements accompanying the trees, and the natural land use of the valley (Tables 4 and 5). According to the classification [30], the hydromorphological condition in the semi-natural section is determined to be good, and in the urbanized section, to be poor. The hydromorphological condition of the whole stream is determined to be moderate (Table 5).

**Table 4.** Components of the Hydromorphological Diversity Score (HDS) and the Hydromorphology Modification Score (HMS).


**Table 5.** Results of hydromorphological evaluation and quality classes of individual examined sections and the whole stream.


#### **4. Discussion**

In the catchment area of the analyzed stream, an increase in the sum of annual precipitation and the average annual temperature is visible. Carpathian streams play an important role in draining the local catchment area, especially during flash floods, which are appearing more and more frequently [45]. The highest flood risks in Poland are characteristics of the Carpathian tributaries of the upper Vistula river, including the catchment area of the studied stream [46]. The negative impact of extreme weather phenomena is already visible, with drought and soil erosion bringing major problems, especially in the Carpathian Region; however, as reported by Ionita et al. [47], the frequency of droughts has not unusually increased in Central Europe when compared to preindustrial drought records. Natural climate change has led to changes in the water regimes of small upland and mountain streams, thus causing the occurrence of flash floods, or the drying up of watercourses in late summer [48]. The results of climate models show that in the area of Poland, in the near and far future, there will be an increase in air temperatures and in precipitation [49–51]. On the other hand, according to Atwood et al. [52], Climate models are quite poor when it comes to rainfall on a regional level. Any future trends in rainfall coming from climate models have to be treated with utmost caution. Attempts to identify climate change in the city micro-environment and improvement through water retention are conducted in Slovakia [53].

The physiographic (geometrical and morphometric) parameters of the catchment indicate a mountainous character of the examined stream; this includes a small area and significant elongation, a mean elevation of >200 m a.s.l., and a large slope [54]. The results of hydromorphological evaluation revealed significant differences between the semi-natural and urbanised part of the stream, indicating regulation and canalisation as the main reasons for hydromorphological degradation. Studying the catchment areas of two Carpathian flysch streams (the Jaszcze and Jamne), Bucała-Hrabia and Wiejaczka [54] concluded that a higher diversity of hydromorphological forms, with a corresponding HDS, occurs in the upper and medium course of the streams in the upland landscape; they also concluded that a lower diversity occurs in the lower part of the streams' catchment areas, which are subjected to anthropopression. A similar regularity of naturality and anthropogenic changes of the habitat between the upper and lower course of streams in Poland, India and China was shown by Kijowska-Strugała et al. [55].

The results of the scoring correspond to the preliminary hypothesis, according to which, despite the small stream, there is a noticeable difference between the analysed sections of the stream. An extensive amount of literature on linkages between the effects of anthropopression in stream habitats and river hydromorphology reports a number of relationships, at many levels. The need is urgent to develop refined and updated hydromorphological assessment systems targeting small ungauged watercourse evaluation, for use by the European Water Framework Directive (WFD) and national water-related policies. The topography and rich vegetation in the river's catchment area may make it difficult to perform field studies. Recent studies shows that Unmanned Aerial Vehicles (UAVs) can be used for inventory and assessment of the shydromorphological status of streams and rivers, especially in hard-to-reach areas. Drones and digital photogrammetry now provide an alternative approach for monitoring river habitats and hydromorphology [56,57].

The aim of the research currently carried out in Ireland is to advance knowledge on the role of small streams in water quality, biodiversity and ecosystem services protection; this will inform policy, measures and management options to meet water quality and other resources protection targets [58].

As a rule, natural watercourses have a higher HDS and the artificial watercourses (ditches, channels) have a higher HMS [21]. Hajdukiewicz et al. [6] indicate sbed regulation and significant grading (canalisation of the river) as the main causes of the hydromorphological degradation of the Biała river. Bryndal et al. [59] prove that high human activity in the catchment area deteriorates its natural drainage/water removal, thus contributing to a more frequent occurrence of high water stage and floods. Bedla et al. [60] and Pietruczuk, et al. [61] indicate that that the excessive piping of a river or stream adversely affects the natural environment in their valleys, deteriorating the aesthetic values of landscape. The diversity of morphological forms (flow type, bottom and slope material) plays a significant role in the biodiversity of fauna and flora in the watercourse bed and valley [28,62,63]. This can be seen in the examined stream, where the vegetation cover of the

watercourse body is much larger in the lower urbanised part of the stream than in the upper, semi-natural part (Table 4), characterised by a great slope and a larger dynamic of the flowing water. In the studied stream, it was demonstrated that, depending on the vegetation diversity, the HDS variation could change by about 25%. Furthermore, Kiraga [64] reported that, depending on the season, the vegetation diversity variation could change, which could lead to a step from one hydromorphological class to another. It was similar in the case of the HMS parameter, where modifications in the riverbed and stream valley were visible in the urbanised section. These changes significantly reduced the hydromorphological condition. It is well known that riparian vegetation plays a crucial role in sustaining river hydromorphological conditions [65]. Forests, trees and shrubs play an important part in the hydromorphological diversity of the zone near the stream banks. They act as biological protection (natural development) of the stream bed and valley, increasing the catchment retention, equalising the flows, and slowing down erosion of the stream slopes, bottom and banks. The root systems of these plants protect the soil from being washed out. Kałuza et al. ˙ [66] reported that, in the forested river reach, bottom and riverbank vegetation was completely absent. The river channel was narrow; however, due to the vicinity of trees and shrubs, considerable accumulation of organic matter was observed in the river channel. The flora acted as basic protection, mainly in the upper part of the catchment, where the land gradient was higher [62,67]. The majority of soils in the Wi´snicz Foothills, including the soils in the catchment of the examined stream, are made of Carpathian flysch (loess) which, without appropriate vegetation cover, are subject to intensive erosion [1,54].

Any transformation of the stream channel (riverbed) and catchment area should be thought out and included in the land use plan, as some of them might prove to be irreversible and have a destructive impact on the quality of the water environment [68–70].

#### **5. Conclusions**

The conclusions of the studies are as follows:

In its upper course, the stream has winding meanders, where the banks are subjected to erosion by flowing water, particularly during high water and high flow speed. In addition, it has low anthropopression and moderate diversity of morphological forms and conditions. Its lower course, on the other hand, has been significantly modified and transformed, in both the bed and the banks. The mouth part of the stream belongs to the intensively modified part of the surface waters. Factors significantly worsening hydromorphological conditions are various forms of anthropopressure, with particular emphasis on urbanization.

The evaluation indicates hydromorphological features of the stream which have been changed significantly at canalised sections, and which will most likely undergo the largest improvement when the bank protections are removed and the free channel migration is made possible, which—as a result of the progressing development of the lower part of the stream—is no longer possible.

This applied method can be used for the hydromorphological assessment and inventory of small, ungauged streams, because its scope includes elements that testify to the naturalness of stream habitats and a radical transformation in stream catchments. A major advantage of the HIR method is that it accounts for the HDS and HMS; this allows one to determine to what extent a given watercourse section is natural, and to what extent it is transformed. On the other hand, the final valorisation results are strongly affected by the presence of vegetation, which varies during the whole year. The type of water flow may also turn out to be an insufficient parameter—especially in periodic streams or those susceptible to rapid fluctuations in the water level after heavy rainfall.

This example of an examined stream indicates the need for more frequent monitoring the catchment areas of small streams not included in any hydrological classifications.

The authors believe that proper management of water resources in small, ungauged catchments lies with the local and regional authorities. All actions taken in the stream bed, and its catchment area, should be in accordance with the principle of sustainable development.

Based on this hydromorphological assessment, the results obtained helped us to evaluate the environmental changes and anthropogenic pressures on the stream sections; however, further research is required on the changes in the hydromorphological status of small watercourses.

**Author Contributions:** Conceptualization, Ł.B. and T.K.; methodology, Ł.B. and T.K.; software, Ł.B.; validation, Ł.B.; formal analysis, Ł.B. and T.K.; investigation, Ł.B. and T.K.; resources, Ł.B. and T.K.; data curation, Ł.B. and T.K.; writing—original draft preparation, Ł.B.; writing—review and editing, Ł.B. and T.K.; visualisation, Ł.B.; supervision, T.K.; project administration, Ł.B.; funding acquisition, T.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Department of Land Reclamation and Environmental Development, Faculty of Environmental Engineering and Land Surveying, University of Agriculture in Krakow.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


## *Article* **Use of Pedotransfer Functions in the Rosetta Model to Determine Saturated Hydraulic Conductivity (Ks) of Arable Soils: A Case Study**

**Łukasz Borek \*, Andrzej Bogdał and Tomasz Kowalik**

Department of Land Reclamation and Environmental Development, Faculty of Environmental Engineering and Land Surveying, University of Agriculture in Krakow, Al. Mickiewicza 24-28, 30-059 Krakow, Poland; andrzej.bogdal@urk.edu.pl (A.B.); tomasz.kowalik@urk.edu.pl (T.K.)

**\*** Correspondence: lukasz.borek@urk.edu.pl

**Abstract:** A key parameter for the design of soil drainage and irrigation facilities and for the modelling of surface runoff and erosion phenomena in land-formed areas is the saturated hydraulic conductivity (Ks). There are many methods for determining its value. In situ and laboratory measurements are commonly regarded as the most accurate and direct methods; however, they are costly and time-consuming. Alternatives can be found in the increasingly popular models of pedotransfer functions (PTFs), which can be used for rapid determination of soil hydrophysical parameters. This study presents an analysis of the Ks values obtained from in situ measurements conducted using a double-ring infiltrometer (DRI). The measurements were conducted using a laboratory permeability meter (LPM) and were estimated using five PTFs in the Rosetta program, based on easily accessible input data, i.e., the soil type, content of various grain sizes in %, density, and water content at 2.5 and 4.2 pF, respectively. The degrees of matching between the results from the PTF models and the values obtained from the in situ and laboratory measurements were investigated based on the root-mean-square deviation (*RMSD*), Nash–Sutcliffe efficiency (*NSE*), and determination coefficient (R2). The statistical relationships between the tested variables tested were confirmed using Spearman's rank correlation coefficient (rho). Data analysis showed that in situ measurements of Ks were only significantly correlated with the laboratory tests conducted on intact samples; the values obtained in situ were much higher. The high sensitivity of Ks to biotic and abiotic factors, especially in the upper soil horizons, did not allow for a satisfactory match between the values from the in situ measurements and those obtained from the PTFs. In contrast, the laboratory measurements, showed a significant correlation with the Ks values, as estimated by the models PTF-2 to PTF-5; the best match was found for PTF-2.

**Keywords:** soil; saturated hydraulic conductivity; pedotransfer function; Rosetta program; irrigation; climate change

#### **1. Introduction**

Water permeability is a key property of soils, especially with regard to the design of soil irrigation and drainage facilities, modelling of surface runoff and erosion phenomena in land-formed areas, and environmental processes occurring in porous media [1–3]. The study of soil physics, and in particular the determination of the filtration coefficient (Ks) based on direct methods, is a very interesting issue; nevertheless, it is both labour-intensive and costly [4–9].

Climate change in a region or environment entail changes in the method of soil cultivation and plant production. Changes in the interactions between agriculture (i.e., soil compaction) and natural environment (i.e., weather, soil conditions) are a key feature of the transitions which scientists are trying to explain. Therefore, learning the landscape is an important tool in the proper management of water resources in rural areas. Agricultural

**Citation:** Borek, Ł.; Bogdał, A.; Kowalik, T. Use of Pedotransfer Functions in the Rosetta Model to Determine Saturated Hydraulic Conductivity (Ks) of Arable Soils: A Case Study. *Land* **2021**, *10*, 959. https://doi.org/10.3390/ land10090959

Academic Editor: Claude Hammecker

Received: 10 August 2021 Accepted: 8 September 2021 Published: 10 September 2021

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

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

development depends on the temperature and rainfall distribution during the growing season and on the fertilities and properties of soils, including their abilities to conduct and retain water. The climatic and soil conditions are often unfavourable for the rational and efficient management of environmental resources. Moreover, global climate change, as characterised by an increased frequency of periods of excesses and shortages of water, may hinder the development of modern agriculture, i.e., by reducing the quantity and quality of crop yields. To counteract these adverse phenomena, water engineering and land reclamation facilities are being used in many regions of the world [10,11]. However, to correctly design the spacing for drains, drainage elements, and irrigation ditches and/or the technical parameters of irrigation devices influencing the intensity of the water supply, it is crucial to identify the water permeability of the corresponding soils, including their infiltration and filtration processes [12,13].

The saturated hydraulic conductivity (i.e., Ks) is a key input parameter for modelling the water flow in a soil [14]. Its value can be determined using in situ and laboratory methods, but this often turns out to be a difficult task, owing to the large spatial and temporal variabilities of soil properties [15], and the need to take special care of samples (e.g., to provide an intact structure [16]. Ks is a basic physical property of the soil that affects all soil-plant-water relationships and processes. However, it is one of the most variable soil properties, as it is related to the soil texture and structure, and is influenced by factors such as the land topography, vegetation, land use, and climate [17]. It is widely accepted that in situ measurements provide the most accurate results for the determination of a soil's physical and hydraulic properties.

There are different ways to determine the movements of water in unsaturated (infiltration) and saturated (filtration) zones that are of interest to the scientific community and engineers in various industries. The high cost and effort of in situ measurements and the growing demand for such data have led scientists to seek alternative indirect methods [18,19]. In recent years, pedotransfer functions (PTFs) have become popular and have been used to estimate time-consuming and difficult-to-measure soil properties, such as Ks or soil water retention (pF) values [1,20–22]. Ahuja et al. [23], Rawls et al. [24], Timlin et al. [25], and Suleiman and Ritchie [26] proposed formulas for estimating Ks values based on the effective porosity. Cosby et al. [27], Wösten et al. [28], Saxton and Rawls [29], and Weynants et al. [30] proposed formulas for determining Ks values based on using the textural composition, density, and other physical properties of soils. Currently, studies demonstrating the possibility of using neural networks to estimate Ks values [5,31] are very popular, along with studies aiming to improve existing formulas with appropriate amendments [32]. Puckett et al. [33] proposed a model to predict Ks based on only claysized particles. The authors showed that fine sand, sand, and clay percentages were highly correlated with Ks. Jabro [4] developed a model that used the site-average dry bulk density and grain size as predictive variables of Ks. A simple Campbell's model, with soil texture data, can be used [34] or Smettem and Bristow's model from soil clay content using a variety of agricultural topsoil samples [35,36] published a model which were used the relationships between soil texture and soil moisture content at saturation and soil texture and Ks.

Pedotransfer functions, as a five hierarchical models, offered by Rosetta are one of the most applied PTFs, as demonstrated by a number of citations in other research papers more than 1000 times (in Google scholar website). One of the reasons for the Rosetta PTFs popularity is its easy availability and an extensive soil database from North America and Europe containing 2134 soil samples with water retention data and measurements for the saturated hydraulic conductivity (Ks) for 1306 of soils samples and continues to be developed [37]. The majority of work on the performance of PTFs in Rosetta with the aim of determining Ks factor has focused on surface layers of the soil (0–30 cm or 0–50 cm) [13,14,38].

In this work concentrated on the full depth soil pits (150 cm) and including an estimate of Ks for surface layers (0–50 cm) and additionally for a deeper levels of soil (50–150 cm). It was considered that for engineering purposes, such as irrigation or drainage planning, it is advisable to know the Ks parameter.

The aim of this research was to test the possibility of using PTFs in the Rosetta program to determine the Ks values of arable soils based on easily accessible input data, i.e., the soil type, granulometric composition, density, and water content. We hypothesize that the Ks values obtained from PTFs in the Rosetta program for arable soils correspond to the values obtained from (i) direct in situ and (ii) laboratory measurements. To test the hypothesis we examined the arable soils in Central Europe (district of Racibórz, south Poland) in a 150 cm depth soil pits.

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

#### *2.1. Description of Study Area*

The study area is located in the southwestern part of the Silesian Voivodeship (Poland) in the district of Racibórz (Figure 1). According to Kondracki's [39] division of Poland into physio-geographical regions, the study area is located in the Central European Plain (31), in the macroregion of the Silesian Lowlands (318.5), and on the border of two mesoregions: the Głubczyce Plateau (318.58) and Racibórz Basin (318.59). The main watercourse flowing through the region is the Oder River, along with its left tributary, the Psina River.

**Figure 1.** Location of the study area.

#### *2.2. Meteorological Conditions*

In terms of climate, the study area is considered as one of the warmest areas in this region. In general, Poland has a mostly temperate climate, in transition between an oceanic climate dominating in the north and west of the country, and a continental climate in the south and east. In the multiannual period from 1971–2000, the mean annual air temperature was 8.5 ◦C, and total precipitation was 616 mm (according to the Institute of Meteorology and Water Management—station in Racibórz). Higher mean temperatures and a decrease in rainfall have been observed over the last decade. The characteristics of the meteorological conditions are listed in Table 1.


**Table 1.** The characteristics of meteorological conditions of the study area.

To assess the meteorological conditions of the region where the study was conducted, Selyaninov's hydrothermal coefficient (*HTC*) values were calculated, as a measure of the precipitation efficiency in a given month. The calculation was performed as follows [40]:

$$HT\mathbb{C} = \frac{10 \cdot P}{\Sigma t} \text{ (-)}\tag{1}$$

In the above, *P* is the monthly sum of precipitation (mm), and Σ*t* is the sum of the mean daily air temperature values in a given month (in ◦C).

Based on the value of the *HTC*, the months of the growing season (April-September) were classified as follows: *HTC* ≤ 0.4—extremely dry; 0.4 < *HTC* ≤ 0.7—very dry; 0.7 < *HTC* ≤ 1.0—dry; 1.0 < *HTC* ≤ 1.3—relatively dry; 1.3 < *HTC* ≤ 1.6—optimal; 1.6 < *HTC* ≤ 2.0 moderately humid; 2.0 < *HTC* ≤ 2.5—humid; 2.5 < *HTC* ≤ 3.0—very humid; *HTC* > 3.0 extremely humid.

Based on the *HTC* values calculated for the months of the growing season from the 2010–2019, and as shown in Table 1, it was determined that 58.4% of the months were dry (8.3% extremely dry, 13.3% very dry, 28.5% dry, and 8.3% relatively dry); in 13.3% of the months, the conditions were optimal; and 28.3% of the months were humid (8.3% moderately humid, 6.7% humid, 8.3% very humid, and 5.0% extremely humid). In general, the variable meteorological conditions in the study area, characterised by periods with excesses and shortages of water, indicate the need to use drainage and irrigation facilities for agriculture. To determine their technical parameters, information on the soil permeability is needed. This confirms the advisability of conducting research on the possibility of using mathematical models to determine the soil Ks values, as an alternative to time-consuming in situ or laboratory tests.

#### *2.3. Field Measurement and Soil Sampling*

The research was conducted on arable lands during the agricultural season from 2012 to 2015. The field soil tests were conducted at 16 measurement-control points (up to a depth of 150 cm) in Wojnowice (2 points), Bojanów (2 points), Owsiszcze (6 points), Strzybnik (2 points), and Tworków (4 points). Three to five genetic levels were identified for each soil profile.

The soil infiltration was measured in situ, at depths of 10 cm (topsoil) and 35 cm (subsoil), using a double-ring infiltrometer (DRI) method (Figure 2). The DRI comprised an inner ring (9.5 cm diameter) and outer ring (19.5 cm diameter) inserted into the ground at a depth of 10 cm. The DRI was inserted by using a falling weight-type hammer striking on a wooden plank placed uniformly on top of the ring, and without undue disturbance to the soil surface. Each ring of the DRI was filled with a constant head of water level, and the outer ring helped when checking the lateral flow from the inside ring, so as to better estimate the infiltration, reducing losses. The Ks value was estimated when the water flow rate inside the inner ring reached a steady state [1], which in the case of the studied soils, lasted approximately 3–4 h. The infiltration rate was calculated for the respective time intervals Δ*T* as follows [41]:

$$\dot{q} = \frac{864 \cdot 4 \cdot V}{3.14 \cdot D\_r^2 \cdot \Delta T} \text{ (m} \cdot \text{day}^{-1}) \tag{2}$$

where *V* is the volume of water (cm3) added to the inner ring at time Δ*T* (s), and *Dr* is the diameter of the inner ring (cm2).

For steady infiltration (as a constant value for the soil), the Land and Water Development Division [42] developed infiltration classes, as follows: very slow—<0.024 m·day<sup>−</sup>1; slow—0.024 ÷ 0.12 m·day−1; moderately slow—0.12 ÷ 0.48 m·day−1; moderate—0.48 ÷ 1.56 m·day−1; moderately rapid—1.56 ÷ 4.20 m·day−1; rapid—4.20 ÷ 5.81 m·day−1; and very rapid—>5.81 m·day<sup>−</sup>1.

Undisturbed soil samples were taken from each genetic horizon, using cylinders with a volume of 100 cm<sup>3</sup> (three replicates). In addition, approximately 1.0 kg of disturbed soil from each genetic horizon was used to determine the soil texture, and other laboratory analyses were performed.

**Figure 2.** The double-ring infiltrometer (DRI) (photo. Ł. Borek).

*2.4. Laboratory Analysis*


$$BD = \frac{M\_s}{V\_t} \text{ (g\text{-}cm}^{-3})\tag{3}$$

In the above, *Ms* is the mass of the dry soil weight (g), and *Vt* is the volume of the total soil sample (cm3).

• The soil water retention was investigated based on determining the soil suction using ceramic plates in a 5/15 bar pressure plate extractor. The pressure plate equipment used in this study was manufactured by the American Soil Moisture Equipment Corporation. In engineering practice, the soil suction is usually calculated in units of pF, as follows:

$$\text{pF} = \log|\text{h}| \,\text{ }\tag{4}$$

Here, h is soil suction pressure (in cm H2O).

In the laboratory, the soil water potentials were measured at 300 cm H2O ≈ 33 kPa ≈ 2.5 pF (representing the field capacity) and 15,000 cm H2O ≈ 1500 kPa ≈ 4.2 pF (representing the permanent wilting point) [45]. The soil volumetric water content (*θV*) at pF = 2.5 and 4.2 was determined based on the *gravimetric method.* In particular, it was calculated as the ratio of the amount of water in the soil sample to the dry weight of the soil, after drying in an oven at 105 ◦C for approximately 18–24 h, as follows [44]:

$$
\theta\_V = \frac{V\_{H\_2O}}{V\_t} \left( \text{cm}^3 \cdot \text{cm}^{-3} \right) \tag{5}
$$

In the above, *VH2O* is the volume of water in the soil sample (cm3), and *Vt* is the volume of the total soil sample (cm3).

• The Ks was measured under laboratory conditions using a laboratory permeameter (Figure 3) and using the Darcy's law [46] with constant head method (Equation (6)) and falling head method (Equation (7) on undisturbed soil samples for three replications. The constant method can be used with virtually any soil, apart from poorly permeable soils such as clay, whereas the falling-head method is used to measure low-permeability soils, such as f.i. clay or peat samples [47,48]. The soil samples were placed in the laboratory permeameter and then saturated in water for 2–3 days. For the purpose of this study, a permeameter produced by the Eijkelkamp with a closed or open system and 25 holders was used. The Ks was calculated using the Darcy's [46] equation, as follows:

$$V = K \cdot i \cdot A \cdot t \tag{6}$$

Here, explanations below;

Determining Ks using the constant head method was calculated from the following equation:

$$\text{Ks} = \frac{V \cdot L}{A \cdot t \cdot h} \tag{7}$$

Here, *V* is the volume of water flowing through the sample (cm3); *K* is the permeability coefficient or 'Ks' (cm·day−<sup>1</sup> or m·day−1); *<sup>h</sup>* is the water level difference inside and outside ringholder or sample cylinder (cm); *L* is the length of the soil sample (cm); *i* is the permeability rise gradient or *h*/*L* (–); *A* is the cross-sectional area of the sample (m2); and *t* is the time used for flow through of water volume (day). During measuring the following parameters have been determined: *L* and *A*—constants, depending on the type of sample ring used; *V*—volume measured in the burette (1 mL = 1 cm3); *t*—length of time lapse; *h*—calculated with the water levels measured with the water level meter.

Determining Ks using the falling head method was calculated from the following equation:

$$\mathcal{K}\_{\mathbf{s}} = \frac{a \cdot L}{A \cdot (t\_2 - t\_1)} \cdot \ln \frac{h\_1}{h\_2} + \frac{\mathbf{x} \cdot a \cdot L}{A \cdot \sqrt{(h\_1 \cdot h\_2)}} \tag{8}$$

Here, Ks is the permeability coefficient (cm·day−<sup>1</sup> or m·day−1); *<sup>h</sup>* is the water level difference inside and outside ringholder or sample cylinder (cm) *h*<sup>1</sup> and *h*<sup>2</sup> water level difference inside and outside the ringholder at respectively *t*1(start) and *t*<sup>2</sup> (end); *A* is the surface of a cross-section of the sample (cm2); *L* is the length of the soil sample (cm); *t* is the time between beginning and end of the measuring *t*<sup>2</sup> − *t*<sup>1</sup> (day); *a* is the cross-section surface of a ringholder or sample cylinder (cm2) for a sample cylinder applies *A* = *a*; *x* is the evaporation factor (literature value): 0.0864 cm·day−<sup>1</sup> or 0.000864 m·day<sup>−</sup>1.

The permeability of the soil is also determined by viscosity of the soil solution. Viscosity depends on the temperature. The laboratory water temperature varies from 18 to 22 ◦C, whereas the average groundwater temperature is 10 ◦C. Therefore, for certain applications the permeability will have to be corrected for the viscosity of the soil solution (usually water). The Ks for viscosity was corrected using the following equation [49]:

$$K\_{10} = K\_T \cdot \frac{h\_T}{h\_{10}} \tag{9}$$

Here, *<sup>K</sup>*<sup>10</sup> is the corrected Ks at 10 ◦C (cm·day−<sup>1</sup> or m·day−1); *KT* is the Ks at the applied temperature (cm·day−<sup>1</sup> or m·day<sup>−</sup>1); *<sup>h</sup>*<sup>10</sup> is the dynamic viscosity of water at 10 ◦<sup>C</sup> (Pa·s); *hT* is the dynamic viscosity of water at T ◦C (Pa·s).

The soil classification was established according to the Polish Soil Classification [50], World Reference Base for Soil Resources (FAO and International Union of Soil Sciences (IUSS) Working Group World Reference Base (WRB), [51]), and USDA soil taxonomy (Soil Survey Staff, [43]).

**Figure 3.** The laboratory permeameter produced by the Eijkelkamp (photo. Ł. Borek).

#### *2.5. Rosetta Description*

The Ks values were calculated using the published neural network program Rosetta, with hierarchical PTFs based on five levels of input data [5]. The first level (PTF-1) used soil textural classes based on a lookup table providing parameter means for each USDA textural class. The second level (PTF-2) used the sand, silt, and clay percentages as inputs, and in contrast to PTF-1, provided hydraulic parameters that varied continuously with the texture. The third level (PTF-3) included the predictors used in the level PTF-2, along with the soil dry-bulk density. The fourth level (PTF-4) used PTF-3 and the soil volumetric water content (*θ*) at a water suction of 33 kPa (2.5 pF). The last level (PTF-5) comprised all of the other parameters, plus the *θ* value at a water suction of 1500 kPa (4.2 pF) [38,52]. The necessary input data for the Rosetta program (e.g., % soil texture group and bulk density) were determined based on the laboratory analyses of soil samples from the different soil horizons (Tables 2–4).

In the Rosetta program, the relationship between *θ* and the water suction (*h*), i.e., the water retention *θ*(*h*), as well as Ks, are described using the well-known Mualem–van Genuchten equations [53] and are given as follows:

$$\theta(h) = \theta\_{\bar{r}} + \frac{\theta\_{\bar{s}} - \theta\_{\bar{r}}}{\left[1 + \left(\alpha |h|^{\bar{n}}\right)^{m}}\right]^{M} \tag{10}$$

$$K(h) = K\_{\mathbf{s}} \frac{\left\{ 1 - (\mathbf{a} \cdot \mathbf{h})^{n-1} [1 + (\mathbf{a} \cdot \mathbf{h})^{n}] \right\}^{m}}{\left[ 1 + (\mathbf{a} \cdot \mathbf{h})^{n} \right]^{m/2}} \tag{11}$$

In the above, *<sup>θ</sup>*(*h*) is the soil volumetric water content (cm3·cm−3) at suction *<sup>h</sup>* (cm H2O); *<sup>θ</sup><sup>s</sup>* and *<sup>θ</sup><sup>r</sup>* are the saturated and residual water contents (cm3·cm−3) at *<sup>h</sup>* = 0 and 15,000 cm H2O, respectively; *α* (>0 in cm<sup>−</sup>1) is related to the inverse of the air entry suction; *n* (>1) is a measure of the pore-size distribution, and *m* = 1 − 1/*n*; and Ks is the saturated hydraulic conductivity (m·day<sup>−</sup>1), as mentioned above.


#### **Table 2.** Selected soil physical properties over all depths.


**Table 2.** *Cont.*

\* soil texture group according to USDA (1999): SL—sandy loam; SCL—sandy clay loam; L—loam; CL—clay loam; SiL—silt loam; SiCL—silty clay loam; C—clay.


**Table 3.** Basic descriptive statistics for soil physical properties for topsoils and subsoils (*n* = 32).

**Table 4.** Basic descriptive statistics for soil physical properties over all depths (*n* = 68).


#### *2.6. Statistical Analysis and Model Performance Evaluation*

The data set comprised the analytical results from the soil samples collected from the 16 soil pits. For the statistical analysis, the procedures provided by the program *Statistica PL Version 12.5* were used, with a 5% significance level. The minimum and maximum values were determined for each physical parameter of the soil, and the arithmetic mean, median, standard deviation (*SD*), and coefficient of variation (*CV*) were computed. Moreover, a Spearman correlation test was conducted for the dataset. The correlation strength and direction of the relationship between two variables were determined based on Spearman's rank correlation coefficient (*rho*), which takes values from −1 to 1. A positive sign of the coefficient indicated the existence of a positive correlation, and vice versa. The closer the values were to −1 and 1, the stronger the correlation. When *rho* = 0, there was no correlation between the examined variables [54]. This statistical method was chosen after finding that most data were not normally distributed (Shapiro-Wilk test). All of data are presented in the tables and graphs, i.e., to produce a visual image that is helpful in interpreting the results.

Based on the available soil data in this study, five widely used PTFs were selected to assess their respective performances in estimating Ks, via comparison with the measured Ks values (via the DRI) for the soil cores collected in the field (using a laboratory permeameter). The calculated values of Ks obtained using Rosetta were compared to the corresponding measured values and evaluated using two statistical parameters: rootmean-square deviation (*RMSD*) and Nash–Sutcliffe efficiency (*NSE*). The *RMSD* gives the mean difference between the measured and calculated values of Ks, and is calculated as follows [38,55]:

$$RMSD = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (m\_i - p\_i)^2} \tag{12}$$

where: *mi* is the measured value of Ks, and *pi* is the corresponding (predicted) value of Ks as obtained using the Rosetta program.

The *RMSD* values are always non-negative and should be as low as possible; *RMSD* = 0 indicates a perfect fit of the model to the measurement data.

The *NSE* compares measured and predicted values (here, for Ks), and is given as follows [56]:

$$NSE = 1 - \left[\frac{\sum\_{i=1}^{n} (m\_i - p\_i)^2}{\sum\_{i=1}^{n} (m\_i - \overline{m}\_i)^2}\right] \tag{1.3}$$

In the above, *n* is the number of observations, *mi* and *pi* are the measured and predicted values of Ks, respectively, and *mi* is the mean of the measured Ks values.

The *NSE* values range from −∞ to 1. A value of *NSE* = 1 corresponded to a perfect match of calculated values to the measured values of Ks. The 'efficiency *NSE'* = 0 indicated that the calculations of Ks using the Rosetta program were as accurate as the mean of the measured data; in contrast, an efficiency *NSE* < 0 occurred when the measured mean was a better predictor than the model or, in other words, when the residual variance, as described by the numerator in the equation above, was larger than the data variance described by the denominator.

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

According to the PTG [50], FAO and IUSS Working Group WRB [52], and USDA soil taxonomy [43], the examined soils were classified as follows.

• Order 3: Brown forest soils (brown earths, PTG—Polish: *Gleby brunatnoziemne*; WRB: Cambisols; USDA: Inceptisols—Udepts), Type 3.1. 'Euthrophic brown soils' (PTG— Polish: Gleby brunatne eutroficzne; WRB: Haplic Cambisol, Haplic, Stagnic, Endogleyic, or Vertic Cambisol (Eutric); USDA: Typic or Humic or Aquic or Oxyaquic or Vertic Eutrudepts)—occurred in Strzybnik.


As shown in Table 2 and Figure 4, the examined soils are heterogeneous in terms of their textures. Seven textural classes are observed in the examined soil profiles: silt loam (SiL; *n* = 42), sandy loam (SL; *n* = 10), loam (L; *n* = 6), clay (C; *n* = 3), sandy clay loam (SCL; *n* = 3), clay loam (CL; *n* = 2), and silty clay loam (SiCL; *n* = 2). The silt loam is the predominant soil texture. The *SD* values (Tables 3 and 4) for most soil properties are large, indicating moderate to strong variability; this is not uncommon for soils [12,57,58]. Typically, the *CV* values are < 45%, indicating the mean variability of the measurement data. The highest *CV* values are observed for sand (57.0%) and clay (57.2%). The sand content in the soil profiles is between 11% and 76%, the silt content is between 16% and 72%, and the clay content is between 8% and 65% (Tables 2 and 4). The high swelling-clay content could be a cause of the poor water flow through the samples, as well as a consequence of abnormal Ks results [59].

**Figure 4.** Textural composition of soils investigated.

With regard to the surface soil horizons (*n* = 32), the mean soil bulk density (ρb) is 1.55 g·cm<sup>−</sup>3, whereas that for overall depths is 1.60 g·cm−<sup>3</sup> (Tables <sup>3</sup> and 4). The soil volumetric water contents at the matric potentials of 33 kPa (2.5 pF) and 1500 kPa (4.2 pF) in both cases have similar mean values (0.350/0.351 cm3·cm−<sup>3</sup> and 0.152/0.155 cm3·cm<sup>−</sup>3, respectively).

Higher values of Ks are found in the upper horizons of the tested soils, owing to the significant macroporosity caused by plant roots, reconsolidation [60], and annual loosening from agrotechnical works and alternating soil freezing and thawing [61]. In addition, zoogenic channels (especially earthworm channels) and plant residues increase soil permeability [62]. For these reasons, the samples taken from arable and sub-arable horizons are less dense than those from lower horizons; the latter are also subject to self-compacting (Tables 2–4).

Table 5 presents descriptive statistics concerning the values of Ks for the topsoils and subsoils (*n* = 32), as obtained from in situ and laboratory measurements, and as estimated with the use of the PTF models in the Rosetta program (variant 1). In contrast, Table 6 presents statistical measures considering the Ks values of all of the soil horizons (*n* = 68), but excludes the in situ tests not conducted on larger and deeper horizons, for practical reasons (variant 2).


**Table 5.** Basic descriptive statistics for soil hydraulic conductivity (Ks) for topsoils and subsoils.

**Table 6.** Basic descriptive statistics for soil hydraulic conductivity (Ks) over all depths.


In variant 1, the characteristic values of Ks (maximum, mean, and median values) as obtained from the in situ measurements are approximately four times higher than the Ks values obtained from the laboratory measurements. Such large discrepancies between the results obtained from both direct methods result from the manner of conducting the measurements; in particular, small-volume soil samples, which do not allow for the inclusion of all of the factors influencing the soil permeability, were taken for the laboratory tests. The mean Ks value obtained from the in situ measurements, i.e., 0.85 m·day−1, is four times higher than those estimated by PTF-1 and PTF-2, almost half lower than that obtained from PTF-3, and twice as high as the values obtained from PTF-4 and PTF-5. In addition, the median value from the in situ tests (0.40 m·day−1) is from 1.2 to 2.2 times higher than the values obtained from the four PTF models—the exception is model PTF-3, whose median value is four times higher than that from the in situ tests. The maximum and mean values of the Ks values obtained from laboratory tests, compared to those from the in situ method, are usually more similar to those obtained from the PTF models (Table 5). The mean Ks values obtained from the in situ measurements and PTF-3 indicate the mean class, whereas the results obtained from laboratory measurements and the other four PTF models indicate a moderately slow basic infiltration class [42].

Most of the subsurface horizons are characterised by low Ks values, usually not exceeding 0.50 m·day−1. In variant 2, covering the analysis of all soil horizons, the mean value of Ks from the laboratory measurements (0.21 m·day−1) is very close to the results obtained from the PTF-1 and PTF-2 models, almost six times lower than the Ks value from PTF-3, and three times lower than the values estimated by PTF-4 and PTF-5 (Table 6). The mean Ks values obtained from the laboratory measurements and PTF-1, PTF-2, and PTF-3 indicate the moderate class, whereas the results obtained from the other three PTF models indicate a moderately slow basic infiltration class [42].

Regardless of the amount of data taken for analysis (*n* = 32 or 68), the values of Ks are usually characterised by high or very high random variability (*CV* = 26.1–230.7%), not only in comparison to the other soil hydrophysical properties, but also compared to the methods used to determine Ks, i.e., field, laboratory, and PTF tests (Tables 5 and 6, Figure 5). This regularity has been described in various other studies, for example, Rezaei et al. [63]. As reported by Merdun et al. [64] and García-Gutiérrez et al. [65], the high spatial variability between the soil horizons with regard to Ks are the result of the soil heterogeneity, which hampers the prediction of Ks using mathematical models.

**Figure 5.** Box plots showing the performance of Ks to different measurement methods.

The Ks values obtained from the laboratory measurements are statistically significant (*p* < 0.05), negatively correlated with the clay fraction and density, and positively correlated with the sand fraction (Table 7). Similar observations were noted by Zhao et al. [13], who conducted studies on the Loess Plateau of China. The interdependence of Ks values and soil textures were well documented by Hillel [66], García-Gutiérrez et al. [65], and other researchers [3]. In the case of in situ measurements, the soil hydraulic conductivity is only significantly positively correlated with the sand fraction (*rho* = 0.29). The PTF-1 model is significantly positively correlated with the soil density, but negatively correlated with the clay fraction and soil volumetric water content (*θ*) at water suctions of 33 and 1500 kPa (2.5 pF and 4.2 pF). The PTF-2 and PTF-3 models show a significantly positive correlation with the sand fraction and negative correlations with the clay fraction and soil volumetric water content (*θ*) at 2.5 and 4.2 pF. The PTF-4 and PTF-5 models are significantly positively correlated with the sand fraction, and are negatively correlated with the dust fraction (Table 7).


**Table 7.** Spearman's rank correlations coefficients (*rho*) between measured and predicted saturated hydraulic conductivity (*Ks*) and other selected soil physical properties.

Red color means, that the determined correlation coefficient is significant (for the significance level α = 0.05).

No significant correlations were detected between the Ks values obtained from the in situ measurements and the PTF models. However, with the exception of the PTF-1 model, the Ks values from the laboratory tests as estimated by the four PTF models are positively correlated (Table 7). For these reasons, only the linear correlations between the values obtained from the laboratory tests and PTF models are presented in Figure 6a–e.

In general, the matching of the Ks results obtained from the direct in situ measurements with the results estimated by the Rosetta PTFs is poorer than that obtained from the laboratory measurements and estimated by the Rosetta PTFs. This is confirmed by the *RMSD* values, which range from 1.40 to 1.64 in the first case, whereas in the second case, they range from 0.19 to 1.58. In addition, the values of the coefficient of determination (*R*2), which indicates the extent to which the model explains the gathered measurement data, are usually much higher when comparing the data obtained from the laboratory tests and that from the PTF models (Table 8).


**Table 8.** Statistical indices for estimation of saturated hydraulic conductivity (Ks).

The *RMSD* values of the Ks results obtained from the laboratory measurements and estimated with PTFs reveal the best match for the case of the PTF-2 model. The *NSE* and *R<sup>2</sup>* values are the best in this case (Table 8), indicating that the granulometric composition (% of sand, silt, and clay) in the PTF-2 model has the greatest influence on the water permeability of the soils.

**Figure 6.** *Cont*.

**Figure 6.** The laboratory measured and predicted saturated hydraulic conductivity Ks for the five pedotransfer function models invistigated (PTF).

#### **4. Conclusions**

The search for alternatives to labour-intensive and cost-intensive in situ tests and laboratory methods for determining the value of Ks, at a time of global climate change, seems to be a priority. The increased frequency of periods with excesses and scarcities of water will require the greater use of soil irrigation and drainage facilities, and the modelling of surface runoff and erosion phenomena in land-formed areas. Correspondingly, knowledge regarding soil permeability is required for all technical and non-technical activities.

In general, when using a descriptive comparative analysis, it is difficult to show a significant relationship between Ks values from in situ tests and the diverse hydrophysical properties of arable soils, which are subject to cyclic loosening, reconsolidation, and/or freezing and thawing. The high sensitivity of this physical parameter to biotic and abiotic factors, especially in the soil upper horizons, hampered a satisfactory adjustment to the Ks values obtained from the PTF models. Therefore, in situ tests, although difficult and timeconsuming, should continue to be conducted to obtain a reliable measurement database that is free of errors resulting from the soil type variability, and that considers the influences of fauna and flora on the physical properties of soil. This is of particular importance for

the modelling of surface runoff and erosion processes and design of surface irrigation systems (e.g., sprinklers), where only 40–50 cm of the upper horizon of the soil is wetted, and the technical parameters of the corresponding drainage facilities are highly dependent on soil permeability.

The in situ measurements of Ks are significantly correlated only with the laboratory tests; nevertheless, the values obtained in situ are much higher, and are influenced by a number of factors, including those of a methodological nature (e.g., the larger diameter of the measuring cylinder). In contrast, the laboratory measurements show a significant correlation with the Ks values estimated using the models PTF-2 to PTF-5. The best match is found (based on the *RMSD*, *NSE*, and *R*<sup>2</sup> values) when using the PTF-2 model in the Rosetta program. This means that the most reliable Ks values can be obtained from the percentages of the sand, dust, and clay fractions. Therefore, there is a rational alternative, in the form of the PTF, for the costly and labour-intensive determinations of Ks conducted in the laboratory. The use of the PTF-2 model allows for the determination of the hydraulic conductivity of deeper soil horizons, where the variability of the hydrophysical properties is lower and there are large methodological limitations in conducting field work or difficulties in taking intact structural samples for laboratory tests. This type of solution can be used to collect data for designing subsurface drainage and irrigation facilities, some of the technical parameters of which are determined by the Ks values.

In conclusion, the research results confirmed the remarkable difference between the performance of PTFs in the Rosetta and the Ks results obtained from field or laboratory measurements. The study thus provided evidence that the Ks values obtained from local dataset give better results performances in predicting the soil saturated hydraulic conductivity, compared to PTFs in the Rosetta derived from large datasets, which in turn confirm the limitation of applying PTFs developed from one region to other regions.

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

**Funding:** This research was funded by the Department of Land Reclamation and Environmental Development, Faculty of Environmental Engineering and Land Surveying, University of Agriculture in Krakow.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to acknowledge the authorities of the Agro-Industrial Company 'AGROMAX' in Racibórz for providing agricultural areas and agricultural equipment for the field study.

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

#### **References**


## *Article* **Concept of Soil Moisture Ratio for Determining the Spatial Distribution of Soil Moisture Using Physiographic Parameters of a Basin and Artificial Neural Networks (ANNs)**

**Edyta Kruk \* and Wioletta Fudała**

Department of Land Reclamation and Environmental Development, Faculty of Environmental Engineering and Land Surveying, University of Agriculture in Krakow, Al. Mickiewicza 24-28, 30-059 Kraków, Poland; wioletta.zarnowiec@urk.edu.pl

**\*** Correspondence: edyta.kruk@urk.edu.pl

**Abstract:** The results of investigations on shaping the soil moisture ratio in the mountain basin of the M ˛atny stream located in the Gorce region, Poland, are presented. A soil moisture ratio was defined as a ratio of soil moisture in a given point in a basin to the one located in a base point located on a watershed. Investigations were carried out, using a TDR device, for 379 measuring points located in an irregular network, in the 0–25 cm soil layer. Values of the soil moisture ratio fluctuated between 0.75 and 1.85. Based on measurements, an artificial neural network (ANN) model of the MLP type was constructed, with nine neurons in the input layer, four neurons in the hidden layer and one neuron in the output layer. Input parameters influencing the soil moisture ratio were chosen based on physiographic parameters: altitude, flow direction, height a.s.l., clay content, land use, exposition, slope shape, soil hydrologic group and place on a slope. The ANN model was generated in the module data mining in the program Statistica 12. Physiographic parameters were generated using a database, digital elevation model and the program ArcGIS. The value of the network learning parameter obtained, 0.722, was satisfactory. Comparison of experimental data with values obtained using the ANN model showed a good fit; the determination coefficient was 0.581. The ANN model showed a minimal tendency to overestimate values. Global network sensitivity analysis showed that the highest influence on the wetness coefficient were provided by the parameters place on slope, exposition, and land use, while the parameters with the lowest influence were slope, clay fraction and hydrological group. The chosen physiographic parameters explained the values of the relative wetness ratio a satisfactory degree.

**Keywords:** soil moisture; physiographic parameters of basins; artificial neural network (ANN); redundancy analysis (RDA)

#### **1. Introduction**

Soil moisture [1,2], especially in the context of climate change, is one of the greatest problems in hydrological responses [3,4], agriculture and on industrial sites [5–9]. The question of the influence of initial soil moisture, despite being mentioned by many studies as a significant erosion factor, has rather rarely been extensively researched [10]. The actual state of soil moisture determines the beginning of surface runoff. Surface runoff is shaped by many phenomena, including exceeding a soil's full saturation and filtration possibilities, subsurface outflow, and groundwater filtration [11,12]. The spatial distribution of soil moisture is realized by various methods [13], such as: direct measurements; the use of models based on physiographic basin parameters [14] as well as remote sensing data [15,16]; and in many cases, a combination of the above methods; and has been the subject of interest of many authors.

Gómez-Plaza et al. [17] carried out investigations of moisture distribution in a small basin located in Spain, using 50 sample points over a 20 m × 20 m regular network.

**Citation:** Kruk, E.; Fudała, W. Concept of Soil Moisture Ratio for Determining the Spatial Distribution of Soil Moisture Using Physiographic Parameters of a Basin and Artificial Neural Networks (ANNs). *Land* **2021**, *10*, 766. https://doi.org/ 10.3390/land10070766

Academic Editor: Krish Jayachandran

Received: 22 May 2021 Accepted: 17 July 2021 Published: 20 July 2021

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

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

They used a TDR (time domain reflectometry) device to measure the soil moisture at depths of 15 cm. They connected the spatial distribution of soil moisture with use of TWI (topographic wetness index). They also pointed out a use for these indicators for rather wet climates. TWI is calculated by the following equation:

$$TWI = \ln\left(\frac{\alpha}{\tan\beta}\right) \tag{1}$$

where: *α* is the specific catchment area defined as the local upslope area draining through a unit contour length, which is equal to the grid cell width in this study; and *β* is the local slope gradient [18].

Tombul [19], simulating surface runoff depending on soil initial moisture, introduced shaping of the spatial distribution in the upper soil layer. The investigations were carried out in the Kurukavak river basin, located in Turkey, in a 4.25 km2 area. Soil moisture measurements were carried out by means of a TDR device. For investigation of spatial distribution, he used the Xinanjiang model [20], based on basic hydrological responses such as: evapotranspiration, runoff generation and flow routing, connecting soil moisture distribution with wetness index (WI):

$$\theta\_{WI} = \theta\_s \cdot \frac{1 - \left(1 - \frac{\theta}{\theta\_\*}\right)^{\frac{1}{1+b}}}{1 - (1 - MI)^{\frac{1}{5}}} \tag{2}$$

where: *θWI* is the critical value of capacity, *θ<sup>s</sup>* is the moisture at full saturation, *θWI* ≤ *θs*, *θ* is the actual moisture, *b* is the shape parameter.

For his investigations, he used high-resolution DEM, constructed based on a topographic map, using the vector model TIN (triangulated irregular network). The introduced distribution was substituted by a topographic index.

Zhang et al. [21] researched the spatial distribution of soil moisture at the 10-cm layer, on a 9.8-ha experimental plot, using a 15 m × 15 m network comprising 250 points in total. They showed that moisture was characterized by the normal distribution. Spatial distribution was examined using kriging, using the module Spatial Analyst connected in the program ArcGIS release 9.0. They showed the significance of investigations of moisture distribution for fertilization and irrigation needs.

Merdun et al. [22] studied the spatial distribution of water retention, for three various initial moistures (dry, mean, and wet), obtained using a rainfall simulator on soil at the experimental centre of the University of Sutcu Imam in Turkey, using 100 sampled points. The structure of the spatial distribution and the semivariogram were executed in the program GS+ 5.1.

Penna et al. [23] investigated moisture distribution in the upper soil layer, in a 1.9 km<sup>2</sup> area located in a mountain basin in the Italian Alps, using 42 sample points. They showed the significance of relief parameters and atmospheric conditions, including slope, TWI, exposition, and solar radiation for the shaping of this property, showing a 42% share of these factors to explain moisture distribution.

Temimi et al. [24] studied the distribution of soil moisture in the Mackenzie River basin, over a 4000 km<sup>2</sup> area located in Canada. In their investigations, they used a modified, not-statistical topographic wetness index (TWI), determined based on DEM:

$$TWI = \ln(a) - \ln\left(a \cdot \tan(\beta)\right) \cdot e^{-u \cdot LAI} \tag{3}$$

where: *a* is the alimentation area, *β* is the local slope, *LAI* is the leaf area index, *μ* is the coefficient of radiation extinction, depending on plant cover, with values between 0.35 and 0.70.

DEM was determined in the framework of the program Shuttle Radar Topography Mission (SRTM). LAI was in turn determined from satellite observations.

Fan et al. [25] carried out investigations on spatial variability of soil moisture, using a regular 50 m network with 20 sample points, on an experimental field in China. They used kriging, realized in the program ArcGIS release 9.0.

Nikolopulos et al. [11] carried out investigations on the influence of initial soil moisture in the Fella river basin, including an area of 623 km2 and their subbasins, located in the Italian Alps, using peak flow and outflow discharge during a rainstorm. For investigation of the spatial distribution of soil moisture, they used the hydrologic model Triangulated Irregular Network (TIN), based on the Real-Time Integrated Basin Simulator. They used a DEM with a resolution of 20 m.

Jia et al. [26] examined the temporary stability of spatial distribution of moisture in loess soil in China, on five experimental plots of dimension 61 × 5 m, in 10 cm intervals, up to a depth of 1 m. They showed that the vertical distribution of moisture does not have a unitary trend; however, the horizontal distribution was similar for the same layers, irrespective of a density increase.

This work aims to present the concept of determination of soil moisture distribution, based on physiographic parameters of a basin, realized by the use of artificial neural networks (ANNs).

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

#### *2.1. Study Sites*

The study was carried out in the Outer Western Carpathians, in the southern region of Małopolska Province, Poland (Figure 1). The M ˛atny Stream (1.47 km2) flows into the Mszanka River in the hamlet of Skiby at 20◦9 2.35 E, 49◦37 30.52 N. The main watercourse starts at 20◦08 28.52 E, 49◦36 25.44 N. The mean annual air temperature is 7.4 ◦C. The lowest temperature of −25.8 ◦C was recorded on 3 February 2012 and the highest, 33.0 ◦C, on 8 July 2013. The long-term annual average total precipitation is 846.63 mm and the total precipitation over the years in the study was 948.1 mm. The highest daily precipitation was 100.7 mm on 15 May 2014. The snow usually starts to melt at the beginning of March. The growing season starts around 10 April, and the winter begins around 30 November.

**Figure 1.** Location of the research area within Poland.

The terrain comprises low- and medium-height mountains with peak heights ranging from 617.6 m above sea level (m.a.s.l.) to 732.0 m.a.s.l. The lowest point is situated 490.0 m.a.s.l. The mean height of the catchment is 582.66 m.a.s.l. [27]. The slope distribution is as follows: <5%, (0.06 km2); 5–10% (0.27 km2); 10–18% (0.67 km2); 18–27% (0.31 km2); and >27% (0.16 km2). The weighted average slope for the entire catchment was 16.28% [28]. The catchment land use structure was dominated by grassland (73.5%); arable lands constituted 14.3% and included the following crops: spring oats (*Avelana sativa*), 7.3%; potatoes (*Solanum tuberosum*), 4.3%; and common wheat (*Triticum aestivum*), 2.7%. Forests accounted for 9.5% and urban areas for 2.7% of the catchment land (Figure 2). The catchment area is cut by a network of dirt roads. Most of them are deeply furrowed and tend to transform into water-carrying streams during and after rain events.

**Figure 2.** Photos showing typical land use of the M ˛atny basin.

Soils in the catchment area are diverse, depending on slope location. The soil cover in the M ˛atny stream catchment is dominated by loamy soils, including sandy clay loam, loam, silt loam, clay loam and sandy loam. Pedological conditions were identified by analysis of a 1:25,000 agricultural soil map and categorized into the respective groups according to USDA standards [29].

#### *2.2. Experimental Design/Models Applied*

Calculations implementing the soil moisture ratio (wetness coefficient) model algorithm were performed using the interface of ArcGIS 10.3.1 software [30]. The spatial data were obtained based on the Polish State Surveying Coordinate System (PUGW). The following parameters were determined: altitude, slope, flow direction, exposition, shape of the slope, situation on the slope, hydrological group, use, and clay fraction. Thematic layers were developed using the sources specified in Table 1.


#### **Table 1.** Description of parameters used in coefficient wetness model.

where: <sup>1</sup> is *S* = <sup>Δ</sup>*<sup>h</sup> <sup>l</sup>* ·100 ; <sup>2</sup> is *SD*<sup>8</sup> <sup>=</sup> *maxi*<sup>=</sup>1,8 *<sup>Z</sup>*9−*Zi <sup>h</sup>*∅(*i*) .

Three sampling time periods were adopted for each of the three growing seasons: dry (period I), medium (period II) and wet (period III) AMC (Antedescent Moisture Conditions), distinguished based on the sum of rainfall from the previous 5 days (mm) and the season category (dormant or growing) (Table 2). In botany and agriculture, the growing season is defined as the portion of each year when native plants and ornamental plants grow, while the dormant season is when growth and development are temporarily stopped [33]. In much of Europe, the growing season is defined as the average number of days with a 24-h average temperature of at least 6 ◦C; in southern Poland, this typically lasts from April to September.

**Table 2.** AMC classes [34].


The soil moisture ratio *Kw,i* in point *i* was determined based on the concept of wetness coefficient [35,36]:

$$K\_{w,i} = \frac{\theta\_i}{\theta\_b} \tag{4}$$

where: *<sup>θ</sup><sup>i</sup>* is the soil volumetric water in the given point *<sup>i</sup>* within the basin area [m3·m−3], and *<sup>θ</sup><sup>b</sup>* is the soil volumetric water content in the basal point [m3·m<sup>−</sup>3].

From this, the soil moisture in a point *i* in a basin can by determined as:

$$\theta\_{\bar{i}} = \theta\_{\bar{b}} \cdot K\_{w,\bar{i}}(p\_1, p\_2 \dots p\_n) \tag{5}$$

where: *Kw,i*(*p*1, *p*<sup>2</sup> ... *p*n) is the relative wetness coefficient determined based on the following parameters: place on a slope, exposition, land use, shape of slope, altitude, flow direction, slope, clay fraction, hydrological group.

The purpose of the wetness coefficient was to assess the distribution of soil moisture in the catchment area, based on the measurement from a basal point located on a hilltop. Soil moisture was measured in 2014, between 11:00 pm and 14:00 am during days without and with rainfall using the TDR device, as a mean value in the 0–0.10 m layer, for 379 measuring points (153 for period I, 100 for period II and 126 for period III). The distribution of points was random, taking into account the specific character of the mountain catchment. The measurements were made, for period I, on 25 July (the sum of rainfall from the previous 5 days was 21.2 mm, in the vegetative growing season); for period II, on 4 October (the sum of rainfall from the previous 5 days was 21.7 mm, in the growing season); and for period III, on 25 September (the sum of rainfall from the previous 5 days was 34.8 mm, in the dormant season) (Figure 3).

**Figure 3.** Location of moisture measurement points.

#### *2.3. Statistical Analysis and Data Procedure*

An ANN (artificial neural network) model in the form of a multilayer perceptron (MLP) was used to generate wetness coefficients, using Statistica software (release 12.5). A total of 70% of all variables were applied for the learning process; 15% were used for validation and 15% to test the model. A quasi-Newton algorithm with a Broyden–Fletcher– Goldfarb–Shanno (BFGS) modification was selected for the learning neural network. The sum of squares (SOS) was treated as the error function. The artificial neural network model was used to establish the association between the wetness coefficient, whereas a physiographic parameter, indicative of the soil and its use, was applied in all data set as an independent variable. The multi-layer perceptron consisted of three layers of neurons: (1) an input layer, (2) an output layer and (3) intermediate (hidden) layers. Each neuron had a number of inputs (from outside to the subsequent layer or out of the network) [37]. In this study, the network system included an input and a hidden layer made of nine neurons (altitude, slope, flow direction, exposition, shape of the slope, situation on a slope, hydrological group, land use and clay fraction) and an output layer with one neuron (wetness coefficient).

Studies on the percentage effect of selected soil, use and physiographic parameters were carried out using ANN based on a global network sensitivity analysis [38–40].

RDA (redundancy analysis) of standardized environmental variables was performed to explain and describe the pattern of variability in a parameter [41] Multivariate analysis showed the presence of two main gradients of environmental variables. Positions of the vectors of independent variables were proportional to the loading factors. The multivariate analysis was carried out with Canoco for Windows version 4.51.

The analysis of empirical model adjustment to experimental data was carried out using the following metrics [42], presented in Table 3.


**Table 3.** Measures of model performance [43,44].

where: <sup>1</sup> is *MEP* = <sup>1</sup> *<sup>n</sup>* · <sup>∑</sup>*<sup>n</sup> i*=1 *cm <sup>i</sup>* − *c i* ; <sup>2</sup> is *RMSE* = 1 *<sup>n</sup>* · <sup>∑</sup>*<sup>n</sup> i*=1 *cm <sup>i</sup>* − *c i* 2 ; <sup>3</sup> is *MPE* = <sup>1</sup> *<sup>n</sup>* · <sup>∑</sup>*<sup>n</sup> i*=1 *i cm i* ·100%; <sup>4</sup> is *ME* <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>∑</sup>*<sup>n</sup> i*=*n*(*c<sup>m</sup> <sup>i</sup>* −*c p i* ) 2 ∑*n <sup>i</sup>*−*<sup>n</sup>*(*c<sup>m</sup> <sup>i</sup>* <sup>−</sup>*c*) <sup>2</sup> ; *c<sup>m</sup> <sup>i</sup>* is the measured values, *c p <sup>i</sup>* is the computed values, and *n* is the number of data.

Spatial variability of the investigated and measured values of the soil moisture ratio was determined using kriging. The kriging method allows to obtain the most probable values in any part of an investigated area and to find the location for new measuring points. To briefly define the technique of kriging, it is a method for optimizing the estimation of the spatially correlated quantity Z, both in stationarity and non-stationarity instances.

#### **3. Results**

Moisture at the base point was 0.146 m3·m−<sup>3</sup> in the dry season (period I), 0.247 m3·m−<sup>3</sup> in the medium season (period II), and 0.433 m3·m−<sup>3</sup> in the wet season (period III). The measured moisture fluctuated between 0.086 m3·m−<sup>3</sup> and 0.226 m3·m−<sup>3</sup> for period I, between 0.146 m3·m−<sup>3</sup> and 0.354 m3·m−<sup>3</sup> for period II, and between 0.145 m3·m−<sup>3</sup> and 0.441 m3·m−<sup>3</sup> for period III. The obtained values were in accordance with the ones introduced in the literature [45] The ANN model was then generated based on the results of the investigation of physiographic parameters (Figure 4).

**Figure 4.** Physiographical parameters of the M ˛atny stream basin.

The best-fitting model turned out to be the MLP 9-4-1, with four perceptrons in the hidden layer. The MLP 9-4-1 model is characterized by high fitting quality and low error, and therefore presents a good adjustment (Table 4). Thanks to the global sensitivity analysis (Table 5), the relative and absolute influence of a number of parameters on the soil moisture ratio Kw were determined.


**Table 4.** Analysis of quality and errors of the MLP 9-4-1.

**Table 5.** Global sensitivity analysis of the MLP 9-4-1 model.

Perceptron activation functions


hidden output

validation 0.004

In the M ˛atny stream basin, the parameters place on the slope (40%) and exposition (12%) had the highest impact (Table 5). Values of the soil moisture ratio fluctuated between 0.75 and 1.85. Based on the simulated values of the soil moisture ratio for every one of the measured points, the map of the spatial distribution of this parameter was generated using Surfer 10 and ArcGIS, based on data generated by ANN 9-4-1 (Figure 5a) and on the measurement date (Figure 5b). The highest values of the soil moisture ratio occurred in the northwest part of the basin. Figure 6 presents a comparison between the ANN simulated values of the soil moisture ratio and the one elaborated based on measured data. Values of the simulated wetness coefficient ranged between 0.89 and 1.13. Model efficiency measures for the MLP 9-4-1 were as follows: MEP: 0.004, RMSE: 0.104, MPE: −0.6%, ME: 0.580 and R2: 0.581 (Table 6).

**Figure 5.** Spatial distribution of the soil moisture ratio in the M ˛atny basin generated using: (**a**) the MLP 9-4-1 model, and (**b**) based on measured data elaborated by kriging technique.



**Figure 6.** The measured versus simulated soil moisture ratio.

Multivariate analysis (RDA) was employed to explain the relative importance of particular explanatory variables and underline differences between parameters. The patterns of the independent variables for soil factors and environmental parameters are shown on the plot (Figure 7). We noted a positive correlation among the flow direction, slope, clay content and wetness coefficient. Generally, multivariate analysis for the studied parameters indicated a strong positive correlation among the soil hydrologic group. Moreover, we saw a correlation between altitude and place on the slope. On the other hand, only exposition N, exposition E and exposition W were strongly negatively correlated to other examined parameters. The first component described 27.43% of the total variation, and the second component explained 42.01% of the total variance among the study parameters.

**Figure 7.** Multivariate RDA analysis for the soil moisture ratio is created by vectors (red squares are for convex; purple squares for concave slopes).

#### **4. Discussion**

In this study, the concept of relative wetness coefficient was developed as a base for determination of the spatial distribution of soil moisture on an area of a mountain basin.

The idea of the coefficient is based on the observations that there is a relation between soil moisture in a point on a slope and on a flat in a basin. Values of relative wetness coefficient were obtained by the ANN model generated based on nine basin parameters: place on a slope, exposition, land use, shape of slope, altitude, flow direction, slope, clay fraction and hydrological group and direct measurements of soil moisture on the area of small Carpathian basin (1.47 km2). Parameters of ANN learning, testing, validation, and model efficiency measures were very satisfactory. The model overestimated prognosed values only of 0.6%. Model was fitted to the experimental data very well. Based on the ANN model, the sequence of significance of the particular basin parameters were arranged to provide the relative wetness coefficient explanation. The highest influence had: place on slope, exposition, and land use. In turn, clay fraction and the hydrological group had no significance, because the investigated basin is not highly differentiated regarding soil. A very interesting model for evaluation of the relative wetness coefficient foe small basin was presented by Svetlitchnyi et al. [33]. It was based on basin parameters such as: the shape of slope, distance from the divide, slope aspect and overall length of a slope. The obtained model errors were very satisfactory. The relative wetness coefficient can be used as an easy way to determine the distribution of soil moisture in the upper layer in small mountain basins.

#### **5. Conclusions**

The soil moisture ratio, as a relation between soil moisture in a given point and in a basal point located on a flat, generated based on basin parameters, by use of artificial neural networks is one of the simplest ways to determine the distribution of soil moisture on an area of small mountain basin. The main basin parameters that can influence the relative soil moisture ratio are the place on a slope, exposition, land use (if it is not uniform on an area of a basin), shape of slope and altitude. Our investigations showed that they controlled the relative soil moisture ratio in sum in 83%. The model for determination of the relative soil moisture ratio, generated by use of the ANN, gave very satisfied results in terms of errors, and explained as many as 58% of cases. The ANN model overestimated prognosis values only of 0.6%. For larger basins and basins with highly differentiated soil, the next studies have to be carried out.

**Author Contributions:** Conceptualization, E.K. and W.F.; methodology, E.K.; software, E.K.; validation, E.K. and W.F.; formal analysis, E.K.; investigation, E.K. and W.F.; resources, E.K. and W.F.; data curation, E.K. and W.F.; writing—original draft preparation, E.K.; writing—review and editing, E.K. and W.F.; visualization, E.K. and W.F.; supervision, E.K.; project administration, E.K.; funding acquisition, E.K. and W.F. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Department of Land Reclamation and Environmental Development, Faculty of Environmental Engineering and Land Surveying, University of Agriculture in Krakow.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


## *Article* **Rate of Fen-Peat Soil Subsidence Near Drainage Ditches (Central Poland)**

**Ryszard Oleszczuk 1, Ewelina Zaj ˛ac 2,\*, Janusz Urba ´nski <sup>3</sup> and Jan Jadczyszyn <sup>4</sup>**


<sup>3</sup> Institute of Civil Engineering, Warsaw University of Life Sciences—SGGW, Nowoursynowska 166 St., 02-787 Warsaw, Poland; janusz\_urbanski@sggw.edu.pl


**Abstract:** This study analyzed design depths (*to*), post-subsidence depths (*t*), shallowing magnitudes (*d* = *to* − *t*) and ratio values (*d*/*t*) of 12 drainage ditches in a fragment of the drained Solec fen-peat (central Poland) over a period of 47 years between 1967 and 2014. A significant decrease of the designed depth of the ditches *to* was shown, from the average designed value of 0.97 m to their average depth after subsidence, *t* = 0.71 m. The ratio (*d*/*t*) of 0.41, which is associated with the degree of organic matter decomposition, indicated medium degree of peat decomposition. The average values of bank and bottom subsidence of the ditches during the analyzed period, 1967–2014, were 0.43 m and 0.17 m, respectively. The values of the average annual rate of land surface subsidence in the vicinity of the ditches were varied and within the range of 0.09 cm year−<sup>1</sup> to 1.70 cm year<sup>−</sup>1, with an average of 0.92 cm year<sup>−</sup>1. Two linear empirical equations were proposed to calculate the amount of subsidence and the average annual rate of subsidence of peat soil surface near the drainage ditch route, based on the knowledge of the initial thickness of the peat deposit. The results of calculations using the equations proposed by the authors were compared with calculations of the same parameters using 10 equations published in the literature. The results obtained using the proposed equations were mostly larger than those calculated with literature-published equations.

**Keywords:** drained peat soils; ditch subsidence; empirical equations of peat subsidence; drainage ditches

#### **1. Introduction**

In excessively humid areas (e.g., river valleys), the high groundwater table and low air content in the soil are conducive to the formation of peat soils, which are characterized by a high organic carbon content. Peat deposits covering only approximately 4% of the world are estimated to accumulate approximately 30% of the global organic carbon of all soils in their deposits [1–3]. In Europe, about 16% of peat soils are used for agricultural purposes (arable land and grassland), including the vast majority in Western European countries [4]. The area of organic soils and soils of organic origin in Poland estimated on the basis of the soil-agricultural map at a scale of 1:25,000 occupies 13.3% of agricultural land, wherein 5.6% are peat soils, 1.6% are sited peat soils, and 6.1% are moorsh soils [5]. These areas, used by many disciplines (including agriculture, horticulture, and forestry), have been drained in numerous cases [6]. This has resulted in the interruption of the existing organic matter accumulation process due to incomplete decomposition under conditions of very high humidity and lack of air access. In their natural state, these soils absorb carbon dioxide as a result of the assimilation and photosynthesis processes of the vegetation living in these areas. The change in air–water relations that accompanies the drainage of peatlands

**Citation:** Oleszczuk, R.; Zaj ˛ac, E.; Urba ´nski, J.; Jadczyszyn, J. Rate of Fen-Peat Soil Subsidence Near Drainage Ditches (Central Poland). *Land* **2021**, *10*, 1287. https://doi.org/ 10.3390/land10121287

Academic Editors: Wiktor Halecki, Dawid Bedla, Marek Ryczek and Artur Radecki-Pawlik

Received: 11 October 2021 Accepted: 18 November 2021 Published: 24 November 2021

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

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

has caused the initiation of mucking and mineralization processes in organic matter, and the conversion of organic forms of nitrogen into mineral compounds and emission of greenhouse gases, such as carbon dioxide and nitrous oxide, into the atmosphere [7–9]. The ongoing climate change, significant increase in air temperature, and reduction in precipitation during the growing season lead to extreme weather phenomena, as well as increased frequency of agricultural drought [10–12]. The foregoing weather phenomena contribute to the intensification of mucking processes and permanent degradation of drained peat soils. Changes occur in their physical properties, i.e., retention, hydraulic properties, and soil compaction due to the shrinkage process, as do changes in their chemical properties, specifically reduction in organic carbon content and increase in ash content. All these phenomena contribute to the lowering of the surface (shallowing of the profile) in drained peat deposits, i.e., subsidence. This leads to a decrease in the thickness, as well as surface area, of these deposits and a permanent reduction in their agricultural production potential. It may even contribute to the total disappearance of peat soils in the natural environment and their transformation into typical mineral soils in the long run. This is an extremely unfavorable phenomenon, taking into consideration the numerous functions that peatlands play in the natural environment (e.g., accumulation of organic carbon, water retention, and biodiversity).

In the past, using peatlands for agricultural purposes (mainly meadows and pastures) required lowering the groundwater level to reduce moisture in the root zone (0–30 cm) in view of the water needs of grassland vegetation. These drainages were done in stages by constructing a network of deep draining channels, usually followed by a dense network of drainage ditches after several years (typically spaced about 100 m apart) and, as needed, an additional network of drainage pipes (spaced about 30 m apart). The foregoing drainage network parameters depend on the drainage depth (usually about 0.8–1.0 m) and the hydraulic properties of the peat soils, specifically the filtration coefficient [13]. Drainage of these soils for agricultural purposes, as opposed to drainage of mineral soils, has resulted in many negative physical, chemical, and biological processes, due to a decrease in the moisture content of the top layers, and a consequent increase in air content. Among the physical processes, the shrinkage process is observed due to a decrease in the moisture content of these soils, which consequently leads to increased compaction, reduced retention capacity, and changes in the hydraulic properties of these soils [14–19]. The disappearance of the buoyancy force and the pressure of the drained peat layers on the underlying layers results in the lowering of the surface of the drained peatland, which assumes the highest annual rate in the initial 15–20 years after drainage and is called the first subsidence phase [20–26]. The magnitude and rate of peatland surface subsidence depends mainly on the physical properties of the soil, peat type, and drainage depth. Despite the stabilization of the soil water level at a certain reduced level, the rate of the soil surface subsidence decreases but does not cease and passes into the second phase of subsidence [20,22,27–29]. Mainly chemical and biological processes occur in Phase II, resulting in partial (humification) and complete (mineralization) decomposition of organic matter. During these processes, the resulting greenhouse gas, carbon dioxide, is emitted from the atmosphere and dissolved in groundwater [2,7,30–36]. Consequently, this leads to a decrease in organic carbon in the surface peat layers, and an increase in nitrogen, both in the soil and in the water into which it is leached. This leads to upsetting the natural C:N ratio and causing it to narrow, which is a measure of the degradation processes of these soils. As a result of degradation, the mineral content, i.e., ash content, increases [16,29,37–39].

The processes of surface subsidence in drained peat soils, the decrease of their thickness, and even disappearance of these areas from the natural environment have been quite well described in the world literature based on field studies around the world at single measurement points. These phenomena have been monitored over many years at single measurement points, along cross sections, or on the surface of entire peatlands [27,28,40–46]. For example, for drained peat soils located in the Groot-Mijdrecht Polder near Amsterdam (The Netherlands), surface subsidence measurements were made at 1423 locations

between 1954 and 1968 [43]. However, less attention has been paid to areas located near drainage channels and ditches where intensification of the surface subsidence process of drained soils is the greatest due to the lowest groundwater table in their vicinity during the drainage process [37,47–49]. In the literature there is much less research related to the effects of described processes on functioning and technical parameters of drainage systems, such as the decreasing depth of ditches and position of drainage pipelines as a result of the subsidence process. Undoubtedly, this creates a threat to the proper functioning of such systems and the necessity of their modernization in terms of conducting proper water management on such areas (irrigation of these areas).

Measurement results for banks and bottoms of 12 ditches located in a drained fen peat are analyzed in this paper and available archival data from the project of this object are used. Based on 1967 (designed) and 2014 (measured) elevation data, the following objectives are undertaken:


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

The study was conducted on a fragment of the Solec fen-peatland (Góra Kalwaria Commune, Piaseczy ´nski Poviat, Masovian Voivodeship, Central Poland (52◦2 16.345 N, 21◦6 14.622 E). The 220-ha peat deposit is made up of quaternary formations up to about 50 m thick. The aquifer consists of sands, and in the upper part there are organic soils, sedge, and sedge-reed peats of medium degree of decomposition [50,51]. In the vicinity of the site boundaries, the peat thickness is approximately 30–40 cm, while in its central part it varies within approximately 130–180 cm, and locally up to 250 cm. The first land draining works at the site were conducted between 1941 and 1943 and consisted mainly of a new section of the Mała River running through the center of the site [52] (Figure 1a). In 1967, a project for a drainage system was created that called for the construction of 62 ditches and about 80 communication and damming structures [52]. Draining works according to the project were carried out in 1967–1971. The site was subdivided into 13 plots in which sub-irrigation system was implemented. The main source of irrigation water was the Mała River and feeder A and R-32 (Figure 1b). Most of the site was in agricultural use as meadows and pastures. Since about 1985, irrigation has no longer been implemented, no damming devices are in place, and the ditches only drain adjacent land, discharging into the Mała River. Most of the area has not been used for agriculture since about 2000.The Mała River drains water in the range of about 100 m on each side. In the growing season the ditches are usually dry, and after short-term intensive rainfall the depth of water in the ditches is about 15–20 cm.

**Figure 1.** Location of the study site, its schematic (**a**), and the analyzed parts of it (**b**).

The study covered a fragment of the Solec fen with an area of about 50 ha, on which there are 12 drainage ditches: R-15, R-17, R-19, R-21, R-22, R-23, R-24, R-25, R-26, R-26a, R-27, and R-29. Ditch spacing in most cases is 90 m, while ditches R-22 and R-25 are spaced 140 m apart, with lengths ranging from 250 m to about 540 m (Figure 1b). Geodetic measurements of the ordinates of both their banks and bottom were made in cross sections located every 100 m on the ditches in 2014, and the thickness of the peat deposit on both banks was measured. Measurements were taken with a NI 050 levelling device from Carl Zeiss Jena, with reference to the existing state geodetic network.

Changes in bank and bottom ordinates of ditches from 1967 to 2014 were analyzed and ditch depths as designed (*to*) (1967), ditch depths after subsidence (*t*) (2014), ditch shallowing magnitudes (*d* = *to* − *t*), ratio (*d*/*t*), and subsidence magnitudes/average annual rate of land surface subsidence near ditches were calculated. The values of the *d/t* ratio were included in ranges depending on the degree of peat decomposition [53,54]. For peat with a low degree of decomposition (*H1*–*H4*), the values of this ratio range from 1.28–0.57, for a medium degree of decomposition (*H5*–*H6*), the ratio (*d/t*) is 0.39, and for peat with a high degree of decomposition (*H7*–*H10*), it is 0.27–0.19. Based on the results, two empirical relationships were developed between the amount of subsidence (*y*) (cm) and the rate of average annual subsidence (*z*) (cm year−1) of the surface of the drained peat deposit and its initial thickness (*x*) near the drainage ditches. The calculation results obtained using the proposed equations were compared with the calculation results of these quantities obtained using 10 empirical equations developed for the second phase of subsidence of large areas (not only for near-ditch locations) of drained fens published in the literature, summarized in Table 1.


**Table 1.** Empirical equations for the second phase of peat surface subsidence based on initial peat deposit depth and soil surface subsidence or annual peat subsidence rates published by various authors.

#### **3. Results**

The example longitudinal profile of the selected ditch R-29 (Figure 2) and the crosssections (Figure 3) show the archival (1967) averaged ordinates of both ditch banks and its bottom, and the results of measurements of the same parameters in 2014, as well as the ordinates of the mineral subsoil. The profiles of 12 ditches included in the drainage system [52] were analyzed. Based on the ordinates of the banks and bottom of the ditches at each hectometer, it was found that the design depth of the ditches (*to*) was within 0.73–1.35 m, with an average of 1.0 m. This indicates a medium-to-high intensity of drainage in the area [55]. The wide range of ditch depths was probably due to the slope of the terrain and the level of the designed bottom. As a result of the dewatering process, it was observed that the depth of the ditches decreased after subsidence (*t*) in 2014 to an average depth of about 0.74 m, and their shallowing averaged about 0.26 m (Table 2 and Table S1). The values of the ratio of the magnitude of ditch shallowing (*d*) to its depth after subsidence (*t*) vary within very large limits from about 0.01 to 1.16, with an average of 0.40. During the analysis period of 47 years (1967–2014), significantly higher values of subsidence were observed on the banks of the ditches compared to the subsidence of their bottoms. The ordinates of the ditch banks decreased by about 0.43 m on average (0.04–0.80 m) compared with the 1967 ordinates, whereas the ditch bottom decreased to a much smaller extent (by about 0.17 m on average), i.e., about 40% in comparison with the amount of their bank subsidence (Table 2 and Table S1). This is also an effect of siltation, the contribution of which is difficult to estimate.

#### **Figure 2.** Longitudinal profile of the R-29 ditch.

**Figure 3.** Cross-sections of the R-29 ditch in the year 1967 (designed) and 2014 (measured).


**Table 2.** Mean values and range (min–max) parameters describing subsidence of drainage ditches (banks and bottom).

The variability of the initial peat thickness in 1967 along the route of the analyzed ditches was quite large and ranged from about 30 cm to about 255 cm, with an average of about 122 cm (Table 3 and Table S2). The deposit was the shallowest along the route of ditch R-25, and the deepest along ditches R-19 and R-21. As a rule, the highest deposit thicknesses were found in the middle part of the trench route, while the lowest thicknesses were found at the ends of the trenches, near feeders A and R-32.

**Table 3.** Mean values and range (min–max) of initial peat depth and average rate of subsidence of drainage ditches.


Explanations: L—total length of the ditch.

Based on the magnitude of bank subsidence of the analyzed 12 ditches over 47 years (1967–2014), it was found that the highest subsidence values were observed in the central part of their route, slightly smaller at the mouth of the ditches feeding into the Mała River, and the lowest near the end of the ditches, which may be related to the variable thickness of the deposit along the route of the ditches (Table 3 and Table S2). The amount of land surface subsidence in the vicinity of the ditches during the analysis period ranged from about 4 cm to 80 cm with an average of 43 cm. This corresponds to a reduction in the thickness of the peat deposit relative to the initial thickness of 4% to 75% (43% on average). The mean annual subsidence rate during the 47-year period analyzed was also variable, ranging from 0.08 to 1.7 cm year<sup>−</sup>1, with an average of 0.92 cm year−<sup>1</sup> (Table 3 and Table S2). The results of land surface subsidence measurements from 1967 to 2014 along the analyzed drainage

ditches are presented in Figure 4a in relation to the initial thickness of the peat deposit. The average annual subsidence rate of the peat surface depended on the same parameter, as shown in Figure 4b.

**Figure 4.** Relationship between the amount of surface subsidence of a drained peat deposit (**a**), and the average annual rate of peat deposit surface subsidence (**b**) on the initial thickness of the peat deposit.

Based on the results obtained (Figure 4a,b), two empirical Equations (11) and (12) were developed, which make the amount of subsidence of the drained peat deposit surface (*y*) (cm) near the ditches and the average annual rate of peatland surface subsidence (*z*) (cm year<sup>−</sup>1) dependent on the initial thickness of the peat deposit (*x*) (cm) (explanation of symbols and units used in the text):

$$y = 0.243\mathbf{x} + 13.71\tag{11}$$

$$z = 0.0052x + 0.292\tag{12}$$

The relationships developed are linear, with correlation coefficient *r* values of 0.65 in both cases. Using the proposed Equations (11) and (12), calculations were performed for the amount of subsidence (*y*) (Figure 4a) and the average annual subsidence rate (*z*) (Figure 5b) depending on the initial thickness of the peat deposit (*x*) within the range of its values from 30 cm to 250 cm (the range of the actual thickness of the deposit in the analyzed fragment of the site). The values of the same parameters (*y*) and (*z*) were calculated using Equations (1)–(10) given in Table 1 (Figure 5a,b). The subsidence values calculated according to Equation (11) were larger over the entire range of thicknesses considered than the subsidence volumes calculated by Equations (1), (2), (4), and (7). Equations (1) and (2) were developed for low (0.4–0.6 m) and medium (0.6–1.0 m) drainage intensity conditions, respectively (Table 1). When calculated with Equations (5) and (6), the results coincide with those obtained using Equation (11) in the thickness range from 30 cm to 100 cm. Next, as the thickness increases, the calculation results of the proposed Equation (11) increase significantly with respect to the other equations. In the case of Equation (3) (high drainage intensity, i.e., 1.0–1.2 m) for thickness changes within 30–100 cm, the calculations exceed the settling calculated by Equation (11). In the range of larger thicknesses of the peat deposit (above 100 cm), the results of the subsidence calculation by Equation (3) are similar to those obtained by Equations (5) and (6) (Figure 5a).

In the case of the equations for the average annual subsidence rate (Figure 4b), the calculation of this parameter with the proposed equation (12) outperforms the results of calculations using Equations (8)–(10) (Table 1) in the prevailing range of soil thickness changes, i.e., from 70 cm to 250 cm. Equations (11) and (12) were developed based on the results of settling measurements of the drainage ditch banks, where the intensity of drainage is locally the highest and consequently causes the greatest settling of the soil surface.

**Figure 5.** *Cont*.

**Figure 5.** Comparison of the results of calculating the amount of surface subsidence using proposed Equation (11) (**a**) and the amount of average annual rate of subsidence using proposed Equation (12) (**b**) with the results of calculating these parameters with equations published in the literature.

#### **4. Discussion of Results**

The lowering of the water table as a result of dewatering causes the buoyancy force to disappear and the pressure of the dewatered layers on the lower lying layers to disappear, resulting in the subsidence of the peat deposit at depth. The land surface, which is composed of the sum of the subsidence of the individual layers and, to a lesser extent, the bottom of the trenches, decreases the most [22,25,26,28]. When conducting long-term studies (1965–1998) on subsidence of peat soils in northern Poland, it was found that drainage pipelines initially located about 1 m below the ground surface decreased their depth by about 40–50 cm due to surface and bottom subsidence. The average annual rate of decrease in their depth was contained in the range of 1.21–1.51 cm year<sup>−</sup>1. In the conditions of the northern Netherlands, according to Van den Akker et. al. [26,56], the position of the water table in the ditches on peat soils used as pastures was mostly maintained at about 60 cm below the ground surface. In these soils, there was a decrease in the water table in the ditches and their bottoms observed as a result of the subsidence process, by about 10 cm on average over a period of 10 years. Maintaining the ditch water table between 0.30 and 1.20 m below ground surface on these soils results in an average annual subsidence rate ranging from 0.50 to 2.20 cm year−<sup>1</sup> [26,56]. Grzywna [57,58], G ˛asowska [37], Oleszczuk et al. [49] observed a decrease in the designed ditch depth of 1.00 m on average to a postsubsidence depth of 0.20–0.60 m over a period of about 40–50 years on this type of site (drained peat soils, grassland use) in central and eastern Poland. In the case analyzed in this work, the design depth of the 12 ditches in 1967 ranging from 0.73 to 1.35 m (average depth of about 0.97 m) decreased over a period of 47 years to post-subsidence depths ranging from 0.40 to 1.05 m, with an average of 0.71 m.

Ostrom ˛ecki [53] and Pierzgalski [54] determined the ranges for the values of the ditch shallowing post-subsidence depth ratio (*d/t*) in relation to the degree of organic matter decomposition in peat soil. In the studied peatland, the values of this coefficient varied over a very wide range, with an average of 0.41, which indicates that most of these soils had an average degree of organic matter decomposition. This is confirmed by independent physical tests conducted under laboratory conditions on soil samples from this site [37,39,45,50,51].

The value of subsidence of the drained organic soils surface is usually expressed in cm and quite often related as a percentage of its initial thickness. However, average annual rate of subsidence expressed in cm year−<sup>1</sup> or mm year−<sup>1</sup> allows the results obtained from this process to be compared across many countries and continents in the world [27].

The average annual subsidence rate in the vicinity of drainage ditches in the study area of Solec peatland was 0.92 cm year−1. Similar studies of the subsidence rates conducted on a section of this peatland by Oleszczuk et al. [45] over a 40-year period (1978–2018) indicated an average annual subsidence rate of 0.62 cm year<sup>−</sup>1. The authors then showed that the subsidence rate also depended on the location of the measurement points in relation to the existing drainage system. Out of the 14 measurement points analyzed, four of them were located in the immediate vicinity of drainage ditches and the Mała River. The average annual rate was relatively varied across these points and amounted to 0.25 cm year<sup>−</sup>1, 0.48 cm year<sup>−</sup>1, 1.20 cm year−1, and 1.45 cm year−1, respectively [45]. It can be concluded that typically the surface settles most along localized drainage facilities (canals, ditches, drains) because the depression curve of the groundwater table near these facilities assumes the greatest depths [37]. This results in a significant decrease in the initial depth of channels and ditches, and a decrease in the depth of drain location. When peatlands are used agriculturally, this poses a threat to the continued proper functioning of drainage systems in these areas [49,59], and requires upgrading such systems by, e.g., dredging the existing ditch network [43,56,60–62]. On the other hand, on sites where agricultural use has already been withdrawn, subsidence and shallowing of this infrastructure may contribute to limiting the functioning of these systems, resulting in a positive effect of inhibiting or blocking water runoff, limiting the degradation of peat soils [62–64]. These topics have been widely reported in the literature. In this context, the authors of this paper emphasize that they focus only on the technical aspect related to the subsidence effects on drainage systems.

The amount of settling of drained peat soils and its rate depends on numerous factors: the depth of drainage of the peat deposit, its initial thickness, the botanical composition of peat and its decomposition degree, physical and chemical properties, the time elapsed since drainage, the type of use, and the climatic conditions in the area [20, 2540]. In order to develop empirical relationships expressing the magnitude/average annual subsidence rate of a drained peat soil surface of practical dimension, it is not possible to consider all these factors. Relationships presented in the literature expressing the magnitude/rate of subsidence have considered the initial thickness of the peat deposit, the length of time since dewatering, the depth/intensity of dewatering, and the position of the groundwater table [20–22,25,27,42,45,58,61,65,66]. Most commonly, simple empirical relationships between the magnitude/rate of annual subsidence and the initial thickness of the peat deposit are found in the literature (see Table 1). Some results of measurements of this type published in the literature show a large variation in the amount of subsidence of peat soils, depending on their initial thickness. For example, in case of an area of fen-peat soils from the Note´c area, Ilnicki [55] showed that with an initial peat deposit thickness of 300 cm, the amount of surface subsidence varied from about 4 cm to about 80 cm over a period of 43 to 66 years. Consequently, the correlation coefficients of this type of linear relationship are *r* = 0.56 for low drainage intensity (0.40 m–0.60 m), *r* = 0.25 for medium drainage intensity (0.60 m–1.00 m), and *r* = 0.58 for high drainage intensity (1.00 m–1.20 m). Based on the results obtained in the peatlands of the Biebrza Marshes, Krzywonos [23] developed a similar linear relationship with a correlation coefficient of *r* = 0.624.

Comparison of the magnitude of subsidence/average annual rate of subsidence over a period of 47 years against the magnitude of the initial thickness of the peat deposit shows a rather large range of field results, which demonstrates the high degree of difficulty in determining this type of relationship. Nevertheless, the correlation coefficients obtained (0.65) are consistent with results published in the literature (e.g., [55]). The analysis of the calculation results for the magnitude and subsidence rate by empirical formulas proposed by different authors indicates that Equations (11) and (12) proposed in this paper can be used to estimate the magnitude of subsidence or the average annual subsidence rate of the drained fen-peat soil surface, with a medium degree of decomposition near drainage ditches in areas with an initial deposit thickness of up to 250 cm.

#### **5. Conclusions**


**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/land10121287/s1, Table S1: Results of field measurements (2014) and archival data (1967) of the analyzed ditches on a part of the Solec peatland; Table S2: Values of initial thickness of the peat deposit (1967) and the size of the subsidence and the average annual rate of subsidence in the years 1967–2014.

**Author Contributions:** Conceptualization: R.O., E.Z., J.U.; methodology: J.U., R.O., E.Z., validation: R.O., E.Z., J.J.; formal analysis: R.O., E.Z., J.U.; investigation: R.O., J.U.; data curation: R.O., J.U., J.J.; writing—original draft preparation: R.O., E.Z., J.J.; writing review and editing: E.Z., J.J.; visualization: J.U., E.Z.; project administration: E.Z., J.J. All authors have read and agreed to the published version of the manuscript.

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

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


## *Article* **Impact of Land Use Changes on the Diversity and Conservation Status of the Vegetation of Mountain Grasslands (Polish Carpathians)**

**Jan Zarzycki 1,\*, Joanna Korzeniak <sup>2</sup> and Joanna Perzanowska <sup>2</sup>**


**Abstract:** In recent decades, in the Polish Carpathians, agriculture has undergone major changes. Our goal was to investigate whether the former management (plowing or mowing and grazing) had an impact on the current species composition, diversity and conservation status of the vegetation of grazing areas. We carried out vegetation studies on 45 grazing sites with traditional methods of grazing (transhumant pastoralism). The survey covered both old (continuous) grasslands and grasslands on former arable land. The most widespread were *Cynosurion* pastures and mesic *Arrhenatherion* grasslands. Wet *Calthion* meadows occurred at more than a half of grazing sites, while nutrient-poor *Nardetalia* grasslands were only recorded at several grazing sites. For each grazing site, we used soil maps from the 1960s to read land use in the past. We mapped present grassland and arable land area. Compared with the 1960s, there was a significant decrease in the area of arable land and an increase in grasslands. Species diversity was greater in grazing sites where grasslands developed on former arable land. However, this diversity was associated mainly with the occurrence of common grassland species. *Cynosurion* pastures and wet *Calthion* meadows had the best conservation status, while nutrient-poor *Nardetalia* grasslands were the worst preserved. We concluded that the conservation status of mesic grasslands and pastures is dependent on the present diversity of land use within a grazing site, rather than the land use history 60 years ago. This is the first study of the natural, not economic, value of pasture vegetation in the Polish part of the Carpathians.

**Keywords:** grazing management; biodiversity; high-nature-value farming; old field grassland

#### **1. Introduction**

Grasslands provide a variety of ecosystem services, such as cultural landscape and environmental values [1]. They are significant reservoirs of biodiversity and belong to the most species-rich plant communities in Europe [2]. Apart from an abundance of plant species, they are also the habitat for a multitude of animal species from different taxonomic groups, such as insects, arachnids, and snails [3]. The biodiversity of grassland is also related to the time since the transformation of arable land into grassland.

The majority of grasslands in Europe originated as a result of the interaction of environmental conditions and farming by grazing and mowing [4]. However, social and economic transformations in the 20th century provoked changes in land use and management practices in agriculture. In Western Europe, this process particularly accelerated in the 1950s/1960s [5], while in Eastern Europe after the political transformation at the beginning of the 1990s [6]. As a result, the area covered by grasslands has been decreasing. On the one hand, wherever land is suitable for agricultural use, grasslands are transformed into arable fields or the intensity of such use increases [7,8], which results in the development of highly productive but species-poor communities [9]. On the other hand, in marginal

**Citation:** Zarzycki, J.; Korzeniak, J.; Perzanowska, J. Impact of Land Use Changes on the Diversity and Conservation Status of the Vegetation of Mountain Grasslands (Polish Carpathians). *Land* **2022**, *11*, 252. https://doi.org/10.3390/ land11020252

Academic Editor: Alexandru-Ionu¸t Petri¸sor

Received: 13 December 2021 Accepted: 1 February 2022 Published: 8 February 2022

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

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

areas less suitable for agriculture, land use is abandoned, which promotes secondary forest succession and the disappearance of grassland communities [10–12]. It is particularly disadvantageous because grassland communities in infertile habitats are characterized by the highest natural value and diversity [9,13]. In the Carpathians, the expansion of forest cover has been progressing for over a century and has especially hastened in recent decades [14]. Despite a significant increase in forest cover, the area of grasslands is not only not reduced but even grows [15], as former arable fields are transformed into meadows and pastures. Similar processes were also observed in other montane regions, e.g., in the Caucasus [16], in Slovakia [17] and in the Apennine mountains [18]. As a consequence, a large portion of grasslands in the Carpathians is situated at present on the former arable fields. They are usually located at lower altitudes and on more fertile soils. The species composition of such communities also formed over a short period. For this reason, these plant communities are usually species-poor and are characterized by a lack of many specialist grassland species [19,20].

A number of studies indicate that the species composition of communities, their diversity and the occurrence of rare species are related not only to the current management but also to the habitat history [21,22]. The historical land use and land use sequences shaped the vegetation of Swedish semi-natural grassland more so than the current management. The highest diversities of grassland plants have been found in pastures continuously grazed since the 18th century [21]. The continuity of extensive mowing and grazing has a positive effect on the occurrence of grassland and forest edge specialist species [22]. The spatial context is also important: the present and historical habitat connectivity influences fine-scale plant species diversity in grazed temperate semi-natural grasslands [23].

Traditional pastoral systems in many mountain areas (Carpathians, Alps, Pyrenees) consisted of transhumance, i.e., grazing livestock at higher altitudes in summer and descending with the herd into valleys for the winter [24]. Such a pastoral system has been used up till now in Romania [25], Ukraine [26], Switzerland, Scandinavia, France and Spain [27]. In Poland, at the end of the 1980s, the period of transition into an open-market economy saw a livestock population crash; for instance, the sheep population was dramatically reduced from 4,837,000 in 1985 to 266,911 in 2018 [28], due to the decreased profitability of agricultural production. However, since the mid-2000s, we have witnessed a revival of transhumant pastoralism in the Polish Carpathians [29] thanks to subsidies under the Common Agricultural Policy (CAP). However, the most common locations of transhumance have changed [29]. New grazing sites are located at lower altitudes and primarily on the site of arable fields (Figure 1). The term grazing sites suggests that they are covered only by pastures. However, in practice, parcels of land used as grazing sites are very diverse and harbor natural and semi-natural types of vegetation, and also arable fields. A sheep herd migrates daily within the grazing sites and in the evening the sheep return to a portable shelter. Biodiversity conservation in rural areas is analyzed in connection with viewing these areas as a patch–corridor matrix, which is a dynamic system linking natural processes with current human impact and its future consequences [30]. This approach, based on a landscape-level perspective, is increasingly often adopted in the conservation of semi-natural grasslands [31,32].

The aim of this work was to characterize the sites used for traditional pastoralism, treated as a kind of functional unit, in terms of (a) the species composition and diversity of plant communities occurring in these areas, (b) the evaluation conservation status of these communities, and (c) the impact of land use in the past and historical and present grassland area to arable land area ratio, on the species diversity and conservation status of grassland.

**Figure 1.** Grazing on former arable land in the Pieniny Mts.

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

#### *2.1. Study Area*

The study was carried out at 45 grazing sites located in 7 physico-geographical regions in the Polish part of the Western Carpathians (Figure 2) extending 45 km N–S and 95 km E–W. The extent of the area from the west to the east in combination with a wide range of altitudes (347–1159 m a.s.l.) caused large variations in climatic conditions. The area of grazing sites was diverse and ranged from 15 to 145 hectares (Table 1). All study sites have been extensively grazed by sheep (livestock density did not exceed 0.5 LU/ha) for at least a few years, but most of them have been used in this way for a long time. The survey covered both old (continuous) grasslands and those on former arable land. Vegetation of the grazing sites primarily comprised different types of grasslands with scattered patches of forest, thickets, fields and fallow lands, wetlands covered by fen and marsh vegetation, and sporadic small orchards.

#### *2.2. Methods*

Vegetation studies were performed in summer 2017. To analyze the species composition of the main types of grasslands, phytosociological relevés were taken according to the Braun-Blanquet method in 25 m<sup>2</sup> plots [33]. Since the grazing sites differed in area, we took 1 relevé per 10 ha of grazing site area. In total, in all 45 sites, 420 relevés were performed. The thickness of the litter layer was also measured for each relevé. In order to assess the conservation status of the grassland types within each grazing site, their structure and function were evaluated (a proxy of habitat quality). The assessment was based on the methodology applied in the Monitoring of Natural Habitats in Natura 2000 areas, which is conducted in Poland by the Chief Inspectorate for Environmental Protection [34]. Since many patches of mesic grasslands were of transitional nature between habitat types 6510 (lowland hay meadows) and 6520 (mountain hay meadows), they were analyzed jointly as mesic *Arrhenatherion* grassland. Nutrient-poor *Nardetalia* grasslands were classified as the Natura 2000 habitat 6230. The assessment of the structure and function of *Cynosurion* pastures and wet *Calthion* meadows, not included in Annex I to the Habitats Directive [35], was based on an analogical formula, as for habitats 6230, 6510 and 6520. Habitat structure and function were assessed according to a three-point scale: favorable, unfavorable

inadequate and unfavorable bad. The composite result of the evaluation of indices was defined individually for each habitat. The following indices were assessed: spatial structure of habitat patches, diagnostic species, invasive alien species, expansive species of herbs, and expansion of shrubs and trees [34]. The nomenclature of vascular plant species was adopted according to Mirek et al. [36]. Field work also included mapping the grassland and arable land area.

The variability of the species composition of grassland vegetation was assessed based on phytosociological relevés. To reduce redundancy due to the oversampling of some areas, the set of 420 relevés was subjected to stratified resampling [37] in JUICE 7.0 [38]. We assumed that the mean geographical distance between all pairs of relevés belonging to the same grassland type should be at least 0.5 km, which corresponds to a geographical grid of 0.25 latitude × 0.4 longitude. The resulting data set used in ordination analyses contained 294 relevés. The unweighted mean Ellenberg Indicator Values (EIV) were calculated for every relevé.

Agricultural soil maps from the 1960s, digitized and georeferenced using QGIS 3.16.7., were used to calculate the area of grassland and arable land at that time in each grazing site. As a quality assessment of grazing sites for agricultural production, indicators of the quality and usefulness of agricultural soils (QUAS) were used [39]. These indicators are the result of calculations made for soil valuation classes and soil-agricultural complexes shown on soil maps. The method of land use (arable vs. grassland) where the relevé was located was also read from the soil maps.

**Figure 2.** Locations of the grazing sites within regions in the Polish Carpathians. The regions are displayed as colored circles. Explanations: (**a**) rivers and water bodies, (**b**) forests, (**c**) localities, (**d**) state border.


**Table 1.** Diversity of grassland types in relation to land use in the past. Land use in the 1960s: G-grassland; A-arable fields.

The diagnostic species [40] for the grassland phytosociological syntaxa (*Calthion*, *Arrhenatherion*, *Polygono-Trisetion*, *Cynosurion*, *Nardetalia*, *Arrhenatheretalia*, *Molinio-Arrhenathereatea*) were distinguished and used to classify relevés into the main types of grassland vegetation. The main diversity measures for sites were determined: average number of vascular plant species in a relevé (α diversity), total number of vascular plant species in all relevés within a grazing site (γ diversity) and the number of diagnostic species for grassland syntaxa.

To explore major gradients in species composition of the main grassland types and their relationship with environmental characteristics and past land use, the Detrended Correspondence Analysis (DCA) was performed using the CANOCO 5.10. For the ordination analysis, species abundance data from 294 plots were square-root transformed with the downweighting of rare species [41]. Because of a lack of normality, the Spearman rank-correlation coefficient for the analyses of relationships between site species diversity, habitat quality and site environmental characteristics was used. Correlations were analyzed using Statistica 13.1.

#### **3. Results**

#### *3.1. Changes in Land Use*

In the 1960s, arable fields dominated in the areas of analyzed grazing sites, and only 3 sites were entirely without them (Figure 3). In 8 sites, 90% of the area was covered by arable fields. In 2017, the area of arable land was marginal and 22 grazing sites were completely devoid of them, while in 16, the share of arable fields to the total grazing site area was lower than 10%. Only in two grazing sites was the share of arable fields relatively high, amounting to 40%. A QUAS calculated for grazing sites were low and as many as 39 grazing sites did not exceed 35 points (on an 18–100-point scale) (Figure 4).

#### *3.2. Species Composition and Diversity of Four Types of Grasslands*

The most widespread vegetation types were pastures of the *Cynosurion* R. Tx. 1947 alliance, mostly *Festuco-Cynosuretum* Büker 1941, and mown and grazed grasslands of the *Arrhenatherion elatioris* (Br.-Bl. 1925) Koch 1926 alliance. Due to a diverse pattern of land use as hay meadows and pastures, many phytocenoses were of a transitional nature, between a pasture and a hay meadow. Wet eutrophic meadows of the *Calthion palustris* R. Tx. 1936 em. Oberd. 1957 alliance occurred at more than a half of grazing sites. *Cirsietum rivularis* Nowi ´nski 1927 was encountered most often, *Scirpetum sylvatici* Ralski 1931 was rarer while *Epilobio-Juncetum effusi* Oberd. 1957 occurred at intensely grazed places. Nutrient-poor grasslands with *Nardus stricta* of the *Nardetalia* Prsg. 1949 order were only recorded at 18 grazing sites. They usually appeared as small patches, often with a large proportion of species typical of the neighboring meadows.

**Figure 3.** Number of grazing sites in classes of arable field area in the 1960s and in 2017.

**Figure 4.** Number of grazing sites in classes of quality and usefulness of agricultural soils index (QUAS).

The main grassland vegetation variability gradient represented by the first DCA axis is associated with soil moisture (Figure 5). The second DCA axis mostly reflects differences in soil fertility, separating nutrient-poor *Nardetalia* grasslands from the remaining grassland types. The occurrence of nutrient-poor grasslands is positively correlated with slope and litter thickness, which suggests a lack of agricultural use. Wet *Calthion* meadows are associated with "old" grasslands. For the remaining types of grasslands, historical land use either as arable fields or grassland is not important. The diagram shows a large degree of similarity between the species composition of *Cynosurion* pastures and mesic *Arrhenatherion* grasslands, and dissimilarity with wet meadows and, to a lesser extent, also *Nardetalia* grasslands. Wet *Caltion* meadows were observed to include a lower proportion of the relevés located on former arable fields compared with the remaining grassland types (Table 2). On the other hand, relevés classified as mesic *Arrhenatherion* grasslands occur more often on former arable fields than on "old" grasslands. In pastures and nutrient-poor *Nardetalia* grasslands, the number of relevés taken on former fields and "old" grasslands is almost equal.

**Figure 5.** The Detrended Correspondence Analysis (DCA) ordination diagram of 294 plots categorized by grassland type with environmental vectors as overlay. Eigenvalues were 0.506 and 0.276 for the first and second axis, respectively. Cumulative percentage variance of species data for the first and the second axis were 4.9% and 7.5%, respectively. Ellenberg Indicator Values: F-moisture, L-light, N-nutrients.


**Table 2.** Diversity of grassland types in relation to land use in the past. Land use in the 1960s: G-grassland; A-arable fields.

The total number of species (gamma diversity) in grasslands situated in grazing sites ranged from 49 to 104, but only three sites were very poor in species (below 60). The average number of species per phytosociological relevé (alpha diversity) was 26.3, ranging from 17.6 to 34.3.

#### *3.3. Conservation Status*

The conservation status assessment of grasslands (Figure 6) showed the worst status of nutrient-poor *Nardetalia* grasslands, which were also rare in grazing sites and occupied the smallest surface areas. The conservation status of the wet meadows was assessed as favorable in one third of grazing areas. Mesic grasslands were present in almost all grazing areas; in more than half of grazing sites, their conservation status was evaluated as inadequate, in nine areas as bad and in eight areas as favorable. *Cynosurion* pastures were recorded in all grazing areas and their status was favorable in almost half, and was assessed as bad only at five sites.

**Figure 6.** Conservation status of four grassland types evaluated on the basis of an assessment of their structure and function.

#### *3.4. The Effect of Land Use on Species Diversity and Conservation Status of Grasslands*

A QUAS index was associated with a greater total number of species in a grazing site and a higher number of species diagnostic of the *Cynosurion*, *Nardetalia*, and other grassland species, but with a lower number of species diagnostic of the *Calthion*. The number of species diagnostic of the *Cynosurion* and species diagnostic for the *Molinio-Arrhenatheretea* class and for the *Arrhentheretalia* order (*M*-*A*+*A*) was positively correlated with the share of arable fields in land use structure in the 1960s, while the number of species diagnostic of the *Calthion* was negatively correlated with this share. The impact of the current share of arable field area in land use structure was much greater and showed a positive correlation with the total number of species in a grazing site, the average number of species in a relevé, and the number of species diagnostic of the *Cynosurion*, grassland species and all other species, but had no effect on the number of species diagnostic of wet meadows (Table 3).

**Table 3.** Spearman rank-correlation coefficients between grazing site diversity metrics, habitat quality and grazing site characteristics. Abbreviations: *M*-*A*+*A*-species diagnostic for the *Molinio-Arrhenatheretea* class and for the *Arrhentheretalia* order.


\* significant *p* < 0.05.

#### **4. Discussion**

#### *4.1. Species Composition and Diversity in Relation to Land Use Change*

The results of the research indicate a process of increasing similarity in the species composition of grassland in areas where transhumant pastoralism is used. This applies in particular to communities of the alliances *Arrhenatherion* and *Cynosurion* because they can occur in the same environmental conditions, while the management type is a differentiating factor. Grazing caused the species composition of different communities to become similar. This is not the case for wet *Caltion* meadows and nutrient-poor *Nardetalia* grassland. The *Calthion* meadows occur in wet places and usually are only sporadically grazed. They are most often found in grazing sites with poor-quality soils at higher altitudes and with a slight proportion of arable fields in the past. Nutrient-poor *Nardetalia* grasslands in the Polish Carpathians are rare [42]. They most often cover small patches in habitat conditions differing from the surrounding grasslands. Thus, they can appear as enclaves in grazing areas among generally more fertile soils.

The highest biodiversity, as measured by the average number of all vascular species, pasture species and, generally, grassland species, was characteristic of grazing sites comprising a large proportion of arable fields. A significant proportion of arable fields in grazing sites in the past and their remains at present resulted in a high heterogeneity of habitats, related to the occurrence of such elements as field boundaries, clearance cairns, and dirt roads with roadsides. They are overgrown by plant species which increase species diversity in the landscape [32] and can easily migrate to grasslands. The heterogeneity of habitats and the related total number of species does not increase the abundance of grassland specialist species, since the latter occur in old grasslands [43,44]. Additionally, in

our studies, in contrast to grassland generalist species (*M-A+A*), no correlation was found between the number of grassland specialist species and the proportion of arable fields in the grazing site.

The large number of pasture species in a grazing site was correlated with a higher soil quality index and with a higher share of arable fields both at present and in the past. It indicates that ex-arable fields create favorable conditions for pasture communities. Grazing enables the dispersion of species by endo- and epizoochory [45], while disturbances caused by trampling and browsing by animals increase the probability of successful recruitments [46]. Similar relationships were observed in the case of grassland generalist species. However, such use may be insufficient for grassland specialist species [47].

The relationship between present species composition and land use type in the 1960s is not strong. The land use data for that period allow only for the conclusion that grasslands on ex-arable fields do not exist for longer than 60 years. However, this period can be much shorter. Transformations occur gradually; thus, grasslands differ in age, which results in a variable species composition because the migration of propagules depends on time [47]. Öster et al. [48] indicated that over a 50-year period, only 50% of species occurring on neighboring semi-natural grasslands migrated to ex-arable fields. Additionally, investigations on dry grasslands demonstrated a significant role of former land use in terms of species distribution [49].

#### *4.2. Conservation Status and Transhumant Pastoralism as a Protective Measure*

Conservation status expressed by structure and function assessment index for all grassland types was not statistically significantly correlated with soil suitability for agriculture or present or past land-use type. Nevertheless, these are dynamic ecosystems in which land use-dependent vegetation changes are relatively fast. In another study [50], significant vegetation changes were found in abandoned mountain meadows only 6 years after the reinstatement of grazing, and after 9 years the proportion of species typical of floristically rich grasslands was increased. The interval of time analyzed by us (ca. 60 years) is probably too long to sustain a statistically significant impact of land use type at that time.

Evaluation of the restoration efficacy of semi-natural grasslands on former croplands should take into account the similarity of the species composition of created communities to the species composition of communities typical of local habitat conditions (referenced community), because it is a better indicator of success than species diversity (number of species) [51]. If no local reference communities have been preserved or the study areas are widely diverse, as in our case, it is necessary to strive for the achievement of conservation status, described as favorable for grassland communities, following the criteria formulated by State Environmental Monitoring [52–54].

The transformation of former arable fields into permanent grasslands is not currently a common practice in Europe, which is related to the extensification of land use, e.g., in protected areas, especially in the mountains [7], or the implementation of agri-environmental programs [55]. The restoration of multispecies communities on formerly intensively used arable fields is difficult and long lasting, and requires the application of diverse measures [56,57]. However, under favorable conditions, this process can be much faster [58]. In the Polish Carpathians, this process proceeds spontaneously through the migration of grassland species from neighboring areas to ex-arable fields. A majority of grazing sites are characterized by a mosaic of small patches of different land use types and by diverse vegetation; thus, an easily accessible source of propagules is available, which is crucial for grassland restoration [57,58]. High doses of fertilizer have never been used in the Polish Carpathians, while high-fertility soils are a significant obstacle to the restoration of species-rich grasslands [56]. A similar mechanism of preservation of high diversity was suggested by Janišová et al. [59] based on studies of semi-natural grasslands in Slovakia.

Traditional transhumance is practiced in the Carpathians, which involves livestock grazing in summer and the migration of herds within a large area covered by a mosaic of semi-natural grasslands, arable fields, ex-arable fields, trees and other small landscape features, allows species diversity to be preserved [59,60] because it supports natural processes, according to the metacommunity theory [61]. Extensive grazing can also produce negative effects from a nature conservation perspective if it is used in communities created by different management types. Grazing is considered to be more advantageous for species diversity than mowing [62]. However, many types of grasslands with high natural value, such as the mesic grasslands occurring in the study area, have developed and currently exist due to a proper mowing regime [63]. Plant diversity and conservation status largely depend on diversified land use, as revealed by the studies of Kun et al. [64] carried out in the Romanian Carpathians and Tölgyesi et al. [65] in Hungary. Hence, the preservation of the mosaic spatial structure of natural and semi-natural features and extensive farming practices in connection with diverse forms of human impact, such as grazing and mowing, is the optimal solution to preserve diversity at the levels of species, community, and landscape in mountain areas.

#### **5. Conclusions**

The transformation of extensively used arable fields into grasslands is relatively fast, but the species composition of communities formed in this way is dominated by less specialized species.

The biodiversity of analyzed grazing sites is dependent on the recent grassland area to arable land area ratio, rather than on the 60-year land use history. Nevertheless, detailed knowledge of past land use should be an integral part of analyses of the current state of vegetation and its dynamics. The spatial context (landscape connectivity, habitat fragmentation) is also worth considering.

Traditional transhumance can have a beneficial effect on the species diversity of grassland communities that have developed on former arable fields. Extensive open grazing over vast areas provides a real opportunity to maintain this type of ecosystem in good condition and should be subsidized by the state.

**Author Contributions:** Conceptualization, J.Z., J.K. and J.P.; methodology, J.Z., J.K. and J.P.; software, J.Z.; investigation, J.Z., J.K. and J.P.; writing—original draft preparation, J.Z., J.K. and J.P.; writing review and editing, J.Z. and J.K.; visualization, J.K. and J.P.; project administration, J.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This manuscript is not funded by a specific project grant.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The study was conducted in the framework of the statutory activity of the Institute of Nature Conservation, Polish Academy of Sciences and with a subsidy of the Ministry of Science and Higher Education for the University of Agriculture in Kraków in 2022. The results from field work carried out for the project "Maintenance of biological diversity of mountain meadows and pastures through pastoral management" were used. The project was realized by the Landscape Parks of the Malopolska Voivodship, co-financed by the European Regional Development Fund (ERDF) under the Regional Operational Program of the Malopolska Voivodship for the years 2014–2020.

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

#### **References**


## *Article* **Influence of Anthropogenic Load in River Basins on River Water Status: A Case Study in Lithuania**

**Laima Cesonien ˇ e˙ 1,\*, Daiva Šileikiene˙ <sup>1</sup> and Midona Dapkiene˙ <sup>2</sup>**


**\*** Correspondence: laima.cesoniene1@vdu.lt; Tel.: +370-37-752224

**Abstract:** Twenty-four rivers in different parts of Lithuania were selected for the study. The aim of the research was to evaluate the impact of anthropogenic load on the ecological status of rivers. Anthropogenic loads were assessed according to the pollution sources in individual river catchment basins. The total nitrogen (TN) values did not correspond to the "good" and "very good" ecological status classes in 51% of the tested water bodies; 19% had a "bad" to "moderate" BOD7, 50% had "bad" to "moderate" NH4-N, 37% had "bad" to "moderate" NO3-N, and 4% had "bad" to "moderate" PO4-P. The total phosphorus (TP) values did not correspond to the "good" and "very good" ecological status classes in 4% of the tested water bodies. The largest amounts of pollution in river basins were generated from the following sources: transit pollution, with 87,599 t/year of total nitrogen and 5020 t/year of total phosphorus; agricultural pollution, with 56,031 t/year of total nitrogen and 2474 t/year of total phosphorus. The highest total nitrogen load in river basins per year, on average, was from transit pollution, accounting for 53.89%, and agricultural pollution, accounting for 34.47%. The highest total phosphorus load was also from transit pollution, totaling 58.78%, and agricultural pollution, totaling 28.97%. Multiple regression analysis showed the agricultural activity had the biggest negative influence on the ecological status of rivers according to all studied indicators.

**Keywords:** pollution; ecological status indicators; water quality

#### **1. Introduction**

Lithuania is committed to achieving the objectives of the EU Water Framework Directive by 2027 and to achieving good water status in inland waters. There are approximately 30 thousand rivers and creeks in Lithuania, with a longer than 200 m, reaching an overall sum of 63,700 km. Although the ecological condition of Lithuanian rivers has been mostly improving over the last few years, it has been determined that only 49% of them correspond to a good ecological state [1], in the period of 2010–2013.

Human activities change the hydromorphological, physicochemical, and biological parameters of surface water bodies, affecting their biodiversity and ecological functioning [2–9]. Pollution caused by anthropogenic activity sufficiently worsens the state of water ecosystems, resulting in water quality degradation, reducing the potential for water use in various fields and endangering people health [10–14].

Sources of pollution are divided as follows: background pollution (forests), diffuse (non-point) pollution from agricultural lands, surface sewage that is not treated in wastewater treatment plants (WWTPs), and concentrated (point) pollution caused by households, urban, municipal, industrial wastewater (wastewater treatment plants), and others. Human activities affect the condition of water bodies differently in rural areas (agricultural activities, livestock) and urban areas (industrial, municipal, domestic wastewater discharges). Changes of land use also have a negative effect on water condition of rivers [15–21]. Landscape changes caused by anthropogenic activities and land cover make a significant

**Citation:** Cesonien ˇ e, L.; ˙ Šileikiene, D.; Dapkien ˙ e, M. Influence ˙ of Anthropogenic Load in River Basins on River Water Status: A Case Study in Lithuania. *Land* **2021**, *10*, 1312. https://doi.org/10.3390/ land10121312

Academic Editors: Wiktor Halecki, Dawid Bedla, Marek Ryczek and Artur Radecki-Pawlik

Received: 25 October 2021 Accepted: 26 November 2021 Published: 28 November 2021

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

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

influence on the state of surface waters [22–24], are closely related to water chemical parameters [25], the diversity of fish and macroinvertebrate species [24,26], and sediment metal concentrations [27]. Fragmented urban land use, with many impermeable surfaces, tends to increase river flows and adversely affects water quality [28–30].

Many research have been conducted to evaluate the ecological status of surface waters [31–36]. The multiple anthropogenic pressures (pollution, hydrological, and hydromorphological alterations) on the ecological status of the European rivers were assessed. In one third of the territory of the European Union, the ecological status of rivers has been found to be good, which is linked to the existence of natural areas, and urbanization is leading to poorer ecological status [37].

A concentrated source of pollution is wastewater from factories and households, and it is insufficiently treated and controlled [38–40]. Joshua et al. (2017) have pointed out that substandard practices of wastewater treatment are utilized in developing countries. In these countries, the major causes for this include the insufficient number of wastewater treatment plants, overloading and ineffective operation of existing WWTPs, and others [41].

Diffuse pollution is one of the main problems in irrigated areas. Intensive irrigation and fertilization of arable land and pastures have the greatest negative impact [42,43].

The sewage from residents and their house-holds' lack of connection to the wastewater collection networks also forms diffuse pollution that enters rivers [44]. Diffuse pollution covers large areas, and all such areas are polluted quite equally [45]. Very important diffuse pollution source affecting the condition of rivers is inadequate farming [46–49]. The riverbanks are a suitable place for agriculture, as they are particularly productive. Diffuse pollution is caused not only by farming on riverbanks, but also by agricultural activities in all the river basin. The use of chemical fertilizers and pesticides is an essential part of diffuse pollution entering surface water bodies [50]. Both organic and mineral fertilizers are used to fertilize crops; however, due to erosion, soil leaching, and runoff, they enter surface water bodies and contaminate them with nutrients [51]. Agriculture was the cause of 48% of the deterioration in water quality of surface water bodies of USA [52]. Approximately 30–35% of nitrogen and 10–15% of phosphorus, which pollute surface waters, have been found to come from agricultural activities [53].

Generally, Lithuanian surface water bodies are impacted by both diffuse and concentrated source pollution. In Lithuania, the effect of pollution on the condition of Venta and Muša-Lielup ¯ e river basins was evaluated [ ˙ 38,54,55]. According to the biochemical oxygen demand (BOD7), human activity accounted for 56% of the influence of pollution on water quality, and 90% of the annual borne total nitrogen (TN) and 78% of the total phosphorus (TP) content in the Merkys River [56]. Of the bodies of water, 46% did not achieve a "good" status in terms of nitrate nitrogen in 2017 [57]. Export coefficients and the retention of biogenic nutrients in Lithuanian river basins were assessed using the MESAW statistical model. The export coefficients of TN and TP showed much higher values from arable land in comparison to forest area, pastures, and meadows from the studied Merkys, Muša, ¯ Žeimena, and Nevežis river basins, and retained from 67% to 78% of the total nitrogen and ˙ from 24% to 63% of the total phosphorus [58].

We performed a study of the effectiveness of diffuse pollution abatement measures in reducing the nutrient pollution of surface water bodies in Lithuania in the context of climate change. The SWAT model was used for this purpose. The results show that climate change is a significant factor in changing the effectiveness of measures to reduce diffuse pollution. The most effective measures to reduce nutrients inputs to water bodies were identified, including pasture/meadow expansion, stubble abandonment for winter, and catch crop cultivation; arable farming was the least efficient method [59].

The aim of this research is to evaluate the impact of anthropogenic load on indicators of the ecological status of rivers.

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

#### *2.1. Study Area*

To determine the risk of water bodies that do not comply with the water quality standards, the physicochemical quality indicators at 94 locations of 24 rivers were studied. Water samples were taken between January and March, April and June, July and September, and October and December in 2014–2020. The investigated river's water sampling areas and their hydrological data are presented in Figure 1 and Table 1.

**Figure 1.** River's water sampling areas.


**Table 1.** Hydrological data of rivers.

#### *2.2. Water Quality Assessment Standards*

The ecological statuses of water bodies at risk were assessed in accordance with the Procedure for Rating the Ecological Status of Surface Water Basin [60]. Research on the physical and chemical quality of elemental indicators have been identified in a laboratory at Vytautas Magnus University. Water samples were collected according to the EN ISO requirement standards: LST EN ISO 5667-14:2016—Water Quality—Sampling—Part 14. The total nitrogen (TN) was tested according to the method of LST EN 13342—2002. The total nitrogen (TN) following determination of nitrogen—determination of bound nitrogen (TNb), following oxidation to nitrogen oxides LST EN 12260:2004. Total phosphorus (TP) analyses were assessed according to LST EN ISO 6878:2004; BOD7—according to ISO 5815-1:2003 ammonium nitrogen (NH4 +-N) LST ISO 7150-1:1998, LST EN ISO 13395:2000 nitrate nitrogen (NO3-N), and LST EN ISO 6878:2004 phosphate phosphorus (PO4-P).

#### *2.3. Presentation of Pollution Sources*

The assessment of contamination sources considered the nature of land use, the nature of cities and settlements, the location of potential sources of concentrated source pollution, the nature and intensity of economic activities in the basin and their potential impact on water bodies, recreational activities, and other economic activities that may not be in good condition according to condition requirements, and so forth.

Diffuse agricultural pollution, consisting of manure and mineral fertilizer loads resulting from agricultural activities and from the load on the population whose households are not connected to sewage collection systems.

The main sources of concentrated pollution are wastewater from cities, settlements, industrial enterprises and rain and surface water, wastewater from urban areas.

High potential for concentrated pollution to enter water bodies directly or through river tributaries.

For the quantification of pollution indicators, the following factors have been assessed:


#### *2.4. Statistical Analyses*

To statistically assess the significant impacts on quality factors related to the ecological status of water bodies, the impacts of anthropogenic load indicators TP and TN from municipal wastewater, surface wastewater, households not connected to sewage networks, agricultural land, background, and transit pollution (t/year); agricultural land, forests, wetlands, meadows, arable, infertile land, and green land (ha) on water quality indicators (Y) for the water in rivers were determined. A multiple linear regression model was applied:

$$\mathbf{Y} = \mathbf{a} + \mathbf{b}\_1 \mathbf{x}\_1 + \mathbf{b}\_2 \mathbf{x}\_2 + \dots \quad + \mathbf{b}\_k \mathbf{x}\_k. \tag{1}$$

The coefficient bj shows how much the value of Y increases (or decreases) by one unit, as xj increases when the remaining xk are fixed. *t* is Student's criterion, according to which we determined whether the bj coefficients differed statistically significantly from zero, and according to this, we decided whether the predicted values depended upon xj. The standardized coefficient beta was used to determine the relative influence of independent variables on the predicted Y. In absolute terms, a higher beta coefficient indicates greater dependence of Y on xj.

The regression model is appropriate due to the following:


#### **3. Results**

*3.1. Ecological Status Classes of the Stretches of Rivers According to the Physicochemical Values of Elemental Indicators*

Studies on the physicochemical quality of element indicators were performed for NO3-N (mg/L), NH4-N (mg/L), TN (mg/L), PO4-P (mg/L), TP (mg/L), and the BOD7 (mg/L). The results are shown in Figure 2.

**Figure 2.** Ecological status classes of rivers according to the values of indicators of the physicochemical quality of elements%.

The results presented in Figure 2 show that, according to the TN, 51% of the studied rivers did not meet the requirements of the "good" ecological class; 19% of the rivers had a "bad" to "moderate" BOD7, 50% had ''bad" to ''moderate" NH4-N, 37% had ''bad" to ''moderate" NO3-N, 4% had ''bad" to ''moderate" TP, and 4% had ''bad" to ''moderate" PO4-P.

#### *3.2. Assessment of Nutrient Loads in River Basins*

Nutrient loads in the river basins were calculated by collected the SWAT model data. Calculations were performed in tons per year for the inflows into the rivers for the total nitrogen and total phosphorus. The TN and TP loads in river basins (t/year) are presented in Figure 3.

**Figure 3.** *Cont*.

**Figure 3.** Amounts of total nitrogen and total phosphorus load in river basins (t/year).

The river basins get the largest amounts of pollution from transit loads located above the research locations, with the total nitrogen equaling 87,599.32 t/year and total phosphorus amounting to 5019.51 t/year. From agricultural activities, the total nitrogen reached 56,030.53 t/year and the total phosphorus was 2474.14 t/year. The amounts from background pollution (urban areas and forests) were 17,941.31 t/year of total nitrogen and 619.09 t/year of phosphorus. The amounts from municipal sewage were 368.66 t/year of total nitrogen and 342.12 t/year of phosphorus. The amounts from surface sewage were 321.64 t/year of total nitrogen and 29.02 t/year of phosphorus. Residents whose sewage was not discharged into sewage treatment systems generated 281.64 t/year of total nitrogen and 57.06 t/year of phosphorus.

Figure 4 shows the percentage distribution of total nitrogen and total phosphorus loads in the studied river basins.

The highest annual total nitrogen load for river basins per year, on average, came from transit pollution, accounting for 53.89%. A total of 34.47% came from agricultural pollution, 11.04% came from background pollution (urban areas and forests), 0.17% came from pollution from residents who were not connected to sewage systems, 0.20% came from surface sewage, and 0.23% came from municipal wastewater.

The highest annual load of total phosphorus in river basins was from transit pollution, accounting for 58.78%. A total of 28.97% came from agricultural pollution, 7.25% came from background pollution (urban areas and forests), 0.67% came from pollution from inhabitants who were not connected to sewage systems, 0.34% came from surface sewage, and 4.01% came from municipal wastewater.

**Figure 4.** Percentage distribution of total nitrogen and total phosphorus loads in studied river basins (names of rivers and water body code).

*3.3. Influence of Anthropogenic Loading on Total Nitrogen, Ammonium Nitrogen, Nitrate Nitrogen, Total Phosphorus, and Phosphate Phosphorus*

The influence of anthropogenic loading on total nitrogen concentration (TN is dependent variable Y) was calculated by multiple regression analysis and results are presented in Table 2.

Multiple regression analysis of the influence of anthropogenic loads on the total nitrogen concentration in the water showed that the total nitrogen value was affected by N from agricultural land, and the total nitrogen amount was generated from agricultural land and arable land (*p* < 0.05). The higher the concentrations of TN were from arable land and agricultural land, the higher the value of TN was in the water (positive function).

The effect of anthropogenic loads on the ammonium nitrogen concentration (NH4-N is dependent variable Y) was calculated by multiple regression analysis. The results are presented in Table 3.


**Table 2.** The influence of anthropogenic loads in rivers basins on the total nitrogen concentration in the water.

Dependent variable: TN; \* significance factor, *p* < 0.05.

**Table 3.** The influence of anthropogenic loads in river basins on the ammonium ion concentration in the water.


Dependent variable: NH4-N; \* significance factor, *p* < 0.05.

Multiple regression analysis of the influence of anthropogenic loads on the ammonium nitrogen concentration in the water showed that the total NH4-N value was affected by the TN from households not connected to sewage networks, the TN from agricultural land and transit pollution, and the NH4-N amount generated from forests, wetlands, meadows, arable land, infertile land, and green land (*p* < 0.05). The higher the concentrations of NH4-N were from households not connected to the sewage networks, agricultural land, transit pollution, wetlands, meadows, and arable land, the higher the value of the NH4-N was in the water (positive function). The higher the NH4-N concentrations were from forests, infertile land, and green land, the lower the NH4-N concentration was in the water (negative function).

The effect of anthropogenic loads on the nitrate nitrogen concentration (NO3-N is dependent variable Y) was calculated by multiple regression analysis. The results are presented in Table 4.


**Table 4.** The influence of anthropogenic loads in basins on the total nitrogen concentration in the water.

Dependent variable: NO3-N; \* significance factor, *p* < 0.05.

Multiple regression analysis of the influence of anthropogenic loads on the nitrate nitrogen concentration in the water showed that the NO3-N value was affected only by arable land (*p* < 0.05). The higher the concentration of NO3-N was from arable land, the higher the value of NO3-N was in the water (positive function).

The effect of anthropogenic loads on the total phosphorus concentration (TP is dependent variable Y) was calculated by multiple regression analysis. The results are presented in Table 5.

Multiple regression analysis of the influence of anthropogenic loads in basins on the concentration of total phosphorus in the water showed that the total phosphorus value was influenced by the discharge of surface wastewater from households not connected to sewage networks, agricultural land, arable land, infertile land, and green land (*p* < 0.05). The higher the TP concentration was in the surface wastewater from households not connected to sewage networks, agricultural land, arable land, the higher the TP value was in the water (positive function). The larger the infertile and green areas were in the river basins, the lower the total phosphorus concentration was in the water (negative function).

The effect of anthropogenic loads on the phosphate phosphorus concentration, (PO4-P is dependent variable Y), was calculated by multiple regression analysis. The results are presented in Table 6.



Dependent variable: TP, \* significance factor, *p* < 0.05.

**Table 6.** The influence of anthropogenic loads in river basins on the total phosphorus concentration in the water.


Dependent variable: PO4-P, \* significance factor, *p* < 0.05.

Multiple regression analysis of the influence of anthropogenic loads in basins on the concentration of phosphate phosphorus in the water showed that the PO4-P value was influenced by the discharge of municipal wastewater from background and transit pollution, agricultural land, forests, meadows, arable land, infertile land, and green land (*p* < 0.05). The higher the PO4-P concentration was in the municipal wastewater from background pollution, agricultural land, meadows, arable land, the higher the PO4-P value was in the water (positive function). The higher the transit pollution, and the larger the forests, infertile, and green areas were in the river basins, the lower the PO4-P concentration was in the water (negative function).

#### **4. Discussion**

Agricultural activity has strict negative impact on condition of surface water bodies, their ecosystems, the degradation of vegetation, and the quantitative and qualitative changes in fish populations in the Mediterranean basin [62]. The main factor affecting the Baltic Sea region environment is the increased amount of nutrients in rivers, mainly from diffuse agricultural sources [63]. The diffused nitrogen of anthropogenic origin account for about 70% of the total load deposited into rivers and lakes of the Baltic Sea basin area. Of the total diffuse load of nitrogen deposited into the Baltic Sea, 80% is from agriculture [64,65]. In Estonia, Latvia, and Lithuania, agriculture was intensified, and the amount of nitrogen fertilizers was increased after the 1990s [66].

Ikauniece and Lagzdinš assessed the status of two rivers, the Slocene and the Age, in Latvia. It was found that the ecological and chemical status of these rivers depended on the following factors: climatic conditions, types of soil and land-use, and human activities. The impact of land-use types and concentrations of total nitrogen, NO3 -N, NH4-N, total phosphorus, and PO4 −-P on the water of rivers was established. The highest concentrations of these substances were determined in the spring. It can be stated that snow melt during the spring period increases losses of biogenic compounds from concentrated sources [36].

An increase in sensitivity was found in basins with more agricultural land and more fertilizer. A change in the use of chemical fertilizers by ±20% affected the NO3-N loads in the water body between zero effect and an increase of ±13%, while a change in manure use by ±20% affected the NO3-N loads in the water body from zero effect to a change of −6% to +7% [63]. Ferrier [67] pointed out that nitrate concentrations in water of the rivers depend on the area of arable land, and there is a relationship between orthophosphate-P, suspended solids concentrations and meadow cover. Studies conducted in the Liswarta basin (Poland) have showed that high concentrations of nutrients in the Liswarta River and its tributaries are closely linked to the agriculture activities in this basin. However, urban wastewater effluents effected the highest concentrations of nutrients set in the Biała Oksza River [34]. Other Polish researchers assessed the influence of land use on the condition of the Dunajec, Czarny Dunajec, Biały Dunajec, and Białka rivers in the Podhale region (southern Poland). The results of their study showed that the concentrated pollution sources, such as effluents from WWTPs or untreated sewage from households, were more important than diffuse sources but agricultural activities significantly affect water quality of rivers [39].

Watershed modelling was used to discover the critical areas of water quality of rivers, as well as to define impacts and identify the most significant pollution sources in the river basins in Lithuania. Regional diffuse pollution leaching patterns were estimated using this model. The largest leaching rates total phosphorus were assessed in the southeastern and western parts of Lithuania. The largest leaching of total nitrogen was determined to occur in the center of the country. It can be seen from the modelling results that agriculture is the dominant pollution source in all Lithuanian river basins. The organic loads from diffuse pollution sources accounted for 60–90% of the annual loads in all of the river basins, excluding the urban catchments of the Neris and Nemunas rivers. The total phosphorus loads from agricultural sources accounted for 50–93% of the annual TP load. The pollution from concentrated sources and non-sewered households had almost no influence on the nitrate loads, and agriculture was the only dominant source of pollution, contributing

90–99% of the annual nitrate load [68]. It was determined that, satisfactorily, 90% of all nitrogen entered the Muša sub-catchment from the diffuse pollution sources, including 87% ¯ from the arable land and just a little more than 3% from the forest territory and pastures. A total of 10% of all nitrogen in the basin came from the concentrated pollution sources. The largest amounts of total phosphorus in the Muša sub-catchment entered the basin from ¯ the concentrated pollution sources (about 49%), arable land (36%) and about 15% from the forest area and pastures [69].

Various sources indicate measures for protection against diffuse pollution. In Poland, the recommendations for the protection of river valleys from biogenic pollution include the activities such as preserving natural vegetation on the banks of rivers, reducing of intensive agriculture activities and others [34]. Scholz [70] introduced diffuse pollution control strategies involving draining the natural wetlands by ditches in Germany.

In Lithuania, the main measures that should be applied to reduce the input of pollution from agricultural activities into rivers and other inland waters are as follows [71]:


#### **5. Conclusions**


**Author Contributions:** L.C. designed the study and performed the experiments; L. ˇ C., D.Š. and M.D. ˇ performed the experiments, analyzed the data, and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

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

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

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data are available in a publicly accessible repository.

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

#### **References**


## *Article* **Attribution Analysis of Climate and Anthropic Factors on Runoff and Vegetation Changes in the Source Area of the Yangtze River from 1982 to 2016**

**Guangxing Ji 1, Huiyun Song 1, Hejie Wei <sup>1</sup> and Leying Wu 2,\***


**Abstract:** Analyzing the temporal variation of runoff and vegetation and quantifying the impact of anthropic factors and climate change on vegetation and runoff variation in the source area of the Yangtze River (SAYR), is of great significance for the scientific response to the ecological protection of the region. Therefore, the Budyko hypothesis method and multiple linear regression method were used to quantitatively calculate the contribution rates of climate change and anthropic factors to runoff and vegetation change in the SAYR. It was found that: (1) The runoff, NDVI, precipitation, and potential evaporation in the SAYR from 1982 to 2016 all showed an increasing trend. (2) The mutation year of runoff data from 1982 to 2016 in the SAYR is 2004, and the mutation year of NDVI data from 1982 to 2016 in the SAYR is 1998. (3) The contribution rates of precipitation, potential evaporation and anthropic factors to runoff change of the SAYR are 75.98%, −9.35%, and 33.37%, respectively. (4) The contribution rates of climatic factors and anthropic factors to vegetation change of the SAYR are 38.56% and 61.44%, respectively.

**Keywords:** runoff variation; vegetation change; attribution analysis; source region of the Yangtze River

#### **1. Introduction**

The source area of the Yangtze River (SAYR) is located in the hinterland of the Qinghai Tibet Plateau, which is an important ecological barrier in China [1]. It has the ecological functions of maintaining the ecological security of the Yangtze River Basin [2]. At the same time, it is also one of the most sensitive and fragile regions in the global ecological environment [3,4]. Due to the joint impact of climate change and anthropic factors (indiscriminating felling of trees, disafforestation, overgrazing, mining, and digging), the ecological environment of the SAYR is deteriorating, and the function of water conservation is weakening [5,6]. Additionally, the problems of grassland degradation, wetland shrinkage, glacier retreat, soil erosion, and biodiversity decline are becoming increasingly prominent. These issues highlight the threats to the water resources and ecological security of the whole basin.

The impact of climate change has been embodied in the components of the hydrological cycle, including precipitation [7], evapotranspiration [8], and runoff [9]. Most of the literature has focused on precipitation and evapotranspiration worldwide [10,11], and runoff has taken less attention. Nevertheless, runoff remains one of the most important water resources and has a great influence on the formation of geomorphology, the development of soil, and the growth of plants [12]. It is an important condition for the regional industrial and agricultural water supply and plays an important role in the development of social economy.

**Citation:** Ji, G.; Song, H.; Wei, H.; Wu, L. Attribution Analysis of Climate and Anthropic Factors on Runoff and Vegetation Changes in the Source Area of the Yangtze River from 1982 to 2016. *Land* **2021**, *10*, 612. https://doi.org/10.3390/ land10060612

Academic Editors: Dawid Bedla, Marek Ryczek, Artur Radecki-Pawlik and Wiktor Halecki

Received: 6 May 2021 Accepted: 4 June 2021 Published: 8 June 2021

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

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

In recent years, under the joint influence of climate change and anthropic factors, the problems of climate, ecology, and hydrology in the SAYR have been a cause for concern. Therefore, numerous programs have been carried out to improve the eco-hydrological environment, among which Grain for Green proposed by the Chinese government in 1998 is the most well-known. In this context, monitoring the dynamic changes of runoff and vegetation and assessing the impacts of anthropic factors and climate change on runoff and vegetation change are of great practical value for formulating rational adaptive management countermeasures [13–19].

Previous studies have focused on the relationship between runoff and vegetation change and its influencing factors in the SAYR. Chen et al. analyzed the temporal and spatial change of vegetation cover in the SAYR from 1982 to 2003 and analyzed the influence of topography and human factors on the vegetation change [20]. Yao et al. analyzed the spatiotemporal variation characteristics of vegetation net primary productivity in the SAYR from 1959 to 2008 and analyzed the relationship between vegetation change and climate factors (temperature, relative humidity, precipitation, wind speed, and sunshine hours) [21]. Liu et al. analyzed the spatiotemporal variation characteristics of vegetation coverage in the SAYR from 1997 to 2012 and believed that the improvement of vegetation status in the SAYR benefited from the joint influences of climate humidification and ecological engineering [22]. Qian et al., Li et al., and Luo et al. studied the variation regularities of runoff in the source area from 1957 to 2009, 1961 to 2011, and 1961 to 2016, respectively, and pointed out that the annual surface water resources in the SAYR all showed an increasing trend [23–25]. Xu et al. analyzed the runoff variation characteristics and the degree of human impact in the SAYR from 1956 to 2004 [26]. Liu et al. separated the influences of climate and anthropic factors on runoff in the SAYR based on the Budyko hypothesis [27]. These research results have important scientific value for understanding the regularities of runoff and vegetation change in the SAYR, but there are few studies on the assessment and quantitative analysis of climate change and the contribution rate of anthropic factors to runoff and vegetation change in the SAYR [27,28].

The objective of this study was to quantify the impact of anthropic factors and climate change on vegetation and runoff variation in the SAYR by following three steps: (1) analyzing the temporal variation regularities of runoff and vegetation in the SAYR; (2) distinguishing the mutation year of runoff and NDVI data with the Mann–Kendall mutation analysis method; and (3) quantifying the contribution rate of climatic and anthropic factors to vegetation and runoff in the SAYR.

#### **2. Study Area and Data**

#### *2.1. Study Area*

The source region of the Yangtze River is located between 90◦43 and 96◦45 E, 32◦30 and 35◦35 N, with an area of about 13.77 × 104 km2 (Figure 1). Its annual average runoff is about 129.17 × <sup>10</sup><sup>8</sup> <sup>m</sup>3, and it affects the interannual fluctuation, long-term evolution trend, and the sustainable utilization of water resources in the Yangtze River Basin. The landform is mainly high plains and hills, with an average altitude of more than 4000 m. Its climate type belongs to the semi-humid and semi-arid region of the plateau frigid zone, with abundant sunshine and large temperature difference between day and night. Temperatures are generally low throughout the year in the Yangtze River source area. The month with the lowest precipitation is December, the month with the highest precipitation is July, and the annual precipitation is concentrated in May–September, which accounts for more than 85% of the annual precipitation [29].

**Figure 1.** The location of hydrological and meteorological stations in and around the source region of the Yangtze River Basin.

#### *2.2. Data Sources*

The data used in this study consist of three parts: (1) The annual runoff observation data of Zhimenda station from 1982 to 2016 were all from the Yangtze River Water Conservancy Commission (http://www.cjw.gov.cn/, accessed on 1 January 2021); (2) the climate station data in and around the source region of the Yangtze River Basin from 1982 to 2016 were obtained from the China Meteorological Administration (http://www.cma.gov.cn, accessed on 1 January 2021). First, the daily potential evaporation of 25 meteorological stations was calculated using the Penman–Monteith equation, and then the monthly precipitation and potential evaporation of 25 meteorological stations were calculated. Finally, the monthly precipitation and potential evaporation were interpolated from data of 25 meteorological stations in and around the SRYR by kriging. Annual precipitation and potential evaporation were obtained by adding monthly scale data. (3) The NDVI data of the source region of the Yangtze River Basin from 1982 to 2016 were obtained from the NOAA CDR AVHRR NDVI: Normalized Difference Vegetation Index, Version 5 (https://www.ncdc.noaa.gov/cdr/terrestrial/normalized-difference-vegetation-index, accessed on 1 January 2021). The time resolution is daily, and the spatial resolution is 0.05◦ × 0.05◦. The monthly NDVI and annual NDVI were extracted with the maximum combination method.

#### **3. Research Methods**

#### *3.1. Trend Analysis Method*

In this study, the univariate linear regression method (*Slope*) was used to analyze the variation trend of runoff, precipitation, potential evaporation, and NDVI in the study period; time t was taken as the independent variable. The calculation equation is as follows [30–32]:

$$Slope = \frac{n\sum\_{i=1}^{n} iX\_i - \sum\_{i=1}^{n} i\sum\_{i=1}^{n} X\_i}{n\sum\_{i=1}^{n} i^2 - \left(\sum\_{i=1}^{n} i\right)^2} \tag{1}$$

where *Slope* is the simple linear regression coefficient of runoff, precipitation, potential evaporation, and NDVI from 1982 to 2016. If *Slope* is positive, it means that the variable has an upward trend in the study period, whereas a negative value represents a downward trend; *n* is the number of years of the study period, *i* is 1 to *n*, *n* is 35, and *X* represents the annual runoff, precipitation, potential evaporation, and NDVI in the study area.

#### *3.2. Mann–Kendall Mutation Analysis Method*

Mann–Kendall (MK) mutation analysis has been widely used to test the abrupt point of hydro-meteorological elements due to its advantages of minimized interference from outliers and simple calculation [33,34].

For *X* with n data, a rank sequence is constructed, as shown in Equations (2) and (3):

$$Sk = \sum\_{i=1}^{k} ri \ (k = 2, 3, \dots, i) \tag{2}$$

$$
\dot{m} = \begin{cases} +1 & \text{xí} > \text{xj} \\ 0 & \text{xí} \le \text{xj} \end{cases} \\ \text{(j = 2, 3, ..., i)} \tag{3}
$$

The *UFk* statistics can be calculated by Equation (4):

$$UF\_k = \frac{[sk - E(sk)]}{\sqrt{Var(sk)}} (k = 1, \ 2, \ \dots, n) \tag{4}$$

where *UF*<sup>1</sup> = 0, *E(sk)* and *Var(sk)* represent the mean and variance of *sk*, respectively. The calculation equations are as follows:

$$E(\text{sk}) = \begin{array}{c} \frac{n(n+1)}{4} \quad (2 \le k \le n) \end{array} \tag{5}$$

$$Var(sk) = \begin{cases} \frac{n(n-1)(2n+5)}{72} & (2 \le k \le n) \end{cases} \tag{6}$$

The inverse order of *X* is calculated again, and at the same time, *UBk* = *UFk* (*k=n*, *n* − 1, ... ,1), *UB1* = 0, then the *UBk* curve can be calculated. When the two lines (*UF* curve and *UB* curve) intersect, and the intersection point is within the range of the 0.05 significance level, the intersection point can be considered as the mutation year.

#### *3.3. Budyko Hypothesis*

There are three assumptions in the attribution analysis of runoff change using the Budyko hypothesis: (1) Human activities and climate change do not affect each other and are independent factors. (2) For multi-year water balance, the change of water storage is usually negligible compared with the average annual precipitation depth. Therefore, the change of catchment storage water for multi-year water balance is assumed as 0; (3) The base period is only affected by climate change. Therefore, except for climate change, other factors affecting runoff change are classified as human activities.

The water balance equation of the basin is expressed as follows:

$$R = P - ET - \Delta S \tag{7}$$

*R* is the runoff depth; *P* is precipitation; *ET* is the actual evapotranspiration; and Δ*S* is the change of water storage. For multi-year water balance, the change of water storage is usually negligible compared with the average annual precipitation depth. Therefore, in the study of long-term hydrological data, it is generally considered that Δ*S* is 0.

The Budyko hypothesis, based on an assumption that the change of catchment storage water for multi-year water balance is considered as 0, is widely applied to quantitatively calculate the impact of climate change and anthropic factors on long-term annual runoff change [35–39]. *ET* can be calculated according to the Budyko hypothesis, the calculation equation of which is as follows [40,41]:

$$ET = \frac{P \times ET0}{\left(P^{\omega} + ET0^{\omega}\right)^{1/\omega}}\tag{8}$$

*ET*<sup>0</sup> is the annual average potential evapotranspiration (mm); *ω* indicates characteristic parameters of underlying surface.

*ET*<sup>0</sup> can be calculated using the Penman–Monteith equation.

$$ET\_0 = \frac{0.408\Delta (R\_{\rm ll} - G) + \gamma \frac{900}{T + 273} \mathcal{U}\_2 (\varepsilon\_a - \varepsilon\_d)}{\Delta + \gamma (1 + 0.34 \mathcal{U}\_2)} \tag{9}$$

Equation (7) can be transformed into Equation (10):

$$R = P - \frac{P \times ET0}{\left(P^{\omega} + ET0^{\omega}\right)^{1/\omega}}\tag{10}$$

where *R*, *P*, and *ET*<sup>0</sup> are known, and parameter *ω* can be calculated by the "fsolve" function of MATLAB.

Based on Equation (10), the elastic coefficients of precipitation (*εP*), potential evaporation (*εET*0), and underlying surface parameters (*εω*) to runoff can be calculated by the following equations [42]. The concept of elastic coefficient is shown in Appendix A.

$$\varepsilon\_P = \frac{(1+\phi^{\omega})^{1/\omega+1} - \phi^{\omega+1}}{(1+\phi^{\omega})\left[(1+\phi^{-\omega})^{1/\omega} - \phi\right]} \tag{11}$$

$$\varepsilon\_{ET0} = \frac{1}{(1 + \phi^{\omega}) \left[1 - (1 + \phi^{-\omega})^{1/\omega}\right]} \tag{12}$$

$$\varepsilon\_{\omega} = \frac{\ln(1 + \phi^{\omega}) + \phi^{\omega}\ln(1 + \phi^{\omega})}{\omega(1 + \phi^{\omega})\left[1 - (1 + \phi^{-\omega})^{1/\omega}\right]} \tag{13}$$

$$
\phi = \frac{ET0}{P} \tag{14}
$$

The runoff variation amount caused by the annual average precipitation change (Δ*RP*), annual average potential evaporation change (Δ*RET*0), and underlying surface parameter change (Δ*Rω*) can be calculated by the following Equations (15)–(17):

$$
\Delta \text{RP} = \varepsilon P \frac{R}{P} \times \Delta P \tag{15}
$$

$$
\Delta \text{RET}0 = \varepsilon ET0 \frac{R}{ET0} \times \Delta ET0 \tag{16}
$$

$$
\Delta R \omega = \varepsilon \omega \frac{R}{\omega} \times \Delta \omega \tag{17}
$$

On this basis, the contribution rate of precipitation (*ηRp*), potential evaporation, (*ηRET*0) and anthropic factor (*ηRH*) to runoff changes can be calculated using Equations (19)–(21):

$$
\Delta \mathcal{R} = \Delta \mathcal{R} P + \Delta \mathcal{R}ET0 + \Delta \mathcal{R} \omega \tag{18}
$$

$$
\eta R\_P = \Delta R P / \Delta R \times 100\% \tag{19}
$$

$$
\eta R\_{ET0} = \Delta RET0/\Delta R \times 100\% \tag{20}
$$

$$
\eta R\_H = \Delta R \omega / \Delta R \times 100\% \tag{21}
$$

#### *3.4. Attribution Analysis of Climate and Anthropic Factors to Vegetation Change*

There are two assumptions in the attribution analysis of vegetation change using this method: (1) Human activities and climate change do not affect each other and are independent factors. (2) The base period is only affected by climate change. Therefore, except for climate change, other factors affecting vegetation change are classified as human activities.

According to the mutation analysis results of NDVI time series data, the NDVI data are divided into two parts: base period (*T*1) and change period (*T*2). On this basis, the change of average NDVI in the two periods (*NDV I*) can be calculated as follows:

$$
\triangle NDVI = NDVI\_{T2} - NDVI\_{T1} \tag{22}
$$

where *NDV IT*<sup>1</sup> and *NDV IT*<sup>2</sup> are the average NDVI values in the base period (*T*1) and change period (*T*2), respectively. This method assumes that NDVI in the base period is only affected by climate change. Therefore, the NDVI difference between the base period and change period can be attributed to climatic factors (*NDV IC*) and anthropogenic factors (*NDV IH*).

The monthly NDVI values are quantitatively correlated with monthly precipitation and potential evapotranspiration [43,44]. Therefore, a multiple linear regression equation between monthly NDVI, monthly precipitation, and monthly potential evapotranspiration in the base period (*T*1) can be established:

$$
\overline{NDVI\_{T1}} = a \ast PT1 + b \ast ET0T1 + c \tag{23}
$$

The monthly precipitation (*PT*2) and potential evapotranspiration (*ET*0*T*2) in the *T2* period are known, and the simulated NDVI in the *T2* period (*NDV IT*2,*S*) can be calculated by Equation (24).

$$\overline{NDVI\_{T2,S}} = a\*PT2 + b\*ET0T2 + c \tag{24}$$

This method assumes that NDVI in the base period (*NDV IT*1) is only affected by climate change. Therefore, *NDV IT*2,*<sup>S</sup>* in the *T*<sup>2</sup> period calculated by Equation (24) is only affected by climate change. *NDV IT*<sup>2</sup> is the average value of NDVI observation in the *T*<sup>2</sup> period, which is affected by the combination of climate change and human activities.

The contribution rates of anthropic factor (*ηNDV IH*) and climate change (*ηNDV IC*) to NDVI can be calculated separately by Equations (25)–(28):

$$
\triangle\,\overline{NDVI\_H} = \overline{NDVI\_{T2}} - \overline{NDVI\_{T2,S}}\tag{25}
$$

$$
\triangle\,\overline{NDVI\_{\mathbb{C}}} = \overline{NDVI\_{T2,\mathbb{S}}} - \overline{NDVI\_{T1}}\tag{26}
$$

$$
\eta \text{NDVI}\_{\text{H}} = \triangle \overline{\text{NDVI}\_{\text{H}}} / \triangle \overline{\text{NDVI}} \tag{27}
$$

$$
\eta N DVI\_{\mathbb{C}} = \triangle \overline{NDVI\_{\mathbb{C}}} / \triangle \overline{NDVI} \tag{28}
$$

#### **4. Results and Analysis**

*4.1. Trends Analysis of Runoff, NDVI, and Climate Factors*

Figure 2 displays the change trend of the annual runoff in Zhimenda hydrological station and NDVI in the SAYR. It can be found that the annual runoff of Zhimenda hydrological station shows an increasing trend from 1982 to 2018, with an average annual increase of 0.8831 × 108 m3. The NDVI value in the SAYR showed a significant growth trend from 1982 to 2016, with an average annual growth of 0.0018.

**Figure 2.** Change trend of annual runoff (**a**) at Zhimenda hydrological station and NDVI (**b**) in the SAYR.

Figure 3 displays the change trend of average annual precipitation and annual potential evaporation in the SAYR. It can be found from Figure 3a that the annual average precipitation in the SAYR showed an increasing trend from 1982 to 2016, with an average annual increase of 1.2394 mm. Figure 3b displays that the average annual potential evaporation in the SAYR also showed an upward trend from 1982 to 2016, with an average annual increase of 0.5965 mm.

**Figure 3.** Variation trend of average precipitation (**a**) and average potential evaporation (**b**) in the SAYR.

#### *4.2. Mutation Analysis of Runoff and NDVI*

In this study, the Mann–Kendall method was used to analyze runoff and NDVI data in the SAYR, and the results are shown in Figures 4 and 5.

It can be found from Figure 4 that the UF and UB curve calculated from the runoff time series data in the SAYR intersected in 2004, and this intersection point was within the range of the 0.05 significance level, indicating that 2004 was the abrupt change year of runoff in the SAYR from 1982 to 2016.

It can be found from Figure 5 that UF and UB curve calculated from the NDVI time series data in the SAYR intersected in 1998, and this intersection point was within the range of the 0.05 significance level, indicating that 1998 was abrupt change year of NDVI in the SAYR from 1982 to 2016. This may be because Grain for Green was proposed by Chinese government in 1998, and the impact of human activities on vegetation change has increased since 1999.

**Figure 4.** Mann–Kendall mutation test result of runoff in the SAYR from 1982 to 2016.

**Figure 5.** Mann–Kendall mutation test result of NDVI in the SAYR from 1982 to 2016.

#### *4.3. Assessment of the Contribution Rates of Climate and Anthropic Factors to Runoff and Vegetation Changes*

Identifying the abrupt change year of runoff and NDVI was the premise for dividing time series data into different periods (i.e., base period and changing period). If there was no abrupt change year of runoff and NDVI data, we did not divide time series data into different periods.

According to mutation analysis results of runoff time series data at Zhimenda hydrological station from 1982 to 2016, the Zhimenda hydrological station 1982–2016 data were divided into two periods: the base period (1982–2004) and the change period (2005–2016), and the eigenvalues of meteorological and hydrological variables for different periods in the Zhimenda hydrological station were calculated (Table 1).

**Table 1.** Eigenvalues of meteorological and hydrological variables in the SAYR.


In order to further evaluate the temporal evolution characteristics of the impact of climate factors and underlying surface parameters on runoff, according to

Equations (11)–(13), the elastic coefficients of precipitation, potential evaporation, and underlying surface characteristic parameters in different years on runoff change of Zhimenda hydrological station were calculated (Figure 6).

**Figure 6.** Change trend of precipitation elasticity, potential evapotranspiration elasticity, and underlying surface parameter elasticity of runoff in the SAYR from 1982–2016.

It can be seen from Figure 6 that the absolute values of *ε*P, *ε*ET0, and *εω* in for Zhimenda hydrological station all showed a decreasing trend, indicating that the sensitivity of runoff change in the SAYR to climate factors and underlying surface parameters is decreasing.

According to Table 1, the differences of average annual precipitation, potential evaporation, and underlying surface parameter in the base period (1982–2004) and change period (2005–2016) were calculated and are displayed in Table 2. Then, according to Equations (15)–(17), the runoff variation amount for Zhimenda hydrological station caused by climate factors (precipitation, potential evaporation) and underlying surface condition change were calculated, as displayed in Table 2. Finally, the contribution rate of climate (precipitation, evaporation) and anthropic factors to runoff change of Zhimenda hydrological station was calculated according to Equations (18)–(21), as shown in Table 2.

**Table 2.** Attribution analysis of runoff variation in the SAYR.


The contribution rates of climate and anthropic factors to runoff variation in Zhimenda hydrological station are 66.63% and 33.37%, respectively. In general, climate change plays a major role in increasing runoff, and precipitation has a more significant effect on runoff increase than the reference evapotranspiration.

This conclusion is similar to that of Liu et al. [27]. However, there are some differences in the specific contribution rate value, and this may be due to the following factors: (1) The time ranges of runoff data are different; (2) different methods are used to calculate the contribution rate of climatic (precipitation, potential evaporation) and anthropic factors to runoff variation.

According to the abrupt analysis results of NDVI data of the SAYR from 1982 to 2016, the NDVI data in the SAYR from 1982 to 2016 can be divided into two periods: the base period (1982–1998) and the change period (1999–2016). Previous studies have found that

monthly NDVI is closely related to monthly precipitation and potential evaporation [43,44], and NDVI change can be used to characterize vegetation change. Referring to the method proposed by Li et al. [45], this study first established the multiple linear regression equation between monthly NDVI, monthly precipitation, and monthly potential evapotranspiration in the period of 1982–1998, and then the contribution rate of climate change and anthropic factors to vegetation change was quantitatively analyzed according to Equations (25)–(28) (Table 3). The results showed that compared to the period of 1982–1998, anthropic factors played a major role in increasing vegetation coverage in the period of 1999–2016, with a contribution rate of 61.44%, and the contribution rate of climate factors was 38.56%.


**Table 3.** Attribution analysis of vegetation change in the SAYR.

#### **5. Discussion**

Although the change trend and mutation characteristics of runoff and vegetation in the SAYR were analyzed, and the impact of climate change and anthropic factors on runoff and vegetation change were calculated, there are still some uncertainties. The data of precipitation and potential evaporation in the study area are interpolated from low-density meteorological station data in and around the SAYR, but there are still some deviations between the interpolation results and their actual distribution. The change of catchment storage water for multi-year water balance is assumed as 0, and we assumed that human activities and climate change do not affect each other and are independent factors. This study also assumed that the base period is only affected by climate change. All of these factors will lead to some uncertainties in the research results.

#### **6. Conclusions**

This study revealed the change trend and mutation characteristics of runoff and vegetation in the SAYR, and the contribution rates of climate change and anthropic factors to runoff change in the SAYR were quantitatively calculated using the Budyko hypothesis and multiple linear regression method, which provides theoretical support for water resource management and ecological protection in the SAYR.

The results showed that the runoff, NDVI, precipitation, and potential evaporation in the SAYR from 1982 to 2016 all showed an increasing trend. The abrupt change years of runoff and NDVI data in the SAYR from 1982 to 2016 were 2004 and 1998, respectively. Anthropic factors play a major role in annual runoff and vegetation change, with contribution rates of 75.98% and 61.44%. In the follow-up study, we will try to quantitatively analyze the contribution rate of climate and human factors to vegetation and runoff changes in different seasons.

**Author Contributions:** Conceptualization, G.J. and L.W.; Data curation, G.J., H.S. and H.W.; Funding acquisition, L.W. and H.W.; Methodology, G.J. and H.S.; Project administration, L.W. and H.W.; Writing—original draft, G.J., H.S. and L.W.; Writing—review and editing, G.J., H.W. and L.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 41901239), Second Tibetan Plateau Scientific Expedition and Research Program (Grant No. 2019QZKK0608), Soft Science Research of Technology Development Projects of Henan Province (Grant No. 192400410085) and China Postdoctoral Science Foundation (Grant No. 2018M640670).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **Appendix A**

Based on the abrupt change years of runoff in the SAYR, the runoff time series data are divided into two periods: the base period and the change period. The average annual precipitation in the base period is expressed as *P*1, and the average annual precipitation in the change period is expressed as *P*2. The variation of annual precipitation can be expressed as:

$$
\Delta P = P2 - P1\tag{A1}
$$

Similarly, the variation of the potential evaporation (Δ*ET*0) and characteristic parameters of underlying surface (Δ*ω*) can be calculated by Equations (A2) and (A3):

$$
\Delta ET0 = ET0\_2 - ET0\_1\tag{A2}
$$

$$
\Delta\omega = \omega 2 - \omega 1\tag{A3}
$$

The elasticity coefficient is the ratio of the change rate of two interrelated indexes in a certain period. The elastic coefficients of precipitation (*εP*), potential evaporation (*εET*0), and underlying surface parameters (*εω*) to runoff can be expressed by Equations (A4)–(A6), respectively.

$$
\varepsilon\_P = \frac{\Delta P}{P} / \frac{\Delta R}{R} \tag{A4}
$$

$$
\varepsilon\_{ET0} = \frac{\Delta ET0}{ET0} / \frac{\Delta R}{R} \tag{A5}
$$

$$
\varepsilon\_{\omega} = \frac{\Delta\omega}{\omega} / \frac{\Delta R}{R} \tag{A6}
$$

#### **References**


## *Article* **Quantitatively Calculating the Contribution of Vegetation Variation to Runoff in the Middle Reaches of Yellow River Using an Adjusted Budyko Formula**

**Guangxing Ji 1, Junchang Huang 1, Yulong Guo <sup>1</sup> and Dan Yan 2,3,\***

<sup>1</sup> College of Resources and Environmental Sciences, Henan Agricultural University, Zhengzhou 450002, China; guangxingji@henau.edu.cn (G.J.); huangjc1014@henau.edu.cn (J.H.); gyl.zh@henau.edu.cn (Y.G.)

<sup>2</sup> Center for Energy, Environment & Economy Research, Zhengzhou University, No. 100 Science Avenue, Gaoxin District, Zhengzhou 450001, China


**Abstract:** The middle reaches of the Yellow River (MRYR) are a key area for carrying out China's vegetation restoration project. However, the impact of vegetation variation on runoff in the MRYR is still unclear. For quantitatively evaluating the contribution rate of vegetation variation to runoff in the MRYR, this paper quantified the relationship between Normalized Difference Vegetation Index (NDVI) and Budyko parameters (*w*). Then, we used multiple linear regression to quantitatively calculate the contribution rate of different factors on vegetation variation. Finally, an adjusted Budyko formula was constructed to quantitatively calculate the influence of vegetation variation on runoff. The results showed that there is a linear relationship between NDVI and Budyko parameters (*w*) (*p* < 0.05); the fitting parameter and constant term were 12.327 and −0.992, respectively. Vegetation change accounted for 33.37% in the MRYR. The contribution of climatic and non-climatic factors on vegetation change is about 1:99. The contribution of precipitation, potential evaporation, anthropogenic activities on the runoff variation in the MRYR are 23.07%, 13.85% and 29.71%, respectively.

**Keywords:** vegetation variation; runoff variation; Budyko hypothesis; attribution analysis

#### **1. Introduction**

In the past few decades, many ecological protection projects carried out in China have led to a significant increase in vegetation [1,2]. The middle reaches of the Yellow River (MRYR) are the key area for the implementation of the vegetation restoration project in China [3,4]. Vegetation index is a simple, effective and empirical measure of surface vegetation, so the Normalized Difference Vegetation Index (NDVI) was used to characterize the vegetation change in the MRYR. Some studies show that NDVI in the MRYR has increased significantly, and its vegetation coverage has been restored [5–7].

Vegetation change can significantly change runoff by changing hydrological processes, such as vegetation transpiration and interception evaporation, and then affect the available amount of watershed runoff [8–10]. The debate on the relationship between "vegetation and water" dates back to at least the mid-19th century [11]. Although most studies show that the increase in vegetation will lead to a decrease in water yield (called negative impact here), only a few studies show that the increase in vegetation has little impact on water yield [12,13] (called no impact here) and even leads to an increase in water yield [14,15] (called positive impact here). Interestingly, the areas with positive and no impact are mostly in humid areas and large watersheds with complex terrain. In the past few decades, the spatial variability and periodicity of precipitation and runoff in the MRYR have changed [16–19]. Therefore, quantitatively calculating the influence of

**Citation:** Ji, G.; Huang, J.; Guo, Y.; Yan, D. Quantitatively Calculating the Contribution of Vegetation Variation to Runoff in the Middle Reaches of Yellow River Using an Adjusted Budyko Formula. *Land* **2022**, *11*, 535. https://doi.org/ 10.3390/land11040535

Academic Editors: Wiktor Halecki, Dawid Bedla, Marek Ryczek and Artur Radecki-Pawlik

Received: 13 March 2022 Accepted: 5 April 2022 Published: 7 April 2022

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

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

vegetation variation on runoff variation is helpful to deeply understand the response mechanism of the water cycle process to the underlying surface condition changes, and has theoretical and practical significance for improving watershed water resources management measures. It is useful for the utilization and development of water resources in the Yellow River, the protection of the human living environment and the sustainable development of the economy.

Many scholars have carried out studies for quantitatively calculating the influence of vegetation variation on runoff variation in the MRYR [20–22]. The runoff decreased with the improvement of vegetation coverage, seen by analyzing the statistical relationship between vegetation coverage and runoff in most basins in the MRYR [23]. Liang et al. [24] quantitatively analyzed the influence of ecological restoration measures for the streamflow variation of most watersheds on the Loess Plateau from 1961 to 2009, based on Budyko theory, and found that ecological restoration measures are the leading factor for the decline in runoff in most watersheds. Li et al. [25] found that the Normalized Difference Vegetation Index (NDVI) in most basins in the Loess Plateau increased significantly after 2000, while the runoff decreased significantly. A higher evapotranspiration rate from restored vegetation is the primary reason for the reduced runoff coefficient. Zhang et al. [26] quantitatively evaluated the influence of vegetation cover variation on streamflow variation in the MRYR over the past 37 years. It was shown that vegetation variation is the leading driving force of streamflow decline from 2000 to 2010. Wang et al. [27] discovered that the vegetation coverage in the MRYR increased by 29.72% from 1982 to 2018. The average runoff at each station showed a downward trend from 1961 to 2018, and ecological restoration played an important role in the decline in runoff in the MRYR. However, the contribution of vegetation variation on streamflow variation in the MRYR has not been calculated quantitatively.

Therefore, in order to quantitatively evaluate the contribution rate of vegetation variation on streamflow variation in the middle reaches of the Yellow River (MRYR), this paper analyzes its hydro-meteorological characteristics from 1982 to 2015, identifies the mutation year of streamflow, and quantifies the relationship between the NDVI and a Budyko parameter (*w*). Then, multiple linear regression is adopted to distinguish the contribution rate of different factors on vegetation variation. Finally, an adjusted Budyko formula is constructed to quantitatively calculate the influence of vegetation variation on runoff in the MRYR. This work is useful to understand the influence of vegetation variation on the evolution of water resources in the MRYR, and also supply important insights for water resources management.

#### **2. Research Region and Data**

#### *2.1. Research Region*

The middle reaches of the Yellow River (MRYR) are in the area of the Toudaoguai-Huayuankou hydrological station (Figure 1). The length of this reach is about 1234.6 km, its area accounts for about 45.7% and its runoff accounts for about 44.3% of the Basin. Guanzhong Plain and Hetao Plain are in the MRYR, which produce a large number of agricultural products every year. The wheat output ranks first, soybean output ranks second, and cotton output ranks fourth in China. The decreasing runoff in the MRYR has brought great harm to agricultural production and food security.

#### *2.2. Data*

(1) The annual runoff observation data of Toudaoguai and Huayuankou hydrometric stations from 1982 to 2015 were obtained from the hydrological statistical yearbook of China. In this study, the annual runoff data in the MRYR was equal to a value, which was the annual measured runoff of Huayuankou hydrometric station minus the annual measured runoff data of Toudaoguai hydrological station.

(2) 82 meteorological stations data (1982–2015) were downloaded from China Meteorological Administration. We adopted the Penman–Monteith formula for calculating the potential evaporation of 82 meteorological stations, and then the annual precipitation and potential evaporation were interpolated by Kriging. Finally, the average annual precipitation and potential evaporation in the MRYR can be obtained.

$$ET\_0 = \frac{0.408\Delta (R\_n - G) + \gamma \frac{900}{T + 273} l l\_2 (\varepsilon\_a - \varepsilon\_d)}{\Delta + \gamma (1 + 0.34 l l\_2)}\tag{1}$$

(3) The NDVI dataset (8 × 8 km2) was sourced from the National Aeronautics and Space Administration. The time span of the dataset is from July 1981 to December 2015, its time resolution is 15 d, and its downloaded files are in NC4 format. The annual NDVI are calculated by the average value composites method, and is equal to the arithmetic mean of all grids in the study area. This method further eliminates some interference of cloud, atmosphere, solar altitude, and sensor sensitivity. Finally, the annual NDVI from 1982 to 2015 can be obtained.

**Figure 1.** Location of the study area and distribution of Hydro-meteorological stations.

#### **3. Methods**

#### *3.1. Trend Analysis Method (Slope)*

The formula of *Slope* method is as follows [28,29]:

$$Slope = \frac{n\sum\_{i=1}^{n} iXi - \sum\_{i=1}^{n} i\sum\_{i=1}^{n} Xi}{n\sum\_{i=1}^{n} i^2 - \left(\sum\_{i=1}^{n} i\right)^2} \tag{2}$$

If *Slope* is greater than 0, the variable shows ever-growing properties. If *Slope* < 0, it implies that the variable gets less with the day; *n* means the number of years; *X* indicates the annual precipitation, and potential evaporation runoff.

#### *3.2. Mann-Kendall (MK) Mutation Detection Algorithm*

For *X* containing *n* data, we first construct a rank sequence. Then, we calculated *Sk*.

$$S\_k = \sum\_{i=1}^k r\_i \quad (k = 2, 3, \dots, i) \tag{3}$$

$$r\_i = \begin{cases} +1 & \text{x}\_i > \text{x}\_j \\ 0 & \text{x}\_i \le \text{x}\_j \end{cases} \quad (j = 2, 3, \dots, i) \tag{4}$$

*E*(*Sk*) represents the average value of *Sk*, *Var*(*Sk*) represents the variance of *Sk*. Finally, *UF*<sup>1</sup> = 0, the *Ufk* statistics can be calculated.

$$E(S\_k) = \begin{array}{c} \frac{n(n+1)}{4} \end{array} \quad \text{( $2 \le k \le n$ )}\tag{5}$$

$$Var(S\_k) = \begin{cases} \frac{n(n-1)(2n+5)}{72} & (2 \le k \le n) \end{cases} \tag{6}$$

$$LIF\_k = \frac{|S\_k - E(S\_k)|}{\sqrt{Var(S\_k)}} \quad (k = 1, 2, \dots, i) \tag{7}$$

Constructing the reverse order of *X*, and making *UBk* = *UFk* (*k=n*, *n* − 1, ..., 1), *UB*<sup>1</sup> *=* 0, we calculated the *UBk* statistics in this way.

#### *3.3. Budyko Hypothesis*

There are two hypotheses in the attribution analysis of runoff change using Budyko hypothesis: (1) assuming that human activities, climate factors and vegetation do not affect each other, and are independent factors; (2) assuming that the base period is only affected by climate factors.

The multi-year scale water balance leads to Formula (8).

$$Q = P - E \tag{8}$$

where, *Q* represents runoff depth; *P* represents precipitation; *E* represents actual evaporation.

According to the Buduko hypothesis, there is a nonlinear function relation between *E*/*P* and dryness ratio (phi = *E*0/*P*) of a specific watershed.

$$\frac{E}{P} = F\left(\frac{E\_0}{P}\right) = F(\varphi) \tag{9}$$

*E*<sup>0</sup> is the potential evaporation (mm). The Fu-type Budyko equation is [30]:

$$F(\varphi) = 1 + \varphi + (1 + \varphi^w)^{1/w} \tag{10}$$

*w* represents the variable representing the underlying surface information of the basin, and is used to characterize human activities. Human activities can affect the runoff change of the Yellow River Basin through many aspects, including vegetation change, water conservancy project construction, domestic water for urban residents and agricultural irrigation water, etc. The larger the value of *w*, the less precipitation is transformed into runoff.

$$Q = P\left( (1 + \varphi^{\rm u})^{1/w} - \varphi \right) \tag{11}$$

Li et al. [31] found that there is a good linear relationship between *NDVI* and Budyko parameter (*w*) in the 26 major global river basins, which makes it possible to quantitatively evaluate the contribution rate of vegetation variation on runoff change.

$$w = a \times NDVI + b \tag{12}$$

*a* and *b* are constants for a specific watershed and can be calculated by the least squares method.

$$Q = P\left(\left(1 + \varphi^{a \times NDVI + b}\right)^{1/(a \times NDVI + b)} - \varphi\right) \tag{13}$$

The elasticity coefficient is used to measure the relative change of another variable caused by the change of one variable. *εP*, *εE*0, *ε<sup>w</sup>* and *εNDV I* represent the elastic coefficient of *P*, *E*0, *w* and NDVI on runoff, respectively.

$$\varepsilon\_P = 1 + \frac{qF'(\varphi)}{1 - F(\varphi)}\tag{14}$$

$$\varepsilon\_{Eo} = -\frac{\varrho F'(\varrho)}{1 - F(\varrho)}\tag{15}$$

$$\varepsilon\_{\mathcal{W}} = -\frac{wF'(w)}{1 - F(q)}\tag{16}$$

$$
\varepsilon\_{NDVI} = \varepsilon\_w \frac{a \times NDVI}{a \times NDVI + b} \tag{17}
$$

$$F'(\varphi) = 1 - \varphi^{w-1} (1 + \varphi^w)^{1/w - 1} \tag{18}$$

$$F'(w) = \frac{(1+q^w)^{1/w} \ln(1+q^w)}{w^2} - \frac{q^w (1+q^w)^{1/w-1} \ln q}{w} \tag{19}$$

Δ*RP*, Δ*RE*0, Δ*Rw*, Δ*RNDV I* and Δ*RH* represent the value of runoff change caused by precipitation, potential evapotranspiration, underlying surface change, vegetation change and human activities.

$$
\Delta R\_P = \varepsilon\_P \frac{R}{P} \times \Delta P \tag{20}
$$

$$
\Delta R\_{E0} = \varepsilon\_{E0} \frac{R}{E\_0} \times \Delta E\_0 \tag{21}
$$

$$
\Delta R\_{\text{\textw}} = \varepsilon\_{\text{\textw}} \frac{R}{w} \times \Delta w \tag{22}
$$

$$
\Delta R\_{NDVI} = \varepsilon\_{NDVI} \frac{R}{NDVI} \times \Delta N DVI \tag{23}
$$

$$
\Delta R\_H = \Delta R\_{\text{w}} - \Delta R\_{NDVI} \tag{24}
$$

Δ*P*, Δ*E*0, Δ*w* and Δ*NDV I* represent the changes of multi-year average *P*, *E*0, *w* and *NDVI* from *T*<sup>1</sup> to *T*<sup>2</sup> period.

Δ*NDV I* is attributed to *NDVI* change value caused by climate factors (Δ*NDV Ic*) and human activities (Δ*NDV IH*):

$$
\Delta NDVI = NDVI\_{T2} - NDVI\_{T1} = \Delta NDVI\_{\mathbb{C}} + \Delta NDVI\_{H} \tag{25}
$$

Average *NDVI* is significantly correlated with *P* and *E*<sup>0</sup> [32–34], so we calculated the multiple linear regression equation between *NDVI*, *P* and *E*<sup>0</sup> in *T*<sup>1</sup> period. There are two assumptions in the attribution analysis of vegetation change using multiple linear regression method: (1) Human activities and climate factors are independent factors. (2) The base period is only affected by climate factors. Therefore, except for climate factors, other factors affecting vegetation change are classified as human activities.

$$NDVI\_{T1} = a \times P\_{T1} + b \times E\_{0T1} + c \tag{26}$$

$$NDVI\_{T2,s} = a \times P\_{T2} + b \times E\_{0T2} + c \tag{27}$$

$$
\Delta NDVI\_H = NDVI\_{T2} - NDVI\_{T2,S} \tag{28}
$$

$$
\Delta NDVI\_{\rm c} = NDVI\_{T2,S} - NDVI\_{T1} \tag{29}
$$

$$
\eta\_{NDVI\_H} = \Delta NDVI\_H / \Delta NDVI \tag{30}
$$

$$
\eta\_{NDVI\_{\mathbb{C}}} = \Delta NDVI\_{\mathbb{C}} / \Delta NDVI \tag{31}
$$

where, *NDV IT*2,*<sup>S</sup>* represents average *NDVI* value under climate change in the *T*<sup>2</sup> period. *NDV IT*<sup>1</sup> and *NDV IT*<sup>2</sup> represent average *NDVI* value of the two periods, respectively. *ηNDV IH* and *ηNDV IC* represent the influence of anthropic and climatic factors on vegetation variation.

$$
\Delta R = \Delta R\_P + \Delta R\_{E0} + \eta\_{NDVI\_C} \times \Delta R\_{NDVI} + \eta\_{NDVI\_H} \times \Delta R\_{NDVI} + \Delta R\_H \tag{32}
$$

$$
\eta\_{Rp} = \Delta R\_P / \Delta R \times 100\% \tag{33}
$$

*ηRE*<sup>0</sup> = Δ*RE*0/Δ*R* × 100% (34)

$$
\eta\_{R\_{NDVI\mathbb{C}}} = \eta\_{NDVI\_{\mathbb{C}}} \times \Delta R\_{NDVI} / \Delta R \times 100\% \tag{35}
$$

$$
\eta\_{\text{R}\_{\text{NDVIH}}} = \eta\_{\text{NDVI}\_{\text{H}}} \times \Delta R\_{\text{NDVI}} / \Delta R \times 100\% \tag{36}
$$

$$
\eta\_{R\_H} = \Delta R\_H / \Delta R \times 100\% \tag{37}
$$

*ηRp* , *ηRE*<sup>0</sup> , *ηRNDV IC* , *ηRNDV IH* and *ηRH* represent the contribution rates of precipitation, potential evaporation, climatic vegetation change, anthropogenic vegetation change and anthropic factors on runoff.

#### **4. Results and Analysis**

#### *4.1. Trend Analysis of Meteorological and Hydrological Elements*

The changes in runoff, precipitation and potential evaporation in the MRYR from 1982 to 2015 are displayed in Figure 2. From Figure 2a, the annual runoff in the MRYR showed a significant decreasing tendency in the period of 1982–2015 (*p* < 0.05), with a slope of −3.3143 × <sup>10</sup><sup>8</sup> m3/a. From Figure 2b, the average annual precipitation in the MRYR showed a tendency to increase from 1982 to 2015 (*p* > 0.05), with a slope of 0.0582 mm/a. The average annual potential evaporation in the middle reaches of the Yellow River displayed a fluctuating growth tendency in the period of 1982–2015 (*p* < 0.05), and the increase rate was 1.6357 mm/a.

**Figure 2.** Changes in annual runoff (**a**), precipitation (**b**) and potential evaporation (**c**) in MRYR.

#### *4.2. Mutation Year Identification*

Because the Mann-Kendall mutation test algorithm is not disturbed by a few outliers and has the advantages of simple calculation, it is widely used to identify the mutation points of hydro-meteorological elements [35–38].

In the Mann-Kendall mutation detection algorithm, UF statistics are a sequence of statistics calculated by the order of time series *X*, and UB statistics are a sequence of statistics calculated by the inverse order of time series *X*. In practical application, it is generally considered that the result is credible at the significance level of 0.05. If the UF statistics line and UB statistics line have intersections and these intersections are within the range of 0.05, this is significant. The year corresponding to the intersection is the mutation year of the time series data.

UF statistics of runoff time series data in the MRYR intersected with UB statistics in 1989 (Figure 3), and this intersection is within the range of 0.05 (significant level), indicating that 1989 is a runoff mutation year in the MRYR.

**Figure 3.** Mann-Kendall catastrophe analysis result of runoff in MRYR from 1982 to 2015.

UF statistics of precipitation time series data in the MRYR intersected with UB statistics in 1984, 2010, 2012 (Figure 4), and these intersections are all within the range of the 0.05 significant level, indicating that 1984, 2010, 2012 are precipitation mutation years in the MRYR.

**Figure 4.** Mann-Kendall catastrophe analysis result of precipitation in MRYR from 1982 to 2015.

#### *4.3. The Relationship between NDVI and w*

In Budyko's hypothesis, the parameter (*w*) reflects the change in the underlying surface, so the parameter (*w*) should change on an annual or multi-year scale. For the Fu-type Budyko equation, if the elements (*Q*, *P*, *E*0) in the equation can be determined, the parameter (*w*) can be solved. Using the hydrological and meteorological data in the MRYR from 1982 to 2015, the parameters (*w*) from 1982 to 2015 are obtained. Figure 5

displayed the variation tendency of NDVI and underlying surface parameters (*w*) in the MRYR from 1982 to 2015. The NDVI and underlying surface parameters (*w*) in the middle reaches of the Yellow River increased significantly (*p* < 0.05), and the slopes were 0.017/10a and 0.268/10a, respectively.

**Figure 5.** Variation trend of NDVI (**a**) and underlying surface parameters (*w*) (**b**) in the middle reaches of the Yellow River from 1982 to 2015.

In this paper, the values of underlying surface parameters (*w*) and NDVI are processed by a 10-year moving average, and then the regression coefficients (a and b) are obtained by least square fitting, a = 12.327, b = −0.992. The R<sup>2</sup> of equation fitting is 0.4465, and is significant (*p* < 0.05), indicating that the fitting result is good (Figure 6).

**Figure 6.** The relationship between NDVI and underlying surface parameter (*w*).

#### *4.4. Quantifying the Influence of Vegetation Variation on Runoff*

Because we want to quantitatively calculate the contribution rate of different factors on runoff change (Figure 7), 1982–2015 is separated into *T*<sup>1</sup> (1982–1989) and *T*<sup>2</sup> (1990–2015), according to the catastrophe analysis results of runoff data in the MRYR. Firstly, we analyzed the correlation between *NDVI*, *P* and *E*<sup>0</sup> (Figure 8), implying that NDVI is significantly correlated with *P* and *E*0. Next, we calculated the linear function relation between *NDVI*, *P* and *E*<sup>0</sup> in *T*<sup>1</sup> period. Finally, we analyzed the contribution rate of different factors on vegetation change, using formulas 27–31 (Table 1). The results showed that the contribution of climatic and non-climatic factors on vegetation variation is close to 1:99.

**Figure 7.** A flow chart for calculating the contribution of different factors on vegetation change.

η

**Figure 8.** The relationship between NDVI and precipitation (**a**) and potential evaporation (**b**).


Shown in Table 2, compared with the *T*<sup>1</sup> period, the potential evaporation in the *T*<sup>2</sup> period increased by 41.70 mm, precipitation and runoff depth decreased by 25.89 mm and 24.49 mm, respectively, underlying surface parameters (*w*) and NDVI increased by 0.52 and 0.023, respectively.

**Table 2.** The value of climate, hydrology and NDVI variables in different periods in MRYR.


Table 3 displays the influence of climatic and anthropogenic factors on the annual runoff of the post-change period. Contrasted to *T*1, runoff variation amount caused by precipitation, potential evapotranspiration, vegetation variation and anthropogenic activities in the *T*<sup>2</sup> is −5.42 mm, −3.25 mm, −7.84 mm, and −6.98 mm. Their contributions were 23.07%, 13.85%, 33.37%, and 29.71%, respectively.

**Table 3.** Attribution analysis of runoff variation in MRYR.


#### **5. Conclusions and Discussion**

Using NDVI, meteorological data and measured runoff data, based on an adjusted Fu-type Budyko equation, this paper quantitatively analyzed the contribution of vegetation variation on runoff in the MRYR from 1982 to 2015. The results showed that the influence of climatic and anthropogenic factors on vegetation variation is about 1:99. Vegetation change played a major role on runoff in the MRYR. Precipitation, potential evaporation and human activities account for 23.07%, 13.85% and 29.71% in runoff variation, respectively.

Although this paper carries out strict quality control on the data and model, some non-determinacy remains. Firstly, precipitation and potential evaporation data are interpolated from the data of meteorological stations in and around the MRYR, which will bring some non-determinacy. Secondly, the precondition of the study is that all variables are independent and do not affect each other. For example, when calculating the partial derivative of multi-year evapotranspiration to runoff, it is required that the factors, such as vegetation change, precipitation and human activities, are independent. However, in practice, the underlying surface and climate system are a whole. These factors interact and affect each other, which will have some non-determinacy on the research results.

This paper quantitatively calculated the influence of vegetation variation on runoff in the MRYR, but in order to reveal the influence mechanism of vegetation on water cycle, it is necessary to analyze how vegetation affects the hydrological process of "precipitationinterception-infiltration-evapotranspiration-runoff". In the follow-up study, we will combine the distributed hydrological model HLMS with the PML model, coupling leaf area index information to build a distributed Hydrothermal Coupling Model [39–41], for clarifying the response of water heat balance to vegetation variation. In addition, in the process of quantitatively calculating the contribution rate of different factors on runoff change, in addition to climate factors and vegetation change, water conservancy project construction, urban residents' domestic water and agricultural irrigation water are classified as human activity factors. We will quantitatively calculate the contribution rate of urban water, industrial water, agricultural water and water conservancy project construction on runoff change in the follow-up study.

**Author Contributions:** Conceptualization: G.J. and D.Y.; Data curation: G.J. and J.H.; Methodology: G.J., Y.G. and D.Y.; Writing—original draft: G.J. and D.Y.; Writing—review and editing: G.J., J.H., Y.G. and D.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China (2021YFD1700900), Soft Science Research of Technology Development Projects of Henan Province (212400410077), Think tanks key research projects of Colleges and Universities of Henan Province (2022ZKYJ07), Young backbone teachers program of Henan Province (2019GGJS048) and the special fund for top talents in Henan Agricultural University (30501031).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** We declare no conflict of interest.

#### **References**


## *Concept Paper* **Water Diversion in the Valley of Mexico Basin: An Environmental Transformation That Caused the Desiccation of Lake Texcoco**

**Carolina Montero-Rosado, Enrique Ojeda-Trejo, Vicente Espinosa-Hernández \*, Demetrio Fernández-Reynoso , Miguel Caballero Deloya and Gerardo Sergio Benedicto Valdés**

> Colegio de Postgraduados, Campus Montecillo, Estado de México, Texcoco 56230, Mexico; montero.carolina@colpos.mx (C.M.-R.); enriqueot@colpos.mx (E.O.-T.); demetrio@colpos.mx (D.F.-R.); mcaballero@colpos.mx (M.C.D.); bsergio@colpos.mx (G.S.B.V.) **\*** Correspondence: vespinos@colpos.mx

**Abstract:** Mexico's basin is one of the most altered in the country, owing to the presence of the megalopolis of Mexico City. Lake Texcoco, which had the basin's biggest extension, dried up almost completely. The basin's evolution over time led in the formation of a megabasin in which water is transported from one source to another to serve the urban region and subsequently drained to prevent flooding. The major hydrotechnical works in Mexico Basin have been interpreted as a solution to the problem of flooding in Mexico City, but they were actually part of a much larger strategy of territorial appropriation by the Spanish colonists. The ecological imbalance that has resulted has sparked a variety of social issues. For the purpose of analyzing the environmental transformation of Lake Texcoco over the last 500 years, actors and processes that influenced specific moments in the country's history were identified; these elements showed the inexorable relationship between the lake and Mexico City. Subsequently, they were grouped by periods with similar trends in terms of the way in which society relates to and appropriates the natural environment of the lake. It was found that the critical moment for the desiccation of Lake Texcoco occurred during the Spanish colonial historical period as part of the redesign of the city; from then on, the same environmental imaginary prevailed century after century, shaped by social and economic factors. This study contributes to the literature on how urbanization affects natural resources by making an original theoretical contribution through an analysis based on political ecology, and it adds to the literature on how people use the prevailing federal area of the lake.

**Keywords:** Mexico City; historical political ecology; urban ecology; water spatial policy; environmental change

#### **1. Introduction**

Lacking natural drainage, in the XV century, the Valley of Mexico Basin was an elevated basin of five interconnected lakes whose levels rose during the rainy season and fell during the dry months [1]. On an island in Lake Texcoco, the capital of the Aztec empire was founded. Through the construction of a complex engineering system, they handled flooding as well as drought while also assuring food production, urban water supply, and navigation [2]. Lake Texcoco originally spanned nearly the whole basin; it shrunk from a surface area of 7868 square kilometers to only 15.83 square kilometers [1].

After the Spanish conquest, the Spaniards took the indigenous Tenochtitlan as their own capital and strove to organize a regional economy modeled on Iberian patterns, which required more extensive modification of the basin's hydrology. Finally, in 1608, they decided to drain the Mexico Basin to the Tula River [3]. After the independence of Mexico from Spain, the draining continued until 1910, when the Grand Drainage Canal was inaugurated and the desiccation of Lake Texcoco was almost complete [4].

**Citation:** Montero-Rosado, C.; Ojeda-Trejo, E.; Espinosa-Hernández, V.; Fernández-Reynoso, D.; Caballero Deloya, M.; Benedicto Valdés, G.S. Water Diversion in the Valley of Mexico Basin: An Environmental Transformation That Caused the Desiccation of Lake Texcoco. *Land* **2022**, *11*, 542. https://doi.org/ 10.3390/land11040542

Academic Editors: Wiktor Halecki, Dawid Bedla, Marek Ryczek and Artur Radecki-Pawlik

Received: 15 March 2022 Accepted: 5 April 2022 Published: 8 April 2022

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

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

To understand the case of the desiccation of Lake Texcoco, it is necessary to analyze the sociopolitical actors and processes that generated or exacerbated the environmental problem in its long-term historical dimension. The purpose is to inform decision makers and improve environmental conditions while, at the same time, promoting nature conservation [5,6]. Bailey and Bryant [7] also mention that environmental problems cannot be fully understood in isolation from the political and economic contexts within which they originate.

The draining of Lake Texcoco is immersed in power relations and decisions about water resource management, which are reflected in historical conflicts over the use and management of lake resources in the Valley of Mexico. From the early Spanish colonial bureaucrats who sought to lower the surface area of the lake, to the nineteenth and twentieth century hydrological experts, the power groups that can be identified as elites in this region have long been at odds with nature. The objectives were frequently contentious which included the eradication of indigenous culture and customs along the lakeshore; the expansion of arable land for agriculture; the promotion of capitalist economic development; the improvement of sanitation; and the parceling out of lots to poor and working-class residents.

Political ecology and historical analysis of decisions that led to environmental changes are essential tools for understanding processes that have occurred over long periods of time [8–13]. The historical approach is especially important for understanding problems involving highly modified environments such as Lake Texcoco, where the natural environment is not just the backdrop for human events, it evolves on its own, both independently and in response to human activity [14].

Finally, the regional context of the society–environment relationship allows us to better understand the scale of a socio-territorial problem, the networks of actors involved and their discourses around it; it also allows us to identify the economic and political dynamics that have an impact on local resource degradation [12]. Figure 1 shows the four stages identified in this research by the policies applied. It can be seen that the critical moment in which the degradation begins occurs during the Spanish colony.

**Figure 1.** Identified stages for environmental policies in Lake Texcoco history. Author elaboration based on critical moments.

Lake Texcoco's history provides a window into the dynamism of the political economy, interactions between society and environment, and the numerous uses and competing developments. This research focuses on identifying the players and processes that contributed to the lake's desiccation throughout its history, with the historical processes grouped into four distinct stages defined by the political backdrop and management decisions made about the lake's environment.

#### **2. Study Area: Former Lake Texcoco**

The study area comprises the area of the former Lake Texcoco basin, which is part of the Mexico Basin located at the highest point of the Mexican plateau (2198 to 5200 m above sea level) with an average area of 9220 square kilometers. It was originally an endorheic basin surrounded by mountain ranges and volcanoes with interconnected lakes that functioned as a system of communicating vessels [1].

At its peak, Lake Texcoco had a surface area of 7868 square kilometers, occupying almost the entire basin. Over time, the lake divided into other lakes, creating a system of fresh and brackish water. In the pre-Hispanic period, the lake complex consisted of four large lakes: Zumpango, Xaltocan and Texcoco to the north with salt water, and Xochimilco to the south with fresh water. Lake San Cristobal, Lake Mexico and Lake Chalco were created by the construction of dams in the lakes of Xaltocan, Texcoco, and Xochimilco, respectively (ibid).

The origin of the soils in the ex-lake basin is related to the mechanical deposition of sediments from alluvial and lacustrine processes; their evolution has been conditioned to the parental material (sediment), the accumulation of organic matter, fine materials and carbonate, which resulted in the formation of entisols to mollisols or vertisols. Artificial drainage and desalination accelerated the pedogenic processes and allowed their incorporation into agriculture [15].

The study area's primary issues are significant environmental degradation caused by deforestation, lake surface reduction, soil erosion, loss of terrestrial and aquatic habitats, river contamination, and uncontrolled urban sprawl. Due to overexploitation and depletion of aquifers, which mostly affects the eastern portion of the Valley of Mexico, differential sinking has resulted in the extraction and pumping of water from the Lerma and Cutzamala rivers to meet Valley of Mexico demand, affecting external basins [16].

Historically, human settlements in the basin had an intimate relationship with water, and the lake's desiccation was closely related to the development of urbanization and Mexico City's expansion.

#### *Overview of Existing Studies*

Water historiography in the Valley of Mexico has been addressed by various authors, prioritizing immediate environmental problems, engineering or political aspects [4,17–22]. Candiani's work on the environmental evolution of Mexico City throughout colonial times [11] is an extremely intriguing and illustrative example. Her research focuses on one of the most difficult and pioneering engineering projects in the Mexico Basin, the Huehuetoca tunnel, which was designed to desiccate the lakes surrounding Mexico City. It investigates the social motivations underlying the drainage project's numerous structures and technological choices but is specific to colonial Mexico.

Iracheta and Dávalos [22] take a descriptive and historical approach to water in Mexico's basin, compiling a wide bibliography on themes such as the lacustrine system, flooding and drainage, hydraulic techniques, environmental sanitation, productive applications, and urban supply. Another historical descriptive approach is that of Lattes [19], who compiles the stories of personalities such as Spanish friars from the colonial era and discuss the floods in Mexico City and the hydrotechnical efforts conducted until the stage of Porfirio Díaz. He also includes historical maps that are critical for the study of drainage in the Mexico basin.

Boyer and Wakild [20] examine President Lazaro Cárdenas's six-year presidency (1934–1940) through the lens of environmental history. The Cárdenas administration is most recognized for implementing the revolution's "promises" through land reform, popular organization, and nationalization. Regardless, such statements obscure nature. Social landscaping, as they call it, was a key component of Cárdenas' ambitious social and political agenda. This is an intriguing approach, but it is again centered on a limited moment in Mexico's history.

Tortolero [21] wrote a historical analysis of water and environmental history in the Mexican basin in which he argues that the history of water management and use can serve as a lens through which to analyze diverse spheres of human activity. It is analyzed from four perspectives: the history of production and consumption, cultural history, economic history, and environmental history. This investigation is critical for gaining a better understanding of what has occurred in the Mexico basin.

In his work entitled *Political Essay on New Spain* [18], Alexander Von Humboldt describes the situation of natural resources in Mexico at the time of the colony. He points out that cultivation of morality can only be done in the groove of freedom; freedom, as a rule of balance in society, also allows the interactions of economic liberalism, noble actions of the new humanitarianism, anti-slavery attitude, and the free desire for the idea of progress. Humboldt frequently demonstrates the absence of freedom in New Spain and the resulting ethical–economic evils: authoritarianism, backwardness, morals, and a lack of culture. This is a seminal work for studies on the transformation of the landscape.

During the Porfiriato period, it was feasible to compile a report detailing all drainage work completed to that point [4]. It is a two-volume study that includes maps and extensive information about the state of the lakes prior to, during, and after the works are completed. While it compiles information about the social and political contexts in which the works were conducted, its primary contribution is technical and administrative data.

From the 1880s to the 1940s, Vitz [23] traces the failure of different conceptions for the engineering and restructuring of the Basin of Mexico, as well as the workers, peasants, and popular sectors that negotiated, fought, or attempted to benefit from these changes. An elite group of engineers, urban planners, and technocrats devised a strategy to push capitalist urban development on a rebellious populace and environment. According to this author, the hydraulic operations necessary to complete the drainage of the basin's lakes, new forest laws that deprived people of their community woods, changes in land use following the Mexican revolution, and the emergence of working-class colonies following 1920 are all part of the same story: Mexico City's capitalist rise.

However, Mexico City's floods and the drainage of the basin encompass a broader range of government policies, lake system descriptions, hydraulic techniques, productive uses, environmental sanitation, legislation, and litigation that can be better understood over a longer period of time than the ones mentioned above.

During the pre-Hispanic period, the dominant social discourse was about the scarcity of drinking water and the need for water conduction works for irrigation and food production; and at the same time, the need to prevent floods in Tenochtitlan [2,24–26].

For the Spaniards and the later independent government of Mexico (16th to 20th centuries), the lake was a problem due to the constant flooding in the city, the insalubrity and the *tolvaneras* (dust storms) caused by the gradual drying up of the lake. As a result, the decision was made to construct works that would allow free circulation of water and drainage of the city out of the Mexican basin towards the Tula River, which was completed in 1910 with the inauguration of the Great Drainage Canal, a work that accelerated the drying up of Lake Texcoco in the following decades [4,17].

#### **3. Research Methodology and Theoretical Approach**

In this study, a historical political ecology approach is used, in which an adaptation of Blaike and Brookfield's decision-making in land management degradation [9] is proposed as a tool to identify the political–environmental relationship, exposing different environmental policies regarding Lake Texcoco management and grouping them by periods (see Scheme 1).

**Scheme 1.** Decision-making in historical Lake Texcoco for alterations to the hydrological regime.

Political ecology originated in the 1980s as an interdisciplinary topic based on the core concept that ecological change cannot be understood without taking into account the political and economic factors that surround it [27]. There are several approaches within this discipline depending on how the notion of power is applied to understand human/society/nature connections. In this research, the environment is regarded as "power-laden rather than politically inert" [28].

According to Zimmerer [29] methodology design and selection of methods should be carried out simultaneously, particularly with the selection of the theoretical framework and conceptual concerns that will be explored. In the case of social–ecological change, such methods must be compatible with the processes, spatial scales, and temporal frames that constitute the study objective.

For social–ecological change, Zimmerer (ibid) proposes spatial scales and temporal framing within discourse and content analysis. To delimit specific timeline stages in this study, two concepts related to environmental change were used as starting points for analysis: critical moment and episodic events.

The first concept has been applied in participatory nature conservation projects [30] and is described by Khan [6] as a:

"Conspicuous and sensitive moment that offers a specific insight into the interplay and autonomy of the actors involved in an event, one that illustrates or informs a political ecology analysis. The critical moment exposes the interplay and autonomy of both endogenous and exogenous actors" (p. 1)

The second concept is mentioned by Bailey and Bryant [7] in historical political ecology analysis to understand politicized environments in terms of physical changes, rate of impact, and political response; the episodic dimension is one that is often sudden and can be prolonged, impacting society in general but with unequal damage exposure, and the political response is focused on disaster relief.

Questions that led this study were:


A review of little-known documentary evidence and historical geographic information on the surface of the lakes, drainage infrastructure, and urban growth was conducted. The data used in this study include archival material from technical documentation about drainage in the Mexico Basin, as well as environmental and socio-political descriptions from each identified stage; the reviewed materials are mainly Spanish-language materials and were first published in Mexico.

In each historical period analyzed, different hydrotechnical works were carried out to prevent the flooding of Mexico City. The type of works built were based on the political moment, the purposes and interests of the actors, the economic resources or the state of emergency in the face of an environmental catastrophe. The impact of this was the reduction in the lake area until it dried up, which in turn allowed the expansion of Mexico City on the uncovered lands of the lake, while the deepest part of the lake became a federal area with exposed saline soils and small artificial lagoons that capture runoff from the basins to the east of the basin.

#### *Limitations of the Study*

Given the extension of the time period studied, information on some other factors that influenced the desiccation of Lake Texcoco may have been inadvertently omitted. Mexico has gone through several moments of transition and change, some of them violent social movements that could reduce the veracity and availability of the information at present; also, there was limited availability of data on the oldest events that were studied.

#### **4. Processes and Actors That Influence the Desiccation of Lake Texcoco**

#### *4.1. Tenochtitlan and Lake Texcoco (1325–1521)*

According to historians, the modern state of Mexico began in the late postclassic period, 625 years ago, with the founding of Tenochtitlan in 1325. Tenochtitlan's founders, the Mexica, were the last Mesoamerican people to build the city in the fourteenth century on one of the western islets of Lake Texcoco in the Valley of Mexico basin [25].

For two centuries, the Mexica altered the lake system, resulting in the development of intricate hydraulic infrastructure systems. The works served a variety of purposes and were constructed on both land and in riparian and lake environments. Palerm [24] identified and defined the most significant hydrotechnical works including *Chinampas* (described by the author as an "inside lagoon", an artificial island built by the Mexica for agriculture and settlement), causeways, *albarradones* (colonial and pre-Hispanic hydraulic works that functioned as dikes), flood defense and drainage works, freshwater conduction through canals, ditches and aqueducts, and formation of lagoons and artificial swamps.

At the middle of the island was the ceremonial center. The population settled around it on chinampas where they built their houses and developed agriculture, all this without damaging its natural lacustrine condition. The causeways that linked the island to the mainland in turn divided the lake. They also built aqueducts to supply drinking water to the island from the springs of Chapultepec and Coyoacán, navigable canals for the transport of people or goods and ditches for the distribution of drinking water and irrigation.

Tenochtitlan did not escape flooding. The community that evolved within that lake ecosystem figured out how to manage the lakes' water levels in order to limit damage. According to historical sources, two main constructions had their origin in environmental events, the Nezahualcóyotl and the Ahuízotl *albarradón* [31].

As a result of the 1449 flood, Nezahualcóyotl (*tlatoani* or emperor of Texcoco from 1402 to 1472) built a 16-km-long *albarradón* that connected the Tepeyac hill to the Santa Catalina Mountain range in order to prevent flooding caused by Lake Texcoco's overflow. This also allowed for the segregation of the Mexico Lagoon's fresh waters from the saline waters of Lake Texcoco.

On the other hand, the flood of 1499 was triggered by Ahuízotl (Tlatoani of Tenochtitlan from 1486 to 1502) when he attempted to build an aqueduct connecting Coyoacán and the city. The city was flooded for 40 days following the inauguration of the aqueduct and, in response to this catastrophe, they built another *albarradón*.

The Mexica developed a great empire in Mesoamerica through hydrotechnical works; by the 16th century, the city of Tenochtitlan had a population of between 170 thousand and 200 thousand inhabitants, and the entire basin had a population of about 1 million [2]. The prosperity of the cities inside the basin demonstrates that they used water as a tool to achieve a successful empire.

#### *4.2. The Spanish Colony and the Drain to Prevent Flooding (1521–1810)*

During the period of the Spanish conquest, there was a significant change in the landscape of Lake Texcoco. It is at this point in history that there is a transformation in the way people relate to the lake environment.

The Spanish explorers and conquistadors led by Hernán Cortés entered Tenochtitlan on 8 November 1519. After almost two years of invasion, the final strategy of the Spaniards was to cut off the supply of drinking water and the entry of food. The Mexica were defeated by using the same natural resource that they had used to establish their empire.

What occurred following the conquest would be significant for the future of Lake Texcoco. The capital of New Spain had to be built, and among the possible locations were Coyoacán, Tacuba, Texcoco, and Tenochtitlan. Hernán Cortés is credited with the idea to establish Mexico City on the ruins of Tenochtitlan, citing political measures to ensure the Spanish conquest was effective and definitive. In this regard, José María Marroquí in his book *La Ciudad de México* published in 1900 [32], mentions that Cortés considered it dangerous to leave the old city free given the possibility that over time the indigenous people could try to recover their temples, palaces and monuments.

Between the demolition of the Mexica city and the creation of the Spanish one, forty years passed. New emblematic buildings were built over the old ones such as the Cathedral of Mexico over Templo Mayor ruins. The pre-Hispanic period moved gradually towards the consolidation of Spanish rule and the beginning of the colonial period. Historians date the end of the conquest stage to around 1560, when tribute began to be expressed in pesos and paid in currency [33].

The new layout of the city broke the natural balance with the lake environment that the Mexica maintained. Furthermore, forests were cut down to use wood for new buildings, to make furniture or as fuel. Furthermore, some rivers were diverted, dried up or polluted. All this altered the physiognomy and productivity of the basin, starting the process of change in the hydrological system [34].

As geographer and naturalist Alexander von Humboldt pointed out in his 1827 publication Political Essay on New Spain [18], "the first conquerors wanted the beautiful valley of Tenochtitlan to resemble in everything the Castilian soil, in the arid and devoid of its vegetation" (translation from Spanish). His observation tells us that, unlike the Mexica civilization, the Spaniards were unfamiliar with floodplains as inhabited spaces, and therefore the conquest altered the interaction between humanity and the basin's natural lake environment.

The changes they made to the pre-Hispanic hydraulic system caused major floods, prompting a series of decisions to mitigate them, including attempts to replicate some of the works built by indigenous people in the past, but without achieving the same results due to the change already present in the environment.

An example of this is the *albarradón* of Nezahualcóyotl, which was demolished to allow the entry of Spanish vessels, causing a great flood in 1555. The Viceroy of New Spain, Luis de Velasco y Ruz de Alarcón, ordered the construction of a dike between the Calzada de Guadalupe to the north and the Calzada de Iztapalapa to the south of the city to substitute for the Netzahualcóyotl *albarradón*. This construction would be known as the *albarradón* de San Lázaro [19].

However, the *albarradón* of San Lázaro was not enough to stop the floods. In 1580, another great flood occurred. At that time, Viceroy Martin Enrique de Almanza was overseeing New Spain, and the idea of draining the lagoon surrounding Mexico City started to be mentioned as a solution to flooding.

Another flood occurred in 1604 and another shortly after, in 1607. Viceroy Luis de Velasco y Castilla ordered the construction of a tunnel that would pass through the northeast corner of the basin with the purpose of draining Lake Zumpango and preventing further flooding. The design would divert the waters of the Cuautitlán River to the Tula River. The person in charge of the design and construction of the Huehuetoca tunnel was the engineer and cosmographer Enrico Martínez [35]. The Valley of Mexico basin ceased to be an endorheic basin as a result of the development of this project.

Perló Cohen and González Reynoso [36] point out that the Huehuetoca tunnel prevented Mexico City from flooding as a consequence of the north's growing rivers and lakes. However, the east, south, and central valley overflows into Lake Texcoco were not controlled, and flooding occurred again in 1615 and 1623.

The flood of 1623 occurred during the government of Viceroy Diego Carrillo de Mendoza y Pimentel, who, according to Lattes [19], had opted to allow the lagoons to once again receive the flow of the rivers due to the high cost of maintaining the works. The entrance of the artificial drainage was blocked so that the waters of Zumpango could reach Lake Texcoco; the result was a gradual increase in water levels in the city and, by 1627, many of the streets were again flooded.

Numerous repairs were required both within the city and on the dikes and canals. The Huehuetoca tunnel was modified to become an open pit, allowing the necessary effluent to flow, but the Nochistongo gorge would not be completed until 1789 [35].

However, the most devastating flooding of the time had not yet occurred. In 1629, a deluge known as the San Mateo Deluge struck Mexico City. It rained continuously for 40 h. For five years, the city was submerged in water due to the flood. Water levels in the city center were recorded as high as two meters, and it would take until 1634 for a drought to put an end to such an overwhelming flood. The cost of the destruction was enormous. At the time, the population was believed to be 150,000, of which about 30,000 died as a result of the flood and its effects, such as the city's insalubrity. Another significant portion of the population emigrated to various cities throughout the country.

Throughout the 18th century, the government of New Spain continued to work on the drainage of the Valley of Mexico basin. Among them were the *albarradón* de San Mateo canal in 1747, the Guadalupe canal in 1796, and the Zumpango canal in 1798 to connect them to the Nochistongo gorge; despite the large lake area drained as a result of these works, the city continued to flood.

Valek Valdés [35] recounts the observations of the renowned geographer and naturalist Alexander von Humboldt on the works that supported the security of the capital at the beginning of the 19th century; among them were the stone causeways that prevented the waters of Zumpango from reaching the Lake San Cristóbal and ending in Lake Texcoco, the causeways and locks of Tláhuac and Mexicaltzingo, in the Nochistongo gorge that took the waters from Cuautitlán River to Tula River and, finally, the two Mier canals that allowed lakes Zumpango and San Cristóbal to drain.

As of now, the administration of natural resources was influenced by the omission of efforts that the country's native inhabitants had made for generations, shifting the concept of a civilization to adapt to the environment to now adapting the environment to the goals of the new civilization. According to Mann [26] in his analysis of the arrival of Spanish in the Americas, the assumption was that the new continent was pristine, and natural adaption mechanisms ignored.

#### *4.3. The Drainage Proceeds during the Independent Period (1810–1910)*

The country's political instability at the turn of the nineteenth century was evident in the condition of decay of the drainage works. During this time period, the government transitioned from a monarchy to a centralist republic to a federalist republic; hence, the drain was abandoned amid armed revolutions and power struggles. According to the engineer Francisco de Garay [3], it was forgotten at those times that personnel guarded the Nochistongo gorge in the north.

Dr. José María Mora, one of Mexico's most illustrious liberalism figures, was commissioned in 1823 to conduct a diagnosis on the state of the drainage work, in which he describes the state of complete abandonment and the pressing need for repairs. Mexico City was in danger of flooding due to the unevenness in relation to the lake caused by the constant deposit of earth carried by the rivers that flowed into the lake. Additionally, he expresses his opinion on the mishandling of past works, which were pricey and brief in duration [37].

In 1847, as the North American army marched towards the Valley of Mexico and Antonio López de Santa Anta was president of the Republic, the eastern sector was flooded as a defensive measure, digging the Viga canal and breaching the Mexicaltzingo locks. This move failed as a military tactic and instead served to exacerbate the risk of flooding in the capital and southern villages. Following that, engineers Francisco de Garay and M.L. Smith repaired the structure.

On 22 April 1853, the Ministry of Development was established, having authority over the canals and drainage of the Valley of Mexico. The proposal to use the city's hydraulic system for navigation was made as early as 1830 by Lucas Alamán, the country's Minister of Interior and Foreign Affairs Relations at the time, and was implemented briefly, but the shallowness of the canals and dikes, as well as the maintenance costs associated with them, increased interest in drainage. [4].

It was not until 1856 that the waters of Lake Texcoco threatened the city again, and in accordance with the established plan of action, a thirty-member General Board was appointed, which would later appoint a Minor Board to oversee the necessary repairs to the city's aging hydrotechnical works, with engineer Manuel Gargollo assigned to the north section, Juan Manuel de Bustillo to the center, and Francisco Garay to the south.

To the north, the Tepotzotlán and San Ignacio outlets were restored to the west bank of the Cuautitlán River, canals were cleaned to restore water to Lake Zumpango, and the Guadalupe canal was reopened. The San Cristóbal dike was reinforced in the center, and it was proposed that the Teotihuacan dike be restored in order to unload the Texcoco reservoir. It was, nevertheless, completely planted and occupied. The landowners presented the solution to facilitate the deposit of the waters by erecting many dams. All of these efforts were insufficient due to their unplanned nature, but in the south, the engineer Garay's efforts were part of a scheme presented years previously that was more successful.

In the south, the diversion in the Mexicaltzingo dike was closed, and the breaches in the Culhuacán and Tláhuac causeways were fixed, while the walls' height was increased to compensate for being totally submerged by the lakes' water. Additionally, the San Lorenzo canal was constructed to cut through the dividing hill between Xochimilco and Texcoco, and a floodgate was constructed in Mexicaltzingo to simply and rapidly manage water

flow. Over the next five months, the Chalco and Xochimilco lakes increased their water levels by 56 cm, while the Texcoco lake decreased.

"The lake of Texcoco being the lowest in the entire valley and with no outlet for its waters, all the salts that the landslides of the mountains pull on their route are deposited there," said Francisco de Garay [3]. Additionally, the salting in lakes was becoming more evident because in the dry season, during the ebb of its waters, salts sprout from the ground, rising to the surface by the capillarity effect or evaporation. The sun evaporates the water, leaving the salt behind, causing the deposit to grow (ibid).

De Garay also mentions that in the time prior to the Spanish conquest, it was common to find a large quantity of fish in Lake Texcoco, but that it had disappeared for at least a century. While the *albarradón* de Nezahualcóyotl was standing, the rivers of the west flowed around the capital, and through the San Lázaro floodgate, there was a freshwater stream at the foot of Peñón Chico, which turned it into a beautiful orchard. With less water coming from the south and no water coming from the north, an important area lost its ability to grow food. Instead, tequezquite, a mineral salt, or ground salt, was harvested.

During his exploratory trips to the American continent in 1803–1804, Alexander von Humboldt [18] describes the distance from the center of Mexico City to Lake Texcoco as 4500 m, and more than 9000 m to Lake Chalco, and notes a depth of 3 to 5 m in some areas, and up to a meter in others.

Thus, during the first two-thirds of the nineteenth century, the valley's drainage and the drying out of Lake Texcoco were effectively halted, but once a degree of stability was achieved toward the end of the century, the drainage was once again promoted primarily by political actors and social elites, despite warnings about the potential environmental consequences that were already looming. As Miranda Pacheco [38] points out, even though drainage works were restarted in the 1870s, the city's long-term growth was due to the "interests and projects" that had been built on private property and public land since the beginning of the century, even before the country became independent (59 p.).

On the other hand, Vitz [23] notes that the environmental practices that dominated the twentieth century were patrimonial and commercial exploitation of natural resources, which benefited from political instability and a colonialist-influenced society. This occurred in a political climate in which decisions affecting the natural environment were made on the basis of what was urgent and convenient. He also emphasized the differences in how different socioeconomic classes relate to the environment.

4.3.1. The Decision That Perpetuated the Desiccation of Lake Texcoco in the Second Mexican Empire

For the first time during the Second Mexican Empire, drainage of the entire valley was considered to combat flooding in the city. Despite Maximilian of Habsburg's liberal thinking and appreciation for Mexican nature, the draining decision was influenced by the conservative group that participated in the monarchy's formation. Miranda Pacheco [38] highlights the participation of politicians, doctors, lawyers, engineers, geographers and the upper class of society.

Shortly after Maximilian of Habsburg's arrival in 1864, the city faced another flood, and to find a solution, national and international experts were invited. The board of experts chose engineer Francisco de Garay's proposal, which included a network of staggered canals for drainage, navigation, and irrigation, with north and south extensions. The objective was to be able to control the level of nearly all vessels by creating additional reservoirs and eradicating stagnant waters. It was a very ambitious project that required a significant expenditure, which is why it was harmed by budgetary constraints following the restoration of the Republic in 1867.

Maximilian chose as his first option the drainage plan proposed by engineer M.L. Smith in 1848, who recommended not drying the lakes but connecting them via a canal, where it would be possible to balance the lakes' waters by utilizing the difference in their heights, thereby avoiding overflow into the city. Others proposed the construction

of artificial containers to repurpose excess water in agriculture. The order of 27 April 1866 reflected this decision; nevertheless, on 13 November the same year, the decree was published to carry out Francisco de Garay's plan, originally offered in 1856 and selected as the winner of the Minor Drainage Board's call for proposals.

This change of decision on the part of Maximilian was partly influenced by his advisers, who convinced him of the obstacles that the waters would present, overflowing or not, to his urban projects, as well as for those people close to him [38]. A letter written on 3 July 1866, by the First Engineer, Miguel Iglesias of the Ministry of Development, explains a different situation to Maximilian:

"Indeed, for our lovely capital, it is a benefit of such transcendence (the drain) that it may be stated that its life is in it, because without it, we will always have it nearly flooded and in risk of becoming completely flooded, and thus seeing it in ruins. It is also an issue for the other populations that dwell in its shadow: how much would their inhabitants lose if that center of business and consumption provided all of their needs? Additionally, even if it were not to protect them from the calamity of a flood, the welfare and health of the valley's inhabitants; the increase in fertility and production that will result from the cultivation of the vast lands currently occupied by the waters of one of the lakes that are so harmful due to the poor conditions in which they are found; the new and convenient communication routes created by the numerous canals that will cross the valley; all of these are enormous benefits" [4] (translation from Spanish)

Notably, M.L. Smith had previously emphasized the dry fields of the San Lázaro canal, an ancient link between Texcoco and Mexico City. He designated them as sterile plains due to their location within a salty lake basin. The flood of 1865 compelled them to act immediately, and as a result, a series of steps for the drainage works were undertaken. Francisco de Garay was appointed "sole and responsible director and inspector of all activities relating to the water issue in the Valley of Mexico" during this conflict. Francisco de Garay, on the other hand, was a devout Republican liberal who accepted the designation without title or compensation.

The lake's water had already reached the streets of Mexico City by October 1865. The situation was critical, and the city's engineers were called in to help. Francisco de Garay [3] indicated that the flood exceeded twice the usual surface of the lakes, invading the entire eastern part of the city, and that "to save the capital from danger, it was necessary to lower the level of Lake Texcoco". The urgent drainage plan presented by Francisco de Garay obtained most of the votes and was defended against Maximilian when he questioned its effectiveness. De Garay was so sure that he bet his own head [3].

The proposal called for isolating Lake Texcoco so that the northern and southern lakes would not continue to discharge their surpluses into the basin's lower reaches. Dikes were built, and help was once again sought from ranches, farms, and adjoining areas in order to construct artificial basins and reduce the risk of flooding in the city's towns and neighborhoods. The water level in Mexico City's historic center was 54 cm lower than in 1629. When the waters receded in 1866, the floors of the flooded streets and squares were ordered to be raised.

In May 1866, Emperor Maximiliano established the Drainage Commission, which included engineers Miguel Iglesias, Aurelio Almazán, Manuel Álvarez, and Jess Manzano, to conduct field studies and verify that the drainage project proposed by Smith was feasible, which they endorsed; work began at the end of the year. He also arranged for Miguel Iglesias to travel to Europe to purchase the necessary machinery [4].

The works conducted during this time period are detailed in the 1868 Development Memory and in the compilation of reports from the Historical, Technical and Administrative Memory of the Valley of Mexico Drainage 1449–1900; nevertheless, the lack of resources meant that most of the works of the drainage were postponed, and government engineers were then mainly concerned with conducting field studies and conservation works.

Maximilian of Habsburg was a European liberal who, upon his arrival, developed a policy in disagreement with the traditional position of the conservative class and the Mexican clergy but, due to his short duration in government, his influence did not reach the entire national territory.

Mexico City's architectural design, which began at this time, would last for decades, and projects for lake drainage would also continue. Vitz [23] points out that during this time, the priority was to safeguard the city from floods and protect private property, even when the interdependence of these with their environment was recognized.

#### 4.3.2. The Porfiriato (1876–1911) and the Basin's Second Opening

The hydraulic work that offers a second outlet to the basin is known as the Grand Drainage Canal. This work was part of several urbanization initiatives in Mexico City under the presidency of Porfirio Díaz, with the distinction that it was begun during Maximilian's reign. During the Porfiriato, Mexico experienced rapid development, which was made possible by foreign investment. The infrastructure of the railway, telegraph, telephone, and electricity were expanded. Despite the city's economic expansion and modernization, it was a time marked by economic disparities among the inhabitants.

The drainage works were cut off from the Ministry of Development in 1867 by the Ejército de Oriente in Mexico City. It was left behind because it was considered as an imperial endeavor. Andres Almazán, José Iglesias, and Jose Manzano came to Porfirio Díaz to ask for help. Porfirio Díaz's reaction, in his letter of 11 May 1867, is as follows [4]:

"Having read with pleasure VV's report on the work being done in Zumpango to aid Valley Drainage [..], I could wish for a few glories in my temporary position, such as adding momentum to these operations [...] As a result, [...] I have directed that the Federal District Treasury Department minister to them the sum of 1500 pesos each month for the conservation of drainage works, while the Supreme Government determines that they continue and carry out with appropriate diligence" (translation from Spanish)

Blas Balcárcel was heading of the Ministry of Development when the Republic was restored in 1867. Together with the Minister of Finance, José Mara Iglesias, he was successful in obtaining approval for a special tax to fund the drainage works. This tax included a 50% rise in municipal taxes collected at customs in the capital, as well as a 20% increase in direct taxes collected in the Valley of Mexico.

The drainage works were restarted in April 1868, but a later decree prohibited the special funds, and these entered a common fund. As a result, the works received an allocation from the general budget. The instability in the public peace influenced the perception of the amount destined for the works, which caused them to be paralyzed in October 1871. The completion of the outlet was the focus of the work done between 1868 and 1871.

The valley's drainage was also viewed as a solution to the city's health and hygiene issues. The Mexican Society of Geography and Statistics issued an opinion on the valley's drainage in 1874, arguing that doing so would prevent the physical and moral degeneration of their race and that the desiccated lands could be sold to benefit commerce, agriculture, hygiene, and the public treasury [38].

Miranda Pacheco [38] reveals the many perspectives on the valley's draining at the time. The Hygiene Commission of the Academy of Medicine in Mexico (at first), Doctor Eduardo Liceaga (president of the Medical Congress, family doctor, and close friend of Porfirio Daz), and Francisco de Garay were among those who supported it.

Those who opposed the drainage included Manuel Orozco y Berra, who believed that the city's drainage should be separated from the valley's drainage so that the Texcoco lake did not become a gigantic sewer. To stop the putrid emanations towards the city and limit the surface of the basin, Leopoldo Río de la Loza recommended reforestation of the area near Lake Texcoco. The doctor José G. Lobato noted, in his 1877 study, Hygienic Study on the Exanthematous Typhus, that preserving the surface of the lakes would reduce outbreaks of typhoid infections, a phenomenon that in pre-Hispanic times did not occur, but the Matlazahuatl pandemic happened when the lake's surface had dropped with the Spanish [38].

Rivalries between drainage officials, particularly between engineers Francisco de Garay and Luis Espinosa, were not new to the project. The tunnel's placement and hydraulic flow were two of the most contentious issues. Initially, the tunnel was designed to follow M.L. Smith's plan via the Acatlán ravine, and construction began in that direction. Francisco de Garay, on the other hand, thought it should go through the Ametlac ravine. When Manuel Fernández de Leal took over the Ministry of Development and designated him acting director of Works, Luis Espinosa's decision was eventually imposed.

As a result, until 1879, it was unclear how much water had to be extracted. The engineers predicted sizes and slopes for unusual and severe rain, which raised the cost of the project. Previously, M.L. Smith estimated a flow of 8 cubic meters per second, Francisco de Garay calculated a flow of 33 cubic meters per second, and Miguel Iglesias (head of the commission that studied and laid out the drain in 1866) calculated a flow of 41 cubic meters per second [4].

Porfirio Díaz accepted the drainage scheme for the Valley of Mexico presented by Luis Espinosa in 1879. It was based on the previous one by Francisco de Garay and took into account the work already completed by Miguel Iglesias but rectified the tunnel's hydraulic flow by 21 cubic meters per second. From 1871 to 1886, Luis Espinosa served as director of drainage works in the Valley of Mexico and as director of the Drainage Board of Directors.

According to Luis Espinosa [39], the lake of Texcoco was destined to disappear since the waters that came from the springs, wells, and drainage of the city were not enough to maintain the extension of the lake. Furthermore, he mentions the health issue commenting that "there are respectable opinions that foresee a danger to health in the complete drying up of the lakes".

With the return of Porfirio Díaz to power in 1884, the project gained enough momentum to be completed. The entire project consisted of three major parts: an open gorge 39.5 km long from Lake Texcoco to the northeast end of Lake Zumpango, a tunnel almost 10 km long from the shore of Lake Zumpango traveling northwest up to the Acatlán ravine, also known as Tequixquiac, and a canal. This would carry water from the Tequixquiac River, which meets the Tula River and, downstream, creates the Moctezuma River, to Pánuco, and finally empty into the Gulf of Mexico [4].

During the Porfirio Díaz administration, a budget close to 16 million pesos was allotted for the construction of the Grand Drainage Canal between 1886 and 1900. It was dedicated on March 17, 1900, in the presence of the President of the Republic, some Secretaries of State, various members of the diplomatic corps, the Drainage Canal Board of Directors, engineers, employees, and visitors from banking, commerce, industry, and the arts. In this regard, Luis Gonzales Obregón, a member of the Drainage Canal's Board of Directors, wrote:

"March 17, 1900, will be a memorable date, because the works inaugurated on this day together with those of sanitation, will make Mexico one of the most pleasant mausions, among the capitals of the American Republics, for its beauty, health and climate" (translation from Spanish)

Mexico City was already far removed from its lacustrine origins at the time of the Grand Drainage Canal's opening. The Hydrographic Commission of the Valley of Mexico estimated that the lakes' maximum extension during the rainy season would be 536.50 square kilometers, with the surface area of Lake Texcoco ranging between 183.28 and 238.54 square kilometers, depending on the amount of water received and the amount of precipitation expected throughout the year (Table 1).


**Table 1.** Lake surface on the Valley of Mexico in 1860—1870.

Source: Luis Espinosa, Technical and administrative report on the Valle de México drain, 1902, in Legorreta, Jorge (2006). Estimates made with data from the Valle de México Hydrographic Chart, by Manuel Orozco y Berra.

In the act issued on the inauguration of the drainage works in the Valley of Mexico, the Grand Drainage Canal is described as a great and beneficial work that will liberate Mexico from floods and improve the hygienic conditions of the capital and the valley [4]. The work was seen as an engineering feat and a symbol of the progress and modernity of the country.

However, only four months later, a second flood revealed that the work was ineffective. The Grand Canal's status would deteriorate over time, as the canal's movement was determined by gravity and Mexico City was sinking. The Grand Drainage Canal displaced about one million cubic meters of water in its first four years [40].

Although the floods continued, this project was critical for future urban expansion because it created new territory that could be salvaged by draining the waters of Lake Texcoco. The decrease in Lake Texcoco to its current extent and the urban area of Mexico City are depicted in Figure 2.

#### *4.4. Drainage of Mexico City in the Twentieth and Twenty-First Centuries: The Battle against the Lake Is Not Over*

Efforts undertaken centuries ago to minimize flooding in Mexico were insufficient. The rapid growth of the city, constructed on a historic lake complex, prompted the development of increasingly complicated infrastructure. This would eventually impact the city's hydraulic infrastructure, necessitating ongoing drainage projects and repair of current ones. Reyes [41] writes about it:

"It covers the drying up of the valley from the year 1449 to the year 1900. Three races have worked in it, and almost three civilizations —that there is little in common between the viceregal organism and the prodigious political fiction that gave us thirty years of august peace. Three monarchical regimes, divided by parentheses of anarchy, are here an example of how the work of the State grows and corrects itself, in the face of the same threats from nature and the same earth to dig. The slogan of drying the earth appears to have run from Nezahualcóyotl to the second Luis de Velasco, and from him to Porfirio Díaz. Our century found us still throwing the last shovel and digging the last ditch"

(translation from Spanish)

The first Tequixquiac tunnel, inaugurated in 1900 by Porfirio Díaz, had a flow of 16 m 37 s; this was enough for Mexico City, which at that time had around 500,000 inhabitants and when the rivers that crossed it emptied into Lake Texcoco. Subsequently, the water discarded by the drainage system began to be used for irrigation of agricultural land, forming part of the irrigation district of Tula, Hidalgo in the Mezquital Valley region.

Only 37 years after the tunnel's completion, it became insufficient to accommodate Mexico City's population increase. By 1940, the city had grown to 1,757,530 inhabitants, and fearful of the first tunnel collapsing, the federal authorities decided to expand the system. Drainage was undertaken via a second tunnel through Ametlac ravine in Tequixquiac. The new tunnel would be the basin's third artificial outlet.

**Figure 2.** Hydrotechnical works that open Mexico Basin to Tula Basin. Own elaboration based on data from CONABIO [42] and Stangl [43], historical maps recovered from Mapoteca Orozco y Berra [44], and photo-interpretation of satellite images from 2021.

The second tunnel began construction on 1 July 1937, initially under the supervision of the Ministry of Communications and Public Works and later transferred to the Ministry of Hydraulic Resources. It was opened in 1946 by President Manuel Avila Camacho and had a flow rate of 60 cubic meters per second, a diameter of 4 m, a length of 11.3 km, and a slope of two thousandths.

However, there was a tremendous ecological cost. These three exits contributed significantly to drainage of the lakes and rivers, hence reducing the risk of floods in the city. The loss of the natural balance of the hydraulic system occurred as a result of the artificial outlets. Since 1608, the drainage flow has increased without being matched by an equal amount of incoming water, compressing the subsoil clays and resulting in differential subsidence in the soil.

The sinking has impacted not only Mexico City, but also other areas of the old lake of Texcoco, with one of the repercussions being the shutdown of the drainage system. When the Grand Drainage Canal was launched, the city's base was 5 m higher than it is now, the collectors had a sufficient slope, and the flow could be dragged to the Tequixquiac tunnel, but by 1950, the situation had changed. Water began to withdraw from the collectors and flood the streets as a result of the unevenness.

On 16 July 1951, the city center was flooded for months due to the drainage system's inefficiency in removing rainwater. The water level reached two meters, and residents had to rely on boats to move around. The Hydrological Commission of the Valley of Mexico Basin recommended installing a system of pumps in the collectors to transport the water to the Grand Drainage Canal, a solution that would be costly due to the amount of electricity consumed.

The interceptor and west emitter were built in 1960 to receive and dislodge waters from the Nochistongo gorge. However, the city's rapid and steady growth rendered it insufficient. Lake Texcoco was 2 m below the city center in 1900, but 5.5 m above it in 1970 [36]. As a result, the federal government (Department of the Federal District) concluded, based on the findings of the General Directorate of Hydraulic Works investigations, that a drainage system was required that was not affected by subsidence and did not require the use of pumps, thereby lowering the cost of its operation. As the Tequixquiac tunnels were full, it was decided to build a fourth artificial outflow from the basin.

The first stage of the Deep Drainage System was built in 1967, during the presidency of Gustavo Díaz Ordaz, and was opened in 1975, with the government of Luis Echeverría. The first stage involved the construction of a 60 km tunnel with a 200 m depth and a discharge rate of up to 200 cubic meters per second. Luis Echeverría called this one of the century's most important works, one that would finally free the city from flooding [36].

An extension of 5.5 km to the central interceptor and a new interceptor of 16 km commenced during José López Portillo's six-year administration (1976–1982). The Deep Drainage System has been extended to more distant parts of the city, as has the Grand Drainage Canal, thus using existing infrastructure and linking the entire city. Under Carlos Salinas de Gortari (1988–1994), the Deep Drainage System covered 125 km.

For a long time, the drainage system was solely meant to remove rainfall and was closed annually for repair. In the mid-1990s, a combination system that dislodges both rainfall and sewage began to be deployed. Drainage work has persisted throughout the first two decades of the twenty-first century since the city continues to have flooding issues connected to expansion and subsidence.

The Deep Drainage System was expanded under Felipe Calderón Hinojosa's presidency. The East Emitter Tunnel (EET) was built in 2008 to increase system flexibility. In the municipality of Atotonilco de Tula in the state of Hidalgo, near the exit of the central emitter, the EET discharges up to 150 cubic meters per second. It is 62 km long, 7.5 m wide, and 55 to 150 m deep. A decade late and 20 million pesos over budget, it began operations in 2020. However, it is one of the world's largest drainage works.

In 2014, under Enrique Peña Nieto's administration, the first stage of the West II Emitter Tunnel was completed, reducing river overflows near Naucalpan, Tlalnepantla, Cuautitlán Izcalli, and Atizapán. It is 5573 km long, 12 to 110 m deep, and has a drainage capacity of 112 cubic meters per second. This project was nominated for Work of the Year in 2019 (see Figure 2).

Lake Texcoco went from covering almost the entire basin, about 7868 square kilometers at its origin [1], to 272.17 square kilometers in the mid-nineteenth century [40] and at present,

the remnants of the lake cover an area of only 16.83 square kilometers within the federal area (see Figure 3).

**Figure 3.** Historical evolution of the desiccation of Lake Texcoco. Own elaboration based on data from CONABIO [42] and Stangl [43], historical maps recovered from Mapoteca Orozco y Berra [44], and photo interpretation of satellite images from 2021.

#### *4.5. Recent Climate Change in Mexico Basin*

Mexico City's paradox is that it has built an artificial mega basin by importing water from the Lerma and Balsas basins to serve the city, while also expelling the waters through the constructed drainage system into the Tula basin (see Figure 4). Additionally, the extraction of water from the subsoil has resulted in the city sinking lower than the lake in respect to the difference from the place where it was located in pre-Hispanic times, leading to the burying of the city drainage deeper and deeper [45].

The delivery record of Cutzamala system (in the Balsas river basin) to Mexico Basin shows that, for the period 1997–2008, the average import was 15,162 cubic meters per second, while the corresponding record of the Lerma system shows that the respective import supply to the basin was 4231 cubic meters per second [46].

The ejected water is directed to irrigation district 003 in Tula, Hidalgo, which covers approximately 50,000 hectares and produces mostly corn (*Zea mays*) and alfalfa (*Medicago sativa*) [47] and, for the years 1989–1992, the yearly amount of exports to the Tula basin was 52,782 cubic meters per second [46]. It should be emphasized that Mexico City's and the metropolitan area's drainage system is of the combined sort (wastewater and rainwater), and that drainage expansion has continued in recent years, as has its water expulsion capacity.

**Figure 4.** Megabasin created by alterations in the Basin of Mexico throughout history. Own elaboration based on data from CONABIO [42].

The Intergovernmental Panel on Climate Change (IPCC) anticipates that both environmental humidity and global temperature will increase [48]. Mexico City's transformation as a result of human-induced changes has resulted in an increase in the frequency and intensity of extreme events such as floods, droughts, and heat waves. According to Romero Lankao [49] between 1980 and 2006, there were 668 floods, 60 storms, 53 showers, 60 hail storms, 85 fires, and 14 droughts.

The National Water Commission's (CONAGUA) records on annual mean precipitation and annual mean temperature over the previous 37 years [50] indicate a decline in the amount of rainfall occurring in Mexico City annually and, in turn, an increase can be observed in the mean annual temperature (Figure 5).

Between 1920 and 1940, the Basin of Mexico had substantial population increase, particularly after 1940, and peripheral areas have been absorbed into the megalopolis. The forests and agricultural fields have vanished entirely, while there still exist portions of the lakes in the basin's south, they are endangered due to the pressures of urban water requirements.

According to Ezcurra et al. [45], numerous aquatic, subaquatic, and halophilic species have become extinct, while other types of vegetation are on the edge of extinction. Forest communities have been deeply affected in general; since 1985 over 9000 hectares of forest have been lost, and those that remain are threatened by pollution and pests. On the other hand, the introduction of exotic species to the basin's forests has resulted in the extinction of animal species that rely on indigenous plants for food and shelter. Eucalyptus trees, in particular, were planted extensively throughout the basin to increase evaporation and accelerate the lacustrine system's desiccation process.

**Figure 5.** Graph of historical data on annual mean precipitation and annual mean temperature of Mexico City.

A study of changes in annual mean precipitation and temperature from 1950 to 2013 [51] indicated that average winter temperature and summer precipitation increased at an average rate of 0.1 degree Celsius per decade and 17.8 mm per decade, respectively. They forecast a temperature increase of between 1 and 3 degrees Celsius over the next 30 years, while rainfall is expected to increase and become more severe at the conclusion of the rainy season.

Clearly, the lake was abandoned in favor of urbanization. Following the Mexican Revolution, and particularly under Lázaro Cárdenas's (1934–1940) land reform, the desiccated territory surrounding Lake Texcoco was targeted for agriculture. They did not realize, however, that the fertility of the flooded land was attributable to the lake's presence.

Additionally, considerable efforts were made to restore what remains on the dry bottom of the former Lake Texcoco. CONAGUA initiated a program in the 1970s to plant halophytes on the lake's ancient bed in order to mitigate dust storms that afflicted the eastern part of Mexico's basin during the dry season. To retain a portion of the lake's ecological functions, a series of artificial lagoons with a combined surface area of 17 square kilometers were constructed. However, the government lost interest in this program, and for the last two decades, there has been a lot of debate about what to do with that area.

#### **5. Discussion and Conclusions**

For this study, segmenting the lake's history into periods marked by a critical moment in the country's social dimension enables the delineation of environmental policies that had a direct impact on Lake Texcoco, allowing for the tracking of changes in the use of this natural resource over time.

The indigenous Mexica were able to control the levels of the lakes by building dikes, roads, gates, and viaducts, which enabled them to live peacefully while taking advantage of the resources of the lake ecosystem through hunting, fishing, agriculture in chinampas, and transportation via the hydraulic system they had developed. While they utilized forests, altered and controlled the level of lakes, diverted rivers, and filled lakes in order to gain territory, their consequences were not as severe as those of the subsequent populations.

The pre-Hispanic era's key actors are the leaders of cities who assume the role of emperor; what stands out in Mexica society's urban design is the influence of their religious beliefs on their way of life, from the location of the city's foundation to the city's layout in relation to the aquatic environment, as Palerm [24] and Matos [2] indicate.

What occurred in Mexico during the Spanish colonization period is a reflection of the environmental approaches used in Europe at the time in countries such as Germany, Italy, France, and England [52], however few had lasting favorable effects.

The Spanish introduced new forms of agriculture, livestock, and forestry that had a negative impact on the natural landscape. The decision to build the capital of New Spain on the ruins of Tenochtitlan resulted in devastating floods and diseases among the local population, leading to the perception that the lakes were a danger.

Despite the fact that an attempt was made during the colonial period to use pre-Hispanic ways to manage the lakes, the transformation of the environment was so profound that they were insufficient, and the idea of drainage eventually won out over the other options. Enrico Martínez is best known for his design of the Huehuetoca tunnel and, later, the Nochistongo gorge, which resulted in the disruption of the basin's water balance and the creation of an outlet to the Tula River, respectively.

The tremendous floods that devastated the city at the colonial period prompted the viceregal government to make this choice, which resulted in the city being submerged under water for years. Making a snap decision in response to imminent natural events was a very risky move.

This analysis confirms Candiani [52], Vitz [23], and Tortolero's [21] assertion that colonization and the introduction of capitalism to Mexico resulted in the drying up of Lake Texcoco. This, combined with an unwillingness to live in a city surrounded by water, transformed the resource from a sanctuary to an unhealthy and dangerous environment for its residents.

As the nineteenth century progressed, the discussion over what to do with the lakes was still going on. During this period, three actors stand out who contributed to the speeding up of the drying process: the engineer Francisco de Garay, who made a number of interventions in the city's hydraulic system, the most notable of which was the design of the city's drainage system in 1856. This was authorized during Maximilian of Habsburg's empire in 1866 after being persuaded by Engineer Miguel Iglesias since this option would benefit, first and foremost, the interests of the group close to him.

Miranda Pacheco [38] and Tortolero [21] identify actors in this stage including doctors, engineers, specialist technicians, and members of Mexico City's top socioeconomic class during the Porfiriato era.

During the Porfiriato, the continuance of imperial works aimed to legitimize the state. It was a symbol of modernity and the power of the state that, far from providing a final solution to the city's drainage difficulties, was a determining element in the emergence of new environmental problems and the hunt for more technical solutions to floods in the years after its construction.

Deep drainage was required as a result of the collapse of the Porfirian drainage system, which has continued to spread to this day. History appears to be repeating itself, but in ever bigger and even more tragic proportions in light of the rapid population development and ecological deterioration of the basin in recent decades.

In modern Mexico, the state and businesspeople emerge as actors. The inclination of political actors, particularly in the post-revolutionary period, to exploit Lake Texcoco to construct a new political–social order is noticed; with Cárdenas, the use of natural resources developed in tandem with social change, what Boyer and Wakild [20] refer to as "Social Landscaping".

When Lake Texcoco dried up, it exposed the lands at its bottom, allowing the urban area to expand. However, in recent years, the remnants of Lake Texcoco have been put on the table for discussion over their future uses, which has sparked opposition from some members of the community. However, before making any further changes to the lake's ecosystem, it is critical to consider the cycle of historical decisions that have resulted in one environmental calamity after another.

Blaikie and Brookfield's study [9] is a seminal work in political ecology for explaining land deterioration. Using this as a starting point, a model was built to depict the cumulative decisions on Lake Texcoco that resulted in environmental degradation. The accumulation of hydraulic projects necessary to "keep the city safe" has degraded the ecosystem. Political economy is critical, since it controls the dynamics of changes in the country's development, which is directly mirrored in Mexico City and its environs.

For example, in post-revolutionary modern times, when floods continued to occur despite the completion of the big drainage canal, the choice to continue developing Mexico City's drainage system and expelling all water from the city, whether rainfall or sewage, prevailed.

Mexico City has experienced flooding continuously, and while it once possessed vast water resources, it now faces a high vulnerability to climate change and the inability to meet the city's water demand. The examination of decisions made on Lake Texcoco and the identification of individuals enables this to be attributed to environmental changes and changes in the hydrological cycle.

The change in the hydraulic cycle has been profound and irreversibly altered the ecosystem, resulting in a paradoxical situation in which floods continue despite the sophistication of the drainage system, while nearly a third of drinking water must be imported from increasingly greater distances. Thus, in addition to impacting the basin's climate, the alteration of the water system has impacted the regional water balance and availability.

Throughout the basin's environmental transition, misinformed decisions, combined with political, economic, and social influences, have had a detrimental impact on natural systems. Mexico City, like cities such as Dhaka, Beijing, and Sao Paulo, relies heavily on groundwater and interbasin water transfers for water supply [53]. Drying, drainage, and transfer projects are excellent examples of how social priorities functioned and evolved over time. These cities are particularly revealing about how humans interact with nature and affect ecology, land use, and water resources.

**Author Contributions:** Conceptualization, C.M.-R., E.O.-T. and V.E.-H.; data curation, C.M.-R. and E.O.-T.; investigation, C.M.-R.; methodology, C.M.-R. and E.O.-T.; Supervision, V.E.-H.; validation, D.F.-R., M.C.D. and G.S.B.V.; visualization, C.M.-R.; writing—original draft, C.M.-R.; writing—review and editing, C.M.-R. and E.O.-T. All authors have read and agreed to the published version of the manuscript.

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

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. Geographic information can be found here: CONABIO. 2022. National Biodiversity Information System. http: //www.conabio.gob.mx/informacion/gis/ (accessed on 14 March 2022). Climatological information can be found here: CONAGUA. 2022. Monthly Temperature and Rainfall Summaries: https://smn. conagua.gob.mx/es/climatologia/temperaturas-y-lluvias/resumenes-mensuales-de-temperaturasy-lluvias (accessed on 14 March 2022). Stangl, W., & (FWF), A. S. F. (2019). Data: Lago Texcoco, 18th century (V1 ed.). Harvard Dataverse. https://doi.org/doi/10.7910/DVN/DOT5EY (accessed on 14 March 2022).

**Acknowledgments:** The authors thank the National Science and Technology Council (CONACYT) for the scholarship awarded for the doctoral studies from which this research is derived.

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

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Land* Editorial Office E-mail: land@mdpi.com www.mdpi.com/journal/land

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-5646-8