**Water Security Assessment of Groundwater Quality in an Anthropized Rural Area from the Atlantic Forest Biome in Brazil**

**Igor Fellipe Batista Vieira <sup>1</sup> , Fernando Cartaxo Rolim Neto <sup>1</sup> , Marilda Nascimento Carvalho <sup>2</sup> , Anildo Monteiro Caldas <sup>1</sup> , Renata Cristina Araújo Costa <sup>3</sup> , Karolyne Santos da Silva <sup>4</sup> , Roberto da Boa Viagem Parahyba <sup>5</sup> , Fernando Antonio Leal Pacheco <sup>6</sup> , Luís Filipe Sanches Fernandes <sup>7</sup> and Teresa Cristina Tarlé Pissarra 3,\***


Received: 13 December 2019; Accepted: 20 February 2020; Published: 25 February 2020

**Abstract:** The exploitation of natural resources has grown mainly due to the high rate of population growth that changed over time around the planet. Water is one of the most needed resources essential for survival. Despite all the efforts made to improve water security, an environmental impact related to anthropogenic influence remains of great concern, which is the alteration of surface and groundwater quality. In many regions around the world, there is limited or no access to rural and urban water supply while there is a need to improve sanitation facilities. This work evaluated the spatial distribution of groundwater and surface water quality as well as their changes in wet and dry seasons of the tropical climate in the Atlantic Forest Biome. The study area is under anthropogenic influence, which is in the municipality of Igarassú, Pernambuco State, Brazil. The analysis of the raw water was based on Standard Methods for Examination of Water and Wastewater, as referenced in the Brazilian Ministry of Health Consolidation Ordinance that sets standards for drinking water. The temporal analyses indicated a variation on water quality from the wet to the dry seasons, whereas the spatial results revealed deviations from the Brazilian's Water Supply Standards for some physicochemical parameters. There was an increase in the values of some parameters during the wet season in some hydrological compartments. The anthropized rural area from the Atlantic Forest Biome is affecting the water quality. It is, therefore, necessary to develop environmental policies and put them into practice by implementing engineering projects that guarantee proper treatment for raw water in order to bring the water quality back to a good status in this region.

**Keywords:** environmental monitoring; water quality; surface water; groundwater; drinking water

### **1. Introduction**

The management of water resources as well as the sustainability of groundwater and surface water systems are topics of great concern [1–3]. Water covers 71% of the Earth's surface, but only 0.3% is available for drinking water [4]. The access to an adequate, reliable, and resilient quantity and quality of water for safe drinking is the main global water security issue related to aquatic ecosystems and human health [5–7] as well as to economic values. The understanding of natural processes and anthropogenic factors [8–10] and their roles for management and sustainability of water security resources require a better comprehension of changes that occur in groundwater and stream water quality [7,11], mainly to develop a governance and an implementation of water and land use policies [12].

Water governance is an excellent alternative for understanding and developing ways for water security [3,13] as well as seeking sustainable ways to exploit the groundwater resource to ensure human development [14] and integrated management [15]. Thus, environmental laws and guidelines are responsible for ensuring groundwater quality standards for consumption [16]. The Brazilian legislation and recommendations related to groundwater and stream quality list the parameters that must meet a certain potability standard directed to human consumption on drinking water [17].

Most of the water in the world goes to irrigation and agriculture, which is estimated at 70%, while industry uses 22%, and domestic use is at 8% [4]. In Brazil, the National Water Agency [18] through the Conjuncture of Water Resources in Brazil estimates that 72% of the country's water is destined for agriculture, 9% is destined for livestock, 6% is meant for industry, and, lastly, 10% is meant for domestic use. Land degradation is also a major source of water pollution related to erosion processes [11,12,19], contamination by heavy metals [20,21], eutrophication [22], and others, which, if unmanaged, can lead to significant economic and environmental costs.

Water quality is influenced by agricultural activities [23,24] and other land uses [11]. However, the consequences for water quality of some activities such as sand extraction or drilling of clandestine wells are still inconclusive or ambiguous. The extraction of sand through dredging is an important activity in the studied area with a need to supply the construction sector. Sand extraction can be an environmental stressor of surface and groundwater quality because the municipality of Igarassu explores a significant number of wells for the supply to local communities [25]. Therefore, the ecosystem integrity is vulnerable to physicochemical, erosive, suspended solid disturbances among others [26].

The Brazilian Northeast faces some inequalities in the access to water resources with certain population groups lacking a water supply system [25,26]. Extensive periods of drought can lead to irreversible socio-environmental impacts, which are related to soil water infiltration, increased runoff, and intensification of erosion. In the state of Pernambuco, the Pernambuco-Paraíba Basin is one of the largest underground water reserves. The Beberibe aquifer is one of the most important public water supplies, and the water company from the state of Pernambuco is responsible for distributing the drinking water to the municipalities of Recife and Igarassú. There are many socioeconomic activities that take advantage of water resources in the region. Additionally, the population in neighborhoods exploit groundwater [26,27].

Igarassú is a municipality of Pernambuco state, Brazil, with an area of 305,560 km<sup>2</sup> , located 28 km away from the capital Recife. The climate is defined as Group As (Tropical savanna climate or tropical wet and dry climate), according to the Köppen climate classification. There are two well-defined dry and wet periods [27]: the wet season begins in the fall (May to October) and the dry season starts in the summer (November to April). The average annual rainfall is 1634 mm.

In a developing municipality such as Igarassú, there are densely populated areas with insufficient water supply, incomplete or inexistent sanitation, and ineffective environmental planning [28]. Under these conditions, water can transmit waterborne diseases and degrade environmental quality [28–30]. It is, therefore, important to systematically identify the factors that lead to water quality deterioration and propose solutions that are likely to secure the path of sustainable development [26,30,31].

This study was evaluated through field visits and laboratory analyses, the spatial distribution of groundwater, and surface water quality as well as the changes occurring from the wet to the dry season in the studied area, which aims to identify environmental issues that interfere within the quality of water consumed by rural communities as well as to propose solutions to impacts caused by human action within the framework of water security. groundwater, and surface water quality as well as the changes occurring from the wet to the dry season in the studied area, which aims to identify environmental issues that interfere within the quality of water consumed by rural communities as well as to propose solutions to impacts caused by human action within the framework of water security.

*Water* **2020**, *12*, 623 3 of 18

This study was evaluated through field visits and laboratory analyses, the spatial distribution of

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

34°53'05.4" W (Figure 1).

The study area is in an anthropized rural area from the Atlantic Forest Biome of Brazil, State of Pernambuco, located in the municipality of Igarassú, with central coordinates 7◦51009.5" S and 34◦53005.4" W (Figure 1). **2. Materials and Methods** The study area is in an anthropized rural area from the Atlantic Forest Biome of Brazil, State of Pernambuco, located in the municipality of Igarassú, with central coordinates 7°51'09.5 "S and

**Figure 1.** The location of the study area is an anthropized rural area from the Atlantic Forest Biome in Brazil, State of Pernambuco, Igarassú Municipality. **Figure 1.** The location of the study area is an anthropized rural area from the Atlantic Forest Biome in Brazil, State of Pernambuco, Igarassú Municipality.

The study area covers approximately 16.1 km2 with the anthropogenic pressures concentrated in the Northeast portion of the municipality. According to the geological map of the municipality of Igarassú (http://rigeo.cprm.gov.br/xmlui/handle/doc/16272), the area belongs to the Borborema Province where rocks from the Salgadinho Complex, sediments from the Beberibe and Gramame formations, and fluvio-marine and alluvial deposits crop out. Among others, the lithologic types include conglomerate and clay siltstone, sandstone with calcareous cement, and phosphorite interdigitated with calcarenites [29]. Topography is characterized by an undulated relief whereas soils are mostly represented by Yellow Latosols and Yellow Argisols with a sandy texture [27,29]. The study area covers approximately 16.1 km<sup>2</sup> with the anthropogenic pressures concentrated in the Northeast portion of the municipality. According to the geological map of the municipality of Igarassú (http://rigeo.cprm.gov.br/xmlui/handle/doc/16272), the area belongs to the Borborema Province where rocks from the Salgadinho Complex, sediments from the Beberibe and Gramame formations, and fluvio-marine and alluvial deposits crop out. Among others, the lithologic types include conglomerate and clay siltstone, sandstone with calcareous cement, and phosphorite interdigitated with calcarenites [29]. Topography is characterized by an undulated relief whereas soils are mostly represented by Yellow Latosols and Yellow Argisols with a sandy texture [27,29].

The water security assessment of groundwater quality was analyzed in three hydrological compartments (HC), as defined in Figure 1. Raw water samples were collected within the HCs and numbered from SN1 to SN16 and geo-referenced (Table 1). The sampling locations were strategically selected, considering the distribution of rural communities and their population. The samples, taken altogether, were representative of the different sources from which drinking water is obtained by the public rural community. The water security assessment of groundwater quality was analyzed in three hydrological compartments (HC), as defined in Figure 1. Raw water samples were collected within the HCs and numbered from SN1 to SN16 and geo-referenced (Table 1). The sampling locations were strategically selected, considering the distribution of rural communities and their population. The samples, taken altogether, were representative of the different sources from which drinking water is obtained by the public rural community.

The sampling of groundwater was conducted in drilled wells from the water table aquifer (SNs 1, 2, 5, 11, 13, and 15) in cacimbas that are small excavations dug near the streams to reach the water table with the purpose to remove water for domestic use or small plantations, which can be lined


**Table 1.** Identification (number), source, and geographic coordinates of the water samples.

The sampling of groundwater was conducted in drilled wells from the water table aquifer (SNs 1, 2, 5, 11, 13, and 15) in cacimbas that are small excavations dug near the streams to reach the water table with the purpose to remove water for domestic use or small plantations, which can be lined with a concrete pipe to prevent the collapse of their walls (SNs 4, 6, 7, 8, 9, 10, 14, 16) in the stream (SN 3) and from a dredging pond, which means an area where silt and sand are removed from the bottom of the water bodies (SN 12). The water well depths ranged from 15 m to 128 m, and the cacimbas depths ranged from 2 m to 8 m. The hydrological compartments (HCs) were drawn using a set of 1/10,000 scale contoured orthophoto cards by considering the water divides.

Campaigns for field data collection and raw water sampling in the study area were carried out during the rainy season between the months of March to July, called Wet Season Samples, and during the season without rain from August to February, which is called the Dry Season Samples. The sampling was done by collecting raw water in a 2-L bottle preserved under refrigeration for analysis, according to Standard Methods for the Examination of Water and Wastewater [32]. The samples were taken to the Minerals, Soils, and Water Analysis Laboratory (LAMSA), located at the Department of Chemical Engineering of the Federal University of Pernambuco (UFPE), for processing and conducting physicochemical and microbiological analyzes. The physicochemical analyses were based on specific protocols and assumed as pre-defined standards, described in PRC No. 5 - Annex XX, Chapter III, section V, Art. 22 of the Brazilian Ministry of Health [33] and in American Public Health Association (APHA) [32]. The samples were collected during the day between 6:00 to 18:00 h at the same hour for each sample number at a given site, by following the numbers from SN1 to SN16, during seven consecutive days in each season. The stream (surface water) was sampled at about 7:30 a.m.

Raw water temperature, turbidity, pH, total dissolved solids, dissolved oxygen, ox-redox potential, and electric conductivity were measured in the field during the sampling campaigns using a multiparameter probe called the HORIBA model U-50 (Table 2).

The heavy metal concentrations (lead, copper, total chromium, zinc, and cadmium) were measured by atomic absorption spectrophotometry. UV-VIS spectrometry determined sulfates, nitrate, and nitrite. Flame photometry was used to determine sodium and potassium. Analyses of total hardness and alkalinity (Mohr method) were performed by volumetric analysis. Total iron, aluminum, ammonia, and color were determined using Merck Millipore spectrophotometric analysis, using Merck Millipore PHARO 100 (VIS) and 300 (UV-VIS) spectrophotometers. The method was based on the use of Millipore filters with pores of 0.6 micron. Microbiological analyses of total and thermotolerant (fecal) coliforms

and heterotrophic bacterial count were based on the methodology of the Standard Methods for the Examination of Water and Wastewater [32].


**Table 2.** Parameters measured by the Horiba U-50 model probe and their measurement units and precisions.

The analyses were based on Directives 98/83/EC and (EU) 2015/1787, which are methods accepted by PRC No. 5 - Annex XX Chapter III, section V, Art. 22 [33] and by the APHA methodology [32].

Geospatial distribution maps were drawn for the pH, color, turbidity, and total iron parameters within the study area. The Nephelometric Turbidity Units (NTU) varies from 0 to 800 NTU. The Hutchinson's Topo to Raster [34] interpolation method was used. This method allowed the use of contours, basin boundaries, and georeferencing to interpolate data, which highlights the areas where the quality of water is not conforming with the legal standards [33]. The interpolation was performed using the weighted sum of squares in the residuals by the surface elevation data [34,35]. Thus, through this interpolation model, vector data were converted into hydrological land models [34–36].

Radar charts were drawn to display the multivariate water sample observations. Each spoke in the chart represents one variable. The length of a spoke is proportional for the magnitude of the variable. Radar charts were drawn for every data point (SN1 to SN16) and were represented in a multi-plot format.

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

### *3.1. Water Parameters' Interaction at the Sample Number and Hydrological Compartments*

The specific studied area was divided in three hydrological compartments (HC) considering the altitude of land that drains the water downslope to the lowest point (HC1, HC2, and HC3) and the anthropized rural area from the Atlantic Forest Biome in Brazil.

Identification of seasonal trends (wet and dry) in raw water-quality constituents was especially important because high or low rates of each parameter have such a substantial effect on analyses of an anthropized rural area [9]. Surrogates monitored on a continuous basis provide resource managers with real-time information on sample number properties that showed different water characteristics and specific land uses as well as distinct reliefs [9,10]. The maximum, minimum, average, and stand deviation measured values at each compartment are shown in the Supplementary Material.

HC1 and HC2 stand out as areas of low dense community occupation within small planting areas, animal husbandry, drinking water, and multiple uses of water. HC3, in turn, represents an area where there are sand extraction sample numbers, a small aerodrome, and coconut processing industries for oil production.

The geospatial distribution of the physicochemical parameters showed that some of the SNs are not in line with the Brazilian Ministry of Health Consolidation Ordinance.

### 3.1.1. pH

The pH is a measure of the hydrogen concentration of the water, which is controlled by chemical reactions and the balance of ions present. According to the pH data, acidic waters are observed in the research sites. The raw water was measured in each Sample Number (SN) by the pH to show how acidic, neutral, or basic is the site. pHs of less than 7 indicate an acidic nature, whereas a pH of greater than 7 indicates an alkaline water, and it is a very important measurement concerning water quality. As the range moves from 0 to 14, with 7 being neutral, the PRC No. 5 recommends that the pH of water for the human supply stays around 6.0 to 9.0 [33]. reactions and the balance of ions present. According to the pH data, acidic waters are observed in the research sites. The raw water was measured in each Sample Number (SN) by the pH to show how acidic, neutral, or basic is the site. pHs of less than 7 indicate an acidic nature, whereas a pH of greater than 7 indicates an alkaline water, and it is a very important measurement concerning water quality. As the range moves from 0 to 14, with 7 being neutral, the PRC No. 5 recommends that the pH of water for the human supply stays around 6.0 to 9.0 [33].

*Water* **2020**, *12*, 623 6 of 18

The pH is a measure of the hydrogen concentration of the water, which is controlled by chemical

As shown on sites, the raw water is changing chemically, and the measured sites (SNs) that showed relative acidic water (4.2 to 6.9) have a higher amount of free hydrogen. Yet, the sites where the raw water had more free hydroxyl ions were considered as basic samples (7.1 to 8.4) (Figure 2). As shown on sites, the raw water is changing chemically, and the measured sites (SNs) that showed relative acidic water (4.2 to 6.9) have a higher amount of free hydrogen. Yet, the sites where the raw water had more free hydroxyl ions were considered as basic samples (7.1 to 8.4) (Figure 2).

**Figure 2.** Spatial distribution of pH on the dry and wet seasons in the three hydrological **Figure 2.** Spatial distribution of pH on the dry and wet seasons in the three hydrological compartments.

compartments. The pH represents a unit related to the activities of H + ions, which indicates in its expression's indices of neutrality, acidity or alkalinity. When comparing this parameter in the three hydrological compartments in the two climatic periods (wet and dry), there is a predominance of acidic pH in the hydrological compartment 1 (HC1), while the hydrological compartment 2 (HC2) presented an index of alkalinity in the studied period, which is similar to work conducted on groundwater as an The pH represents a unit related to the activities of H + ions, which indicates in its expression's indices of neutrality, acidity or alkalinity. When comparing this parameter in the three hydrological compartments in the two climatic periods (wet and dry), there is a predominance of acidic pH in the hydrological compartment 1 (HC1), while the hydrological compartment 2 (HC2) presented an index of alkalinity in the studied period, which is similar to work conducted on groundwater as an alternative source to an irregular surface water in Namaqualand, South Africa [3].

alternative source to an irregular surface water in Namaqualand, South Africa [3]. The pH of groundwater is influenced by salts, acids, and bases present in the environment. In the studied area, the pH may result from natural geologic-soil-water interactions, or trace the environment's quality, which includes the water source, land degradation, or deforestation [2,3,6], among other factors that occur in the region. In the environment, aquatic systems showing low pH values may be related to weathering processes [3]. Some lithological structures, when weathered, The pH of groundwater is influenced by salts, acids, and bases present in the environment. In the studied area, the pH may result from natural geologic-soil-water interactions, or trace the environment's quality, which includes the water source, land degradation, or deforestation [2,3,6], among other factors that occur in the region. In the environment, aquatic systems showing low pH values may be related to weathering processes [3]. Some lithological structures, when weathered, contribute by releasing acid-forming elements [6].

contribute by releasing acid-forming elements [6]. The minimum values of pHs < 6 listed in the Supplementary Material were shown to be more prominent on sites SNs at HC1 during the dry season (Figure 2). In addition, one of the reasons for pH values stays less than 6 and is the higher concentration of clay minerals, which dissolve releasing silica and aluminum in the waters [6]. This parameter directly influences the distribution of elements The minimum values of pHs < 6 listed in the Supplementary Material were shown to be more prominent on sites SNs at HC1 during the dry season (Figure 2). In addition, one of the reasons for pH values stays less than 6 and is the higher concentration of clay minerals, which dissolve releasing silica and aluminum in the waters [6]. This parameter directly influences the distribution of elements and chemical compounds in their free and ionized forms, which gives water an ability to increase or reduce its potential solubility relative to substances, including those with a degree of toxicity [1]. One of the factors that also contribute to acidic pH indices in water is the concentration of organic acids from dissolution resulting from the decomposition of organic matter, which may be happening in the hydrological compartment 1 (HC1), where there is a stream with visible contribution of organic matter. Additionally, the pH of rain varies between 5.0 and 6.0, but also "acid rain" may reach pH values as low as 4.3 [1,33]. The highest pH values (e.g., <7.5) are most likely related to weathering of the carbonate rocks that are represented in the study area [29]. from dissolution resulting from the decomposition of organic matter, which may be happening in the hydrological compartment 1 (HC1), where there is a stream with visible contribution of organic matter. Additionally, the pH of rain varies between 5.0 and 6.0, but also "acid rain" may reach pH values as low as 4.3 [1,33]. The highest pH values (e.g., <7.5) are most likely related to weathering of the carbonate rocks that are represented in the study area [29].

*Water* **2020**, *12*, 623 7 of 18

pH values stays less than 6 and is the higher concentration of clay minerals, which dissolve releasing silica and aluminum in the waters [6]. This parameter directly influences the distribution of elements

of the factors that also contribute to acidic pH indices in water is the concentration of organic acids

pH is a very sensitive component to changes and variations in water resources, and may oscillate, according to the dissolution of salts, decomposed organic matter, leaching processes, lithological soil types [29,37], and, above all, due to temperature [6,38]. pH is a very sensitive component to changes and variations in water resources, and may oscillate, according to the dissolution of salts, decomposed organic matter, leaching processes, lithological soil types [29,37], and, above all, due to temperature [6,38].

#### 3.1.2. Turbidity 3.1.2. Turbidity

The appearance of water with a turbidity less than 5 NTU is acceptable by the Brazilian Standards. The turbidity can be caused by particulate matter that may be present in the water source by resuspension of sediment along the flow path, or by the presence of inorganic particulate matter in groundwater [39,40]. Observing the turbidity in water (Figure 3) in both climatic periods (values in the Supplementary Material), there is an average value as a high turbidity during the dry season in SN 2, 7, 8 (HC1), 4, 9 (HC2), and 11, 13 (HC3) and on wet season in SN 10, 11, and 13 (HC3). The higher values of the parameter during the dry season may be related to the permanence of the suspended material that confers water turbidity through the lower rainfall of the dry season. Usually, clastic suspended sediments, such as sand and silt, give high turbidity in the waters [39–41], which are characteristic of the region where the study area is inserted. The appearance of water with a turbidity less than 5 NTU is acceptable by the Brazilian Standards. The turbidity can be caused by particulate matter that may be present in the water source by resuspension of sediment along the flow path, or by the presence of inorganic particulate matter in groundwater [39,40]. Observing the turbidity in water (Figure 3) in both climatic periods (values in the Supplementary Material), there is an average value as a high turbidity during the dry season in SN 2, 7, 8 (HC1), 4, 9 (HC2), and 11, 13 (HC3) and on wet season in SN 10, 11, and 13 (HC3). The higher values of the parameter during the dry season may be related to the permanence of the suspended material that confers water turbidity through the lower rainfall of the dry season. Usually, clastic suspended sediments, such as sand and silt, give high turbidity in the waters [39–41], which are characteristic of the region where the study area is inserted.

**Figure 3.** Spatial distribution of turbidity in the dry and wet seasons in the three hydrological compartments. **Figure 3.** Spatial distribution of turbidity in the dry and wet seasons in the three hydrological compartments.

Considering the maximum values observed, the turbidity was higher than 5 NTU in almost all sites evaluated in both seasons except on SN 15 and SN 16 (HC2) during the dry season and SN9 (HC2) during the wet season (Supplementary Material). Anthropic actions have a direct relationship under high turbidity rates (SN11 and SN12) (Figure 3). Leaching and hauling of particles from soil exposed by mining activities, such as sand extraction, provoke upturning, which increases the Considering the maximum values observed, the turbidity was higher than 5 NTU in almost all sites evaluated in both seasons except on SN 15 and SN 16 (HC2) during the dry season and SN9 (HC2) during the wet season (Supplementary Material). Anthropic actions have a direct relationship under high turbidity rates (SN11 and SN12) (Figure 3). Leaching and hauling of particles from soil exposed by mining activities, such as sand extraction, provoke upturning, which increases the availability of suspended particles. Environments where turbidity values are high are difficult for light to penetrate in water, which impairs the action of photosynthetic organisms. Some microorganisms are physically protected by the turbidity particles, which reduces the efficiency of water treatment [39–42]. Microbes

and other colloidal particles can be physically removed from water by various processes. The sizes of the microbes are especially important for their removal by sedimentation and filtration. Such methods are described in report APHA methodology [32]. by various processes. The sizes of the microbes are especially important for their removal by sedimentation and filtration. Such methods are described in report APHA methodology [32].

water treatment [39–42]. Microbes and other colloidal particles can be physically removed from water

*Water* **2020**, *12*, 623 8 of 18

availability of suspended particles. Environments where turbidity values are high are difficult for

### 3.1.3. Color 3.1.3. Color

The color in water can be caused by dissolved and/or suspended materials, and a brown shade often comes from rust in the water pipes. The physico-chemical characterization of color of sampling sites that exceeded the thresholds of PRC n◦5, with color > 100 UH, were SN 3, SN 11, and SN 13 during the wet season and SN 3 and SN11 with color 67 UH and 80.9 UH, respectively, during the dry season. Studies refer that the presence of organic matter, metals, and other chemical and biological components can cause changes in color values in surface water and in groundwater [43]. To understand color, it is important to deepen the characterization of sampling sites where the anomalous values occurred. The color in water can be caused by dissolved and/or suspended materials, and a brown shade often comes from rust in the water pipes. The physico-chemical characterization of color of sampling sites that exceeded the thresholds of PRC nº5, with color > 100 UH, were SN 3, SN 11, and SN 13 during the wet season and SN 3 and SN11 with color 67 UH and 80.9 UH, respectively, during the dry season. Studies refer that the presence of organic matter, metals, and other chemical and biological components can cause changes in color values in surface water and in groundwater [43]. To understand color, it is important to deepen the characterization of sampling sites where the anomalous values occurred.

The characteristics of site 3 are illustrated in Figure 4a. There is vegetation on the banks and a substantial amount of decomposing organic matter is suspended in the stream surface. The site number 11 (Figure 4b) is a lagoon where massive sand exploration occurs, and high turbidity occurred. Sample number 13 (Figure 4c) is a shallow cacimba without any fence and is situated in a contaminated area in a sanitation-free community surrounded by vegetation. Water is used for multiple uses and there is a makeshift toilet and laundry facility located within meters. The characteristics of site 3 are illustrated in Figure 4a. There is vegetation on the banks and a substantial amount of decomposing organic matter is suspended in the stream surface. The site number 11 (Figure 4b) is a lagoon where massive sand exploration occurs, and high turbidity occurred. Sample number 13 (Figure 4c) is a shallow cacimba without any fence and is situated in a contaminated area in a sanitation-free community surrounded by vegetation. Water is used for multiple uses and there is a makeshift toilet and laundry facility located within meters.

**Figure 4.** Sampling sites with anomalous color: (**a**) in the creek (SN3), (**b**) in the drainage pond of an area with sand exploration (SN 11), and (**c**) in an open cacimba (SN13). **Figure 4.** Sampling sites with anomalous color: (**a**) in the creek (SN3), (**b**) in the drainage pond of an area with sand exploration (SN 11), and (**c**) in an open cacimba (SN13).

According to Figure 5, the hydrological compartments HC1 and HC3 possibly presented the highest color indices, which is above what is allowed by the comparison norms (PRCn.5 threshold > B15UH). The results indicate that a probable source of color in these compartments is, as well as turbidity, responsible for the dissolution of organic substances that confer color and natural According to Figure 5, the hydrological compartments HC1 and HC3 possibly presented the highest color indices, which is above what is allowed by the comparison norms (PRCn.5 threshold > B15UH). The results indicate that a probable source of color in these compartments is, as well as turbidity, responsible for the dissolution of organic substances that confer color and natural pigmentation.

pigmentation. Color as well as turbidity is a parameter influenced by natural factors such as decaying organic matter and substances dissolved in water. By analyzing the geospatial distribution of this parameter in the three hydrological compartments, it was shown that a similar profile in the two climatic periods Color as well as turbidity is a parameter influenced by natural factors such as decaying organic matter and substances dissolved in water. By analyzing the geospatial distribution of this parameter in the three hydrological compartments, it was shown that a similar profile in the two climatic periods occurred and was analyzed (Figure 5).

occurred and was analyzed (Figure 5). At the same sampling sites that showed color values outside the Brazilian Standards, turbidity was also out of the standard range. High turbidity was measured at the SN11 (> 100 UT), where sand extraction occurred for years. The emerging issue of sand extraction and solutions to address potential environmental impact is of great concern. Some authors discuss environmental stressors related to the exploitation of sand in which one of them is the alteration of the surface and groundwater quality since this activity is capable of dissolving, suspending, and transporting organic and inorganic substances, which changes several quality parameters [41–43]. Turbidity is one of the At the same sampling sites that showed color values outside the Brazilian Standards, turbidity was also out of the standard range. High turbidity was measured at the SN11 (> 100 UT), where sand extraction occurred for years. The emerging issue of sand extraction and solutions to address potential environmental impact is of great concern. Some authors discuss environmental stressors related to the exploitation of sand in which one of them is the alteration of the surface and groundwater quality since this activity is capable of dissolving, suspending, and transporting organic and inorganic substances, which changes several quality parameters [41–43]. Turbidity is one of the most affected parameters, and very little attention has been given to this parameter related to sand extraction.

extraction.

**Figure 5.** Spatial distribution of color on the dry and wet seasons in the three hydrological **Figure 5.** Spatial distribution of color on the dry and wet seasons in the three hydrological compartments.

### compartments. 3.1.4. Total Iron

3.1.4. Total Iron Total iron was non-standard in four sampling sites. During the dry season, iron values were above the legal thresholds in sample numbers 1, 3, and 11 while, during the wet season, the affected sample numbers were SN3 and SN10. The SN1 presented relatively high values of iron where the soil and geology may be an influencing factor when considering the highest value of the well depth [6,43]. During the dry season, SN3 presented iron values of 4.5 mg/L and a pH of 5.41. It is common in waters that present this kind of pH values that the occurrence of Fe2+ iron indicates that water with a certain acidity presents high iron concentrations [38,44,45]. Iron in groundwater may contain ferrous iron at concentrations up to several milligrams per liter without discoloration or turbidity in the water when directly pumped from a well [45]. On exposure to the atmosphere, however, the ferrous iron oxidizes to ferric iron, which gives a reddish-brown color to the water. At levels above 0.3 mg/liter (SN1, SN3, and SN11), iron stains laundry and plumbing fixtures. There is usually no Total iron was non-standard in four sampling sites. During the dry season, iron values were above the legal thresholds in sample numbers 1, 3, and 11 while, during the wet season, the affected sample numbers were SN3 and SN10. The SN1 presented relatively high values of iron where the soil and geology may be an influencing factor when considering the highest value of the well depth [6,43]. During the dry season, SN3 presented iron values of 4.5 mg/L and a pH of 5.41. It is common in waters that present this kind of pH values that the occurrence of Fe2+ iron indicates that water with a certain acidity presents high iron concentrations [38,44,45]. Iron in groundwater may contain ferrous iron at concentrations up to several milligrams per liter without discoloration or turbidity in the water when directly pumped from a well [45]. On exposure to the atmosphere, however, the ferrous iron oxidizes to ferric iron, which gives a reddish-brown color to the water. At levels above 0.3 mg/liter (SN1, SN3, and SN11), iron stains laundry and plumbing fixtures. There is usually no noticeable taste at iron concentrations below 0.3 mg/l even though turbidity and color may develop [24].

noticeable taste at iron concentrations below 0.3 mg/l even though turbidity and color may develop [24]. The physico-chemical characterization of total iron of sampling sites that exceeded the The physico-chemical characterization of total iron of sampling sites that exceeded the thresholds of PRC n◦5, > 0.3 mg/L were SN3 (0.6 mg/L) and SN 10 (0.53 mg/L) during the wet season, and SN1 (0.44 mg/L), SN3 (4.5 mg/L), and SN 11 (0.61 mg/L).

thresholds of PRC nº5, > 0.3 mg/L were SN3 (0.6 mg/L) and SN 10 (0.53 mg/L) during the wet season, and SN1 (0.44 mg/L), SN3 (4.5 mg/L), and SN 11 (0.61 mg/L). Observing the values of the total iron in both wet and dry seasons (Figure 6), an increase in concentration during the dry season is noticeable, especially in the SN 3, in hydrological compartments of HC1. In this compartment, sample number 10 showed values that did not comply with the compared norms. Hydrological compartments HC2 and HC3 did not show a significant Observing the values of the total iron in both wet and dry seasons (Figure 6), an increase in concentration during the dry season is noticeable, especially in the SN 3, in hydrological compartments of HC1. In this compartment, sample number 10 showed values that did not comply with the compared norms. Hydrological compartments HC2 and HC3 did not show a significant change in iron concentrations in both climatic periods even though this element is not in accordance with the Brazilian Ministry of Health ordinance [33].

change in iron concentrations in both climatic periods even though this element is not in accordance with the Brazilian Ministry of Health ordinance [33]. Iron can also arise from corrosion of ferrous pipework and chemicals used in treatment processes (coagulation). Iron suspensions cause aesthetic problems including metallic taste and discoloration of water fittings and laundry. The Brazilian drinking water quality regulations include national Iron can also arise from corrosion of ferrous pipework and chemicals used in treatment processes (coagulation). Iron suspensions cause aesthetic problems including metallic taste and discoloration of water fittings and laundry. The Brazilian drinking water quality regulations include national standards for iron (0.3 mg/L), which can be removed from water by filtration, oxidation, coagulation, and sedimentation [32,33].

standards for iron (0.3 mg/L), which can be removed from water by filtration, oxidation, coagulation, and sedimentation [32,33]. Generally, in groundwater, iron levels derive from minerals and sediments that can be present in particulate or dissolved forms [1,38]. In surface waters, iron levels increase in the wet season as a result of soil runoff and erosion due to higher precipitation. Iron dissolved in water sets color, odor, and taste.

Iron concentrations in the dry season are higher than during the wet season, mainly in HC1. In both seasons, most of the area is higher than 0.3 mg/L (standard value for drinking water), but, during the dry season, much higher concentrations (up to 0.58 mg/L) can be encountered in groundwater from water wells, cacimbas, and surface water. result of soil runoff and erosion due to higher precipitation. Iron dissolved in water sets color, odor, and taste. Iron concentrations in the dry season are higher than during the wet season, mainly in HC1. In both seasons, most of the area is higher than 0.3 mg/L (standard value for drinking water), but, during the dry season, much higher concentrations (up to 0.58 mg/L) can be encountered in groundwater from water wells, cacimbas, and surface water.

*Water* **2020**, *12*, 623 10 of 18

in particulate or dissolved forms [1,38]. In surface waters, iron levels increase in the wet season as a

**Figure 6.** Spatial distribution of total iron during the dry and wet seasons in the three hydrological compartments. **Figure 6.** Spatial distribution of total iron during the dry and wet seasons in the three hydrological compartments.

The sample number SN3 during the wet season (0.6 mg/L) and the sample numbers SN3 (4.5 mg/L) and SN11 (0.61mg/L) during the dry season showed a common phenomenon where high iron levels that occurred showed high color concentrations [38]. It is usually possible to find high levels of iron in groundwater whose pH has acidity, and, in surface water, that has organic matter [44,45]. The sample number SN3 during the wet season (0.6 mg/L) and the sample numbers SN3 (4.5 mg/L) and SN11 (0.61 mg/L) during the dry season showed a common phenomenon where high iron levels that occurred showed high color concentrations [38]. It is usually possible to find high levels of iron in groundwater whose pH has acidity, and, in surface water, that has organic matter [44,45].

With this distribution, a higher concentration of acidity pH in the hydrological compartment 1 (HC1) is present (Figure 2). This same compartment presented levels of turbidity, color, and total iron out of the PRC threshold [33] (Supplementary Material and Figures 3, 5, and 6). With this distribution, a higher concentration of acidity pH in the hydrological compartment 1 (HC1) is present (Figure 2). This same compartment presented levels of turbidity, color, and total iron out of the PRC threshold [33] (Supplementary Material and Figures 3, 5 and 6).

### 3.1.5. Correlation Analysis 3.1.5. Correlation Analysis

The chart on 7, 8, and 9 contains the star plots of seven water parameters. The variable list for the sample star plot is Total Dissolved Solids (TDS), Temperature (Temp), pH, Oxi-redox potential (ORP), Electric Conductivity (Cond), Turbidity (Turb), and Dissolved Oxygen (OD). The plots were analyzed individually to identify clusters of the water parameters with similar features. The star plot of the water quality parameters compares the variables during dry and wet seasons at the three hydrological compartments (HC1, HC2, and HC3). The star plot in Figure 7 predicts the concentration of the parameters in HC1 and their comparison with drinking water standards in each sample number. In the hydrological compartment 1 (HC1), there is a significant correlation in STD and Cond in SN1 and SN2. The SN3 was higher during the dry season than during the wet season. The chart on 7, 8, and 9 contains the star plots of seven water parameters. The variable list for the sample star plot is Total Dissolved Solids (TDS), Temperature (Temp), pH, Oxi-redox potential (ORP), Electric Conductivity (Cond), Turbidity (Turb), and Dissolved Oxygen (OD). The plots were analyzed individually to identify clusters of the water parameters with similar features. The star plot of the water quality parameters compares the variables during dry and wet seasons at the three hydrological compartments (HC1, HC2, and HC3). The star plot in Figure 7 predicts the concentration of the parameters in HC1 and their comparison with drinking water standards in each sample number. In the hydrological compartment 1 (HC1), there is a significant correlation in STD and Cond in SN1 and SN2. The SN3 was higher during the dry season than during the wet season.

The dissolved oxygen (OD) was higher on SN1, SN3, SN6, and SN8 during the wet season. This is likely a consequence of gas exchange with the atmosphere, agitation, rainwater recharge, and increased movement of the water stream, which improves the aeration capacity of this ecosystem. Dissolved oxygen concentration is one of the most important factors for maintaining biodiversity in surface waters [46] such as rivers and streams. The dissolved oxygen (OD) was higher on SN1, SN3, SN6, and SN8 during the wet season. This is likely a consequence of gas exchange with the atmosphere, agitation, rainwater recharge, and increased movement of the water stream, which improves the aeration capacity of this ecosystem. Dissolved oxygen concentration is one of the most important factors for maintaining biodiversity in surface waters [46] such as rivers and streams.

A constant relationship in HC1, regardless of the season, is Total Dissolved Solids (TDS) with Electrical Conductivity (Cond). There is a similar response of both parameters, and, when this does

not occur, a significant reduction in turbidity is noted. This can be explained as described in water quality manuals where the dissolution of salts in water can result in electrolytes capable of conducting certain electrical current [23,46,47]. The temperature increase corresponds to a gradual increase in conductivity with a proportionality relationship [47]. All sample numbers (SN) show similar ecosystem performance, especially during the dry season. During the wet season, the correlation of the dissolved solids and conductivity became more evident. Electrical Conductivity (Cond). There is a similar response of both parameters, and, when this does not occur, a significant reduction in turbidity is noted. This can be explained as described in water quality manuals where the dissolution of salts in water can result in electrolytes capable of conducting certain electrical current [23,46,47]. The temperature increase corresponds to a gradual increase in conductivity with a proportionality relationship [47]. All sample numbers (SN) show similar ecosystem performance, especially during the dry season. During the wet season, the Electrical Conductivity (Cond). There is a similar response of both parameters, and, when this does not occur, a significant reduction in turbidity is noted. This can be explained as described in water quality manuals where the dissolution of salts in water can result in electrolytes capable of conducting certain electrical current [23,46,47]. The temperature increase corresponds to a gradual increase in conductivity with a proportionality relationship [47]. All sample numbers (SN) show similar ecosystem performance, especially during the dry season. During the wet season, the correlation of the dissolved solids and conductivity became more evident.

A constant relationship in HC1, regardless of the season, is Total Dissolved Solids (TDS) with

A constant relationship in HC1, regardless of the season, is Total Dissolved Solids (TDS) with

*Water* **2020**, *12*, 623 11 of 18

In hydrological compartment HC2 (Figure 8), in general, the correlated comportment between dissolved solids and conductivity was similar for both wet and dry seasons. The parameters showed the similar concentration mainly during the wet season. correlation of the dissolved solids and conductivity became more evident. In hydrological compartment HC2 (Figure 8), in general, the correlated comportment between dissolved solids and conductivity was similar for both wet and dry seasons. The parameters showed the similar concentration mainly during the wet season. In hydrological compartment HC2 (Figure 8), in general, the correlated comportment between dissolved solids and conductivity was similar for both wet and dry seasons. The parameters showed the similar concentration mainly during the wet season.

**Figure 7.** Correlation of water quality parameters during the wet and dry seasons in the hydrological compartment (HC) 1. Red line (dry season), blue line (wet season). **Figure 7.** Correlation of water quality parameters during the wet and dry seasons in the hydrological compartment (HC) 1. Red line (dry season), blue line (wet season). **Figure 7.** Correlation of water quality parameters during the wet and dry seasons in the hydrological compartment (HC) 1. Red line (dry season), blue line (wet season).

**Figure 8.** Correlation of water quality parameters during wet and dry seasons in a hydrological **Figure 8.** Correlation of water quality parameters during wet and dry seasons in a hydrological compartment (HC) 2. Red line (dry season), blue line (wet season). **Figure 8.** Correlation of water quality parameters during wet and dry seasons in a hydrological compartment (HC) 2. Red line (dry season), blue line (wet season).

compartment (HC) 2. Red line (dry season), blue line (wet season).

The pH at all sample numbers from HC2 remained in the same ranges. However, the Oxi-Redox water Potential (ORP) had a different value when compared to HC1. The potential for oxir-reduction with the tendency of substances to receive electrons can be evaluated, and it is possible to determine the potential by the possibility of the microorganisms grown in the area by the values of the redox potential, as shown on research on ORP in environmental research [48]. The ORP values in this compartment were similar to the values studied [48], within this premise, in a sample of natural water at a pH of 7. Oxygen should be the main electron receptor when the measured redox potential is close to (and above) + 400 mV. When the ORP value is between +100 and +300 mV, all oxygen must have been consumed and the main electron receptors will be NO<sup>3</sup> − and Mn, respectively, with the most abundant products being nitrogen and ammonia in addition to solubilizing manganese in the form of Mn2+. In more drastic anoxic conditions, ranging from 0 to –300 mV, the electron receptors will be Fe3+, then SO4, and, finally, organic matter and CO2, which are generated as reduction products iron (II), sulfide, and methane, respectively. It can be characterized as stable for both seasons [44,48]. Where there is a greater availability of dissolved oxygen (OD), there is an increase of ORP, which is a clear correlation in both climatic periods in this compartment. *Water* **2020**, *12*, 623 12 of 18 potential, as shown on research on ORP in environmental research [48]. The ORP values in this compartment were similar to the values studied [48], within this premise, in a sample of natural water at a pH of 7. Oxygen should be the main electron receptor when the measured redox potential is close to (and above) + 400 mV. When the ORP value is between +100 and +300 mV, all oxygen must have been consumed and the main electron receptors will be NO3− and Mn, respectively, with the most abundant products being nitrogen and ammonia in addition to solubilizing manganese in the form of Mn2+. In more drastic anoxic conditions, ranging from 0 to –300 mV, the electron receptors will be Fe3+, then SO4, and, finally, organic matter and CO2, which are generated as reduction products iron (II), sulfide, and methane, respectively. It can be characterized as stable for both seasons [44,48]. Where there is a greater availability of dissolved oxygen (OD), there is an increase of ORP, which is

In hydrological compartment 3 (HC3) (Figure 9), there is a similar result of the correlations that occurred on HC1 and HC2. All sample numbers have a low correlation to turbidity except on sample number 11, which is a lagoon from which sand has been extracted for years. At this point, it is possible to notice what differs from all the sample points, with low correlation on turbidity (Turb) and high correlation during the wet season for STD and pH. Similar results were found in a natural water reservoir [49,50]. a clear correlation in both climatic periods in this compartment. In hydrological compartment 3 (HC3) (Figure 9), there is a similar result of the correlations that occurred on HC1 and HC2. All sample numbers have a low correlation to turbidity except on sample number 11, which is a lagoon from which sand has been extracted for years. At this point, it is possible to notice what differs from all the sample points, with low correlation on turbidity (Turb) and high correlation during the wet season for STD and pH. Similar results were found in a natural water reservoir [49,50].

**Figure 9.** Correlation of water quality parameters during wet and dry seasons in the hydrological compartment (HC) 3. Red line (dry season) and blue line (wet season). **Figure 9.** Correlation of water quality parameters during wet and dry seasons in the hydrological compartment (HC) 3. Red line (dry season) and blue line (wet season).

In all three hydrological compartments, there was a correlation between STD, pH, Cond and OD, mainly during the wet season. Similar results were verified in natural and anthropic influences on groundwater quality of public wells in São Paulo State, Brazil [28,47,50]. It was observed in all studied sample numbers an increase of the concentrations of the parameters evaluated, mainly during the wet season. In the summer, there is a shorter dwell time of water in the water table [50,51]. This is explained by the lower percolation of rainwater. Thus, it is understandable to increase the values of the physicochemical parameters evaluated during the wet season. In all three hydrological compartments, there was a correlation between STD, pH, Cond and OD, mainly during the wet season. Similar results were verified in natural and anthropic influences on groundwater quality of public wells in São Paulo State, Brazil [28,47,50]. It was observed in all studied sample numbers an increase of the concentrations of the parameters evaluated, mainly during the wet season. In the summer, there is a shorter dwell time of water in the water table [50,51]. This is explained by the lower percolation of rainwater. Thus, it is understandable to increase the values of the physicochemical parameters evaluated during the wet season.

rural area. The view toward increasing inequalities in the access to water and the concentration of

Considering all the parameters evaluated, the inequalities in access to water is verified and

Considering all the parameters evaluated, the inequalities in access to water is verified and should aggregate contextual aspects of the ecosystem during dry and wet seasons and in the demographic perspectives, which reflect the intrinsic characteristics of the Atlantic Forest Biome dynamics and the relationships that are established daily in the water resource use and anthropized rural area. The view toward increasing inequalities in the access to water and the concentration of the deficit in certain rural communities corroborates the research in the Brazilian northeast area in underground flows of nutrients and trace metals [21,26,43].

The water supply system in the area must improve, mainly in treatment for drinking water, and from the aquifer [52]. Considering the inequalities in accessing better quality drinking water, the water security assessment is not present and significant progress is expected from the government to determine the vulnerability to water access [7,18,23]. Assessments involving the treatment of water or a better land use management [53,54] to ensure better water quality is necessary to include basic services such as the water supply for multiple uses.

The contributions of the paper results in (1) providing further evidence on similar investigations by water monitoring those rural communities and by an environmental education service, (2) the paper introduces the land use problems of the anthropized rural areas in the Atlantic Forest Biome, (3) the results show an improvement on the performance measurements of a set of water parameter results and should encourage the view that the subsequent water quality mapping that is necessary for a water security system.

For instance, the Pernambuco government must be alert of the boundaries between unconfined/confined aquifers that mark the interfaces for recharge and, thereby, contaminants that must be included in a study of hydro chemical investigation of the aquifers to determine the trace and minor element anomalies and possibly delineation of the interfaces and its vulnerability [52,53].

The anthropogenic activities from the rural area at the Atlantic Forest Biome in Pernambuco contribute significantly to groundwater contamination and to a water security assessment of groundwater quality. The government must investigate the influences from anthropogenic activities on deep groundwater (SN11) reflected in two possible scenarios: mixing with deep geological features and vertical leakage of shallow groundwater. Additionally, intensive groundwater contamination treatment must be done and controlled by water supply companies. The public drinking-water supply using groundwater should be tested prior to being used for general consumption.

### *3.2. Microbiological Analyses*

### 3.2.1. Coliforms

Presence/absence tests were performed for total and thermotolerant (fecal) coliforms in all samples studied in the established climatic seasonal periods. Ministry of Health Consolidation Ordinance No. 5 establishes water as unfit for human consumption in the presence of coliforms, which requires a maximum permissible (MPN)/100 mL of thermotolerant coliforms [33]. All studied sample numbers (SN1–SN16) presented fecal coliforms in their waters, which make them unsuitable for direct human consumption. The waters studied in this region should only be used for drinking water after a previous disinfection water treatment.

Although water does not naturally provide conditions for the proliferation of pathogenic organisms, they survive long enough in the environment to occur in the water transmission [55]. Therefore, the presence of pathogenic microorganisms in water comes from the contamination by the anthropized area [55,56] from animal or human feces, which results from infiltration and wastewater among others [10].

The studied region does not have basic sanitation structures in the area of influence of this study. The study area receives impacts with the exploitation of sand. The population uses septic tanks as a sewage system including many without adopting the safe distances for the installation of such systems. It is also emphasized that, in areas near the collection sampling sites, there is livestock breeding, and, in some of them, there is no proper sealing of the wells.

The contamination of water in an anthropized rural area from the Atlantic Forest Biome is not favorable since, in the studied area, risks of contamination can occur from the reservoirs or even in the distribution networks created by the residents since the urban supply system does not reach these areas [55–57]. Thus, it is common for distinct sources of contamination to occur [57].

Pathogenic agents have several properties that distinguish them from other drinking water contaminants. If infection is established, pathogens multiply in their host. Certain pathogenic bacteria are also able to multiply in food or beverages, and, thereby, perpetuate or even increase the chances of infection. The water studied must have a water quality verification that complements operational monitoring and assessments of contamination risks such as through auditing of treatment works and evaluating the process control and sanitary inspection. Water intended for human consumption should contain no indicator organisms. In most cases, monitoring for indicator bacteria provides a high degree of safety because of their large numbers in polluted waters [55–57].

### 3.2.2. Heterotrophic Bacterial Count

The Brazilian Ministry of Health Consolidation Ordinance No. 5 performed bacterial counting test analyses, which established a maximum permitted number of less than 500 CFU/mL [33]. From all sample numbers, six presented adequate values within the potability standards. They are the sample number 1, 6, 7, 8, 10, and 15. All other sample numbers presented inadequate values for human consumption. This is an important result because contaminants indicate a high probability of pathogenic microorganisms occurring in water [55]. Microorganisms responsible for gastroenteritis may occur, which are infections that have diarrhea as their main symptom.

### *3.3. Heavy Metals*

Heavy metals are harmful for freshwater ecosystems such as fish communities, and, therefore, can also be harmful to humans if they enter the food chain through the fish [53,54]. The Brazilian Standards for drinking water are established in the maximum permissible physicochemical parameter's concentrations [33]. Based on those values of the analyzed sample numbers (SN 1 to SN16), the heavy metal (Lead (0.01 mg/L −1 ), copper (2.00 mg/L), total chromium (0.05 mg/L), zinc (5.00 mg/L), cadmium (0.005 mg/L)), sulfates (250.00 mg/L), nitrate ( 10.00 mg/L), nitrite (1.00 mg/L), total hardness (500.0 mg/L), aluminum (0.20 mg/L), and ammonia (1.5 mg/L) had established an acceptable range for drinking water [21,33].

### **4. Conclusions**

The water quality in all studied sample numbers and hydrological compartments was not according to drinking water standards related to the parameters established by the PRC n◦5.

Drinking water sources in the rural area are contaminated, which can cause sickness and disease from pathogens. The presence of fecal coliforms is the most serious alteration in the quality of these waters since it configures contamination by pathogens that can lead to gastrointestinal diseases.

Drinking water sources are subject to contamination and require appropriate treatment to provide safe drinking water for their communities.

The color, pH, turbidity, and total iron parameters are the physicochemical parameters that classify, together with the microbiological ones, the waters studied as unfit for human consumption.

None of the collection sample numbers presented alterations for the analyzed heavy metals. There was a significant change in water quality between the dry and wet seasons, which comes with a similar quality pattern throughout the year. There was an increase in the values of some of the specific physicochemical parameters during the wet season.

In all three hydrological compartments, there was a concentration of different chemicals of interest not suitable for drinking water. Assessing water quality must involve monitoring the chemical concentrations in a government baseline concentration in a guideline established to protect human health or ecological communities.

The anthropized rural area from the Atlantic Forest Biome affects water quality. The Pernambuco water supply system must implement a regional simplified water treatment, so that people in the rural area can have access to better quality water. This would solve many of the water quality changes assessed in this paper.

For a water security assessment of groundwater quality, practices can be adopted for effective improvement of water quality evaluated in this research regarding the development of environmental public policies that monitor the quality of water used by the population in the region or the increase of the water and sanitation distribution network. Environmental education and the regulation of wells and cacimbas used by communities in the area are also key practices to ensure water quality within the parameters.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/3/623/s1, Table S1: Raw water maximum, minimum, average, and stand deviation measured values at each compartment (HC).

**Author Contributions:** Conceptualization: I.F.B.V., F.C.R.N., R.d.B.V.P., A.M.C.; Methodology: I.F.B.V., F.C.R.N., M.N.C., K.S.d.S.; Software: I.F.B.V., R.C.A.C.; Validation: I.F.B.V., F.C.R.N., R.C.A.C.; Formal analysis: I.F.B.V., F.C.R.N., A.M.C., M.N.C., R.d.B.V.P., T.C.T.P.; Investigation: I.F.B.V., F.C.R.N., A.M.C., T.C.T.P., F.A.L.P., L.F.S.F.; Writing—original draft preparation: I.F.B.V., F.C.R.N., A.M.C., T.C.T.P.; Writing—review and editing: F.C.R.N., F.A.L.P., L.F.S.F., A.M.C., T.C.T.P.; Supervision: F.C.R.N., A.M.C.; Project administration: F.C.R.N., A.M.C.; Funding acquisition: F.C.R.N., F.L.P., L.F.S.F., T.C.T.P. All authors have read and agree to the published version of the manuscript.

**Funding:** The present study was carried out within the framework of the Post-Graduation Research Programme of Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES); Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). For the author integrated with the CITAB Research Centre, it was further financed by the FEDER/COMPETE/POCI—Operational Competitiveness and Internationalization Programme, under Project POCI-01-0145-FEDER-006958, and by the National Funds of FCT—Portuguese Foundation for Science and Technology, under the project UIDB/AGR/04033/2020. For the author integrated in the CQVR, the National Funds of FCT—Portuguese Foundation for Science and Technology, under the project UIDB/QUI/00616/2020, supported the research. Land Use Policy Research Group – PolUS/UNESP/CNPq/Brazil.

**Acknowledgments:** We acknowledge administrative and technical support for sampling and water analyses and materials used for experiments from the Universities Laboratories, and Polo Regional Centro Norte, Pindorama in name of the Antonio Lúcio Mello Martins for helping in field sampling.

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

### **References**


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

### *Article* **Influence of Small Hydroelectric Power Stations on River Water Quality**

### **Xana Álvarez \* , Enrique Valero, Natalia de la Torre-Rodríguez and Carolina Acuña-Alonso**

Campus A Xunqueira s/n., Forestry Engineering School, University of Vigo, 36005 Pontevedra, Spain; evalero@uvigo.es (E.V.); nataliad@tcd.ie (N.d.l.T.-R.); carolalonso12@gmail.com (C.A.-A.)

**\*** Correspondence: xaalvarez@uvigo.es; Tel.: +34-986-801-959

Received: 29 November 2019; Accepted: 16 January 2020; Published: 21 January 2020

**Abstract:** Hydropower electricity generation is considered one of the cheapest technologies regarding electricity generation costs, and it is the most traditional, clean, renewable energy source. However, despite the environmental benefits offered by hydropower, they also can have negative impacts and consequences in the environment affecting water quality and disrupting river ecology. We investigated the environmental effects of four small hydropower plants (SPH) in north-west Spain by looking at the water quality of the four river stretches where the SPH plants are located. The physicochemical and biological characteristics of the water streams were analyzed, as well as the riparian ecological quality. Results showed that the presence of the hydropower plants did not significantly influence the physical and chemical characteristics of the water. There were no alterations of the benthic macroinvertebrate community at any of the plants except for one, and the riparian habitat was in general classified as good quality or close to natural conditions for all plants.

**Keywords:** hydroelectric power; physicochemical parameters; Iberian Bio-monitoring Working Party index; riparian forest quality

### **1. Introduction**

Between 2010 and 2040 energy demand is expected to increase by 56% worldwide [1], particularly, oil demand which is estimated to increase in countries, such as Brazil, China and India [2]. This is mainly due to the world's population unprecedented growth rate. According to the United Nations, by 2050, the world population will reach nine billion [3]. Fossil fuel combustion is currently the main source of energy production.

Two-thirds of the world's greenhouse gas emissions are produced by burning coal, natural gas and oil [4]. Anthropogenic greenhouse gas emissions are linked to climate change and human health problems [5]. Each year, more than 30 billion tons of carbon dioxide (CO2) are released into the atmosphere [6], and around 65% of the world's excess mortality—The number of deaths caused by a specific condition or exposure to harmful circumstances is directly associated with fossil fuel-related emissions [5]. According to The Lancet Commission on Pollution and Health [7], 9 million premature deaths were associated with pollution in 2015, representing 16% of all deaths worldwide. Environmental problems, such as ozone depletion, forest destruction and acid precipitation are direct consequences of the use of fossil fuels [8]. As suggested by the 2015 Paris Agreement, the average global temperature increase should be limited to 1.5–2 ◦C above pre-industrial levels [9] to avoid dangerous climate change. This will require measures, such as the phase-out of fossil fuels and geoengineering methods, such as CO<sup>2</sup> removal [10]. These options, however, raise issues, such as the costs involved, the environmental consequences and the process effectiveness [11,12]. In addition, Raftery et al. (2017) [13], estimated that global temperature increase would be between 2.0 and 4.9 ◦C by 2100, with a 5% chance that it will be less than 2 ◦C, and 1% chance that it will be less than 1.5 ◦C.

In the European Union, where electricity demand is growing, 40% of total energy consumption comes from the electricity sector [14,15]. In Europe, there are 450 GW of fossil-fuelled power plants, 200 GW each of coal and gas, and 50 GW of oil. These plants supply 40% of Europe's electricity, but release 1.4 GT of CO<sup>2</sup> emissions each year, which account for 30% of Europe's total emissions [15]. In Spain, 70% of the total energy consumption comes from abroad, and half of the energy generated in the country is supplied by imported combustible fuels [16]. This high energy dependency can cause problems of energy supply and affect wholesale electricity prices as these are linked to international combustible prices.

Renewable energy has been of interest since the oil crisis in the early 1970s [8]. Renewable energy technologies produce energy by converting natural resources into different types of energy. These technologies use the energy inherent in sunlight and its direct and indirect impacts on the Earth, such as wind, falling water, heating effects, and plant growth, gravitational forces, such as tides, and the heat of the Earth's core as the resources from which they produce energy [17]. Renewable energies have a positive impact in issues, such as the depletion of non-renewable energy sources, environmental problems (e.g., acid rain, ozone depletion) and the increase of energy use in developing countries [18]. They can supply-two thirds of global energy demands, and in combination with higher energy efficiency, could deliver 94% of the greenhouse emissions reduction that is necessary to comply with the Paris Agreement [19]. Moreover, the ability to access affordable energy at any time in the future is one of the main drivers for sustainable energy [20]. Because of this, renewable energies play an important role by contributing to mitigate climate change and is regarded as a potential solution to the environmental problems caused by the use of fossil fuels. In Europe, renewable energy development has been part of energy policy since 1986 when the Council of Europe established as one of its main targets to promote renewable energy sources direct [21]. The TERES II report, which aims to provide an update on renewable energy in Europe, estimated that by 2020 renewable energy could contribute as much as 14% of Europe's primary energy, creating around 500,000 jobs [22].

Hydroelectric power, known as hydropower, which is obtained by harnessing the power released when water passes through a vertical distance, has its origins in the pre-industrial revolution [23]. Hydropower electricity generation is considered one of the cheapest technologies regarding electricity generation costs [20], and it is the most traditional, clean, renewable energy source [24]. Hydropower plants have the capacity to respond to different energy demand fluctuations much faster than other systems, and they are able to convert direct mechanic work into electricity, making this technology highly efficient [20]. Hydropower is the leading source of renewable energy accounting for up to 70% of the world's supply [25]. However, despite the environmental benefits offered by hydropower, they also can have negative impacts and consequences in the environment, such as deforestation, aquatic and terrestrial biodiversity loss, affecting water quality and disrupting river ecology [26,27]. Yüksel [28] show that under the Renewables Scenario, 7250 MW of gas-fired capacity is substituted for 19,250 MW of wind and 1107 MW of small hydro over 2000–25. By 2025, all renewables combined (including large hydro) amount to more than 54 GW or 35% of installed capacity.

Although the energy industry in Galicia (Northwest of Spain) contributes to 8% of the region's gross domestic product (GDP) and is one of the main sources of employment, 77% of the region's primary energy is imported [29]. Due to its geographical location, Galicia has the potential to harness renewable energy through wind power, hydropower and biomass. In this work, we focus on hydropower and its effects in the ecosystem. Despite the environmental benefits of the small hydropower plants (SHP), they could negatively affect the ecosystem altering the natural hydrologic regime and the water quality, as well as the fish passage [24]. The aim of this study is, therefore, to investigate the environmental effects of four small hydropower plants (SHP) in north-west Spain by looking at the water quality of the four river stretches where the SHP plants are located. The physicochemical and biological (Iberian Bio-monitoring Working Party index) characteristics of the water streams of the SHP are analyzed, as well as the riparian ecological quality (QBR index). Based on the study of these two indices, it is intended to analyze how SPH affects the whole ecosystem in which they are found.

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

them in Galicia (Figure 1).

which they are found.

### *2.1. Study Areas and Parameters Studied*

The SHP plants are located in the towns of San Xusto (UTM: 29T 539378.405; 4708387.927), Hermida (UTM: 29T 537226.719; 4716799.905), Touro (UTM: 29T 562005.977; 4742137.597), and Gomil (UTM: 29T 580284.87; 4786188.633), in the Lerez, Umia, Ulla and Mandeo rivers, respectively, all of them in Galicia (Figure 1). *2.1. Study Areas and Parameters Studied*  The SHP plants are located in the towns of San Xusto (UTM: 29T 539378.405; 4708387.927), Hermida (UTM: 29T 537226.719; 4716799.905), Touro (UTM: 29T 562005.977; 4742137.597), and Gomil (UTM: 29T 580284.87; 4786188.633), in the Lerez, Umia, Ulla and Mandeo rivers, respectively, all of

*Water* **2020**, *12*, x FOR PEER REVIEW 3 of 19

(**b**)

**Figure 1.** *Cont.*

*Water* **2020**, *12*, x FOR PEER REVIEW 4 of 19

**Figure 1.** Location of the sample points. (**a**) Lerez river (**b**) Umia river (**c**) Ulla river (**d**) Mandeo river. Environmental characteristics assessed, parameters and index used, and location of the sampling points (grey circles). Riparian habitat sections, where the QBR (riparian forest quality index) was calculated (red). **Figure 1.** Location of the sample points. (**a**) Lerez river (**b**) Umia river (**c**) Ulla river (**d**) Mandeo river. Environmental characteristics assessed, parameters and index used, and location of the sampling points (grey circles). Riparian habitat sections, where the QBR (riparian forest quality index) was calculated (red).

The SHP plants consist of a diversion dam that blocks the river and diverts the water, a pipeline that draws the water from a higher level to the powerhouse, a penstock and a powerhouse building. The total power production of the San Xusto, Hermida, Touro and Gomil plants is 11.81 MW, 2.916 MW, 12.81 MW and 9 MW, respectively. The total production of the San Xusto, Hermida, Touro and Gomil plants is 37.81 GWh/year, 9.14 GWh/year, 30.89 GWh/year and 22.82 GWh/year, respectively. The ecological status of the stream water around San Xusto was previously studied and analyzed in Valero [24], and it is included in this work for comparison with the other three SHP plants. The SHP plants consist of a diversion dam that blocks the river and diverts the water, a pipeline that draws the water from a higher level to the powerhouse, a penstock and a powerhouse building. The total power production of the San Xusto, Hermida, Touro and Gomil plants is 11.81 MW, 2.916 MW, 12.81 MW and 9 MW, respectively. The total production of the San Xusto, Hermida, Touro and Gomil plants is 37.81 GWh/year, 9.14 GWh/year, 30.89 GWh/year and 22.82 GWh/year, respectively. The ecological status of the stream water around San Xusto was previously studied and analyzed in Valero [24], and it is included in this work for comparison with the other three SHP plants.

The Directive 2000/60/EC by The European Parliament [30], which establishes the framework in the field of water policy within the EU, was used in this work to assess the water status affected by the SHP plants. Thus, some of the water quality elements described in the directive were considered for this study. The physicochemical parameters taken into consideration were water temperature (T, °C), conductivity (mS/cm), dissolved oxygen (DO, mg/L) and pH (Figure 2). The Directive 2000/60/EC by The European Parliament [30], which establishes the framework in the field of water policy within the EU, was used in this work to assess the water status affected by the SHP plants. Thus, some of the water quality elements described in the directive were considered for this study. The physicochemical parameters taken into consideration were water temperature (T, ◦C), conductivity (mS/cm), dissolved oxygen (DO, mg/L) and pH (Figure 2).

The Iberian Biological Monitoring Working Party index (BMWP) score system was introduced in 1980 to provide an index of river water quality for England and Wales based on aquatic macroinvertebrates, was used to assess the biological quality of the water [31,32] (Figure 2). To do this, the taxa of macho-benthonic fauna found in the sampled areas were identified to the family level, a predefined score was allocated for each family, and the total BMWP score for a sample was the summation of the scores of all the families found. The scores of the BMWP (0–>100) were grouped in 5 quality classes [31]. In Ulla river, the Iberian Bio-monitoring Working Party index (IBMWP index) was applied in 1 and 3 sample sites because 2 and 4 were very deep waters, due to the embalming effect. The conditions are not suitable for collecting macroinvertebrates. IBMWP scores were then compared with conditions in Cantabric-Atlantic quality elements at high ecological status [33]. This index obtained optimal results in Thailand [34] United Kingdom [35], among others. The Iberian Biological Monitoring Working Party index (BMWP) score system was introduced in 1980 to provide an index of river water quality for England and Wales based on aquatic macroinvertebrates, was used to assess the biological quality of the water [31,32] (Figure 2). To do this, the taxa of macho-benthonic fauna found in the sampled areas were identified to the family level, a predefined score was allocated for each family, and the total BMWP score for a sample was the summation of the scores of all the families found. The scores of the BMWP (0–>100) were grouped in 5 quality classes [31]. In Ulla river, the Iberian Bio-monitoring Working Party index (IBMWP index) was applied in 1 and 3 sample sites because 2 and 4 were very deep waters, due to the embalming effect. The conditions are not suitable for collecting macroinvertebrates. IBMWP scores were then compared with conditions in Cantabric-Atlantic quality elements at high ecological status [33]. This index obtained optimal results in Thailand [34] United Kingdom [35], among others.

The Riparian Forest Quality index (QBR) was used to assess the ecological quality of the riparian habitats [36] (Figure 2). This index varies between 0 and 100, and it is based on four components of riparian habitat: Total riparian vegetation cover, cover structure, cover quality and channel alterations for the different geomorphology of the river from its headwaters to the lower reaches. The Riparian Forest Quality index (QBR) was used to assess the ecological quality of the riparian habitats [36] (Figure 2). This index varies between 0 and 100, and it is based on four components of riparian habitat: Total riparian vegetation cover, cover structure, cover quality and channel alterations for the different geomorphology of the river from its headwaters to the lower reaches. After completing

After completing the analysis, the sum of scores for the four components gives the final QBR index. The QBR index, despite being developed for Mediterranean forests can be adapted to other types of

riverside forests, as they have done in other studies [37–39].

the analysis, the sum of scores for the four components gives the final QBR index. The QBR index, despite being developed for Mediterranean forests can be adapted to other types of riverside forests, as they have done in other studies [37–39]. *Water* **2020**, *12*, x FOR PEER REVIEW 5 of 19

**Figure 2.** Schematic overview of the methodology. **Figure 2.** Schematic overview of the methodology.

### *2.2. Sampling*

*2.2. Sampling*  Four sample points each were selected at San Xusto and Hermida, and three sample points each were selected at Touro and Gomil. The points were selected in areas properly accessible and free of operational risks. A GPS GPSMAP 60CSx (Garmin, Olathe, KS, USA) was used to register the position of the sampling points. Four sample points each were selected at San Xusto and Hermida, and three sample points each were selected at Touro and Gomil. The points were selected in areas properly accessible and free of operational risks. A GPS GPSMAP 60CSx (Garmin, Olathe, KS, USA) was used to register the position of the sampling points.

The YSI 556 Handheld Multiparameter Instrument (YSI, Yellow Springs, OH, USA), was used to measure the water's physicochemical parameters in each sampling point. Different sensor types were used: The YSI Temperature Precision™ thermistor, the 4-electrode cell with auto-ranging, the steady-state polarographic, the platinum button, and the glass combination electrode. It was first submerged in the watercourse for twenty minutes for stabilization, and later left submerged (5 m deep) for another twenty minutes, while registering in order to obtain representative data. Measurements were made at the same time of the day. The YSI 556 Handheld Multiparameter Instrument (YSI, Yellow Springs, OH, USA), was used to measure the water's physicochemical parameters in each sampling point. Different sensor types were used: The YSI Temperature Precision™ thermistor, the 4-electrode cell with auto-ranging, the steady-state polarographic, the platinum button, and the glass combination electrode. It was first submerged in the watercourse for twenty minutes for stabilization, and later left submerged (5 m deep) for another twenty minutes, while registering in order to obtain representative data. Measurements were made at the same time of the day.

To analyze the biological quality of the water using the IBMWP index, a 50 m radius from the points were sampled. Point 1 at San Xusto, and point 2 at Touro were not sampled, because it was not easy to obtain enough representative data. Pictures were taken of all the sampled areas and the locations were registered with a GPS. The different types of habitats were identified according to depth, flow velocity, substratum, and vegetation. Then, each habitat was sampled for benthic invertebrates, and captures were identified over a white tray on the field using an atlas and identification key. After the identification, specimens were released back into the river. Benthos which could not be identified in the field were kept in bottles with alcohol at 70% to be identified at the laboratory. To analyze the biological quality of the water using the IBMWP index, a 50 m radius from the points were sampled. Point 1 at San Xusto, and point 2 at Touro were not sampled, because it was not easy to obtain enough representative data. Pictures were taken of all the sampled areas and the locations were registered with a GPS. The different types of habitats were identified according to depth, flow velocity, substratum, and vegetation. Then, each habitat was sampled for benthic invertebrates, and captures were identified over a white tray on the field using an atlas and identification key. After the identification, specimens were released back into the river. Benthos which could not be identified in the field were kept in bottles with alcohol at 70% to be identified at the laboratory.

To calculate the QBR index, the Munné protocol [36] was followed, but with some modifications to be adapted to the tree species composition of our study area. A river stretch of 400 m around the dam, and another river stretch of 650 m around the SHP plant and tail race were selected. These reaches were divided into sections with similar characteristics, and the QBR was calculated for each one. Points between two sections with different QBR were located with GPS. We calculated one QBR index for each river bank. Based on the index, it is not necessary to apply it to all the banks, but it is recommended to choose strategic points according to their characteristics. In this study, the riparian vegetation was evaluated in the points where some direct affections were made (because of the To calculate the QBR index, the Munné protocol [36] was followed, but with some modifications to be adapted to the tree species composition of our study area. A river stretch of 400 m around the dam, and another river stretch of 650 m around the SHP plant and tail race were selected. These reaches were divided into sections with similar characteristics, and the QBR was calculated for each one. Points between two sections with different QBR were located with GPS. We calculated one QBR index for each river bank. Based on the index, it is not necessary to apply it to all the banks, but it is recommended to choose strategic points according to their characteristics. In this study, the riparian vegetation was

construction phase). Intermediate points were considered in these cases where topographic

evaluated in the points where some direct affections were made (because of the construction phase). Intermediate points were considered in these cases where topographic conditions allow access to the river. In this last case, the embalming of the waters, due to the installation of the plants, will make the vegetation to change. However, it will be necessary for a considered period of time to evaluate these changes. The sampling period was carried out between December 2007 and November 2010. The physicochemical parameters and the biological quality of the water were sampled every three months, in March, June, September and December. Directed sampling has been performed—simple samples were collected in triplicate, and at the points already described. The samples were collected in amber borosilicate glass bottles previously washed with sample water. The QBR index was calculated every six months, in March and in September.

### *2.3. Analyses and Statistics*

The software Ecowatch was used to download and process all the physicochemical data and mean, standard deviation, maximum and minimum were calculated. At the laboratory, the macroinvertebrates collected were identified using a glass Motic ST-37 20–80 × (Motic, Xiamen, China). Mean values of each physicochemical variable of the series registered in each sampling point, as well as the IBMWP index, were plotted over time to assess if their values stabilized themselves across time, as well as to observe differences between points. Differences between years and points, as well as differences between the four SHP plants, were analyzed using Kruskal-Wallis tests. Kruskal-Wallis is a nonparametric statistical test that assesses the differences among three or more independently sampled groups on a single, non-normally distributed continuous variable [40].

### **3. Results**

The Kruskal-Wallis tests show no significant differences for the physicochemical parameters of pH, temperature, dissolved oxygen and conductivity between sampling points and years for any of the four plants, except for the temperature at the Touro plant (Table 1). No significant differences were obtained for the four water quality parameters, since the data do not present irregularities over time.


**Table 1.** Kruskal-Wallis results for each of the physicochemical parameters for each plant.

This is in agreement with other studies, such as those from [41,42]. Jesus et al. evaluated the impact of an SHP plant over a period of two years on the water quality of the Ardena river in Portugal, and found the presence of the hydropower plant did not significantly influence the physical and chemical characteristic of the water. Santos et al. also studied the effects of SHP plants in central and north Portugal and found no significant difference between the physicochemical parameters at the different study sites.

### *3.1. Temperature*

Mean water temperature for the study period was 12.41 ◦C, 12.26 ◦C, 13.24 ◦C and 11.06 ◦C at San Xusto, Hermida, Touro and Gomil, respectively. Temperature values were in concordant with those from the Reference Condition in Cantabric-Atlantic siliceous rivers from IPH 2008 [43] (Table 2), where it is specified that the river's elements of this study belongs to the siliceous Cantabrian-Atlantic rivers.


Dissolved oxygen (mg/L) 9.00 6.7

**Table 2.** Reference Condition in Cantabric-Atlantic siliceous rivers.

These results also comply with the requirements of Directive 2006/44/EC [44] on the quality of freshwaters needing protection or improvement in order to support fish life, by which thermal discharges must not cause the temperature downstream of the point of thermal discharge to exceed 21.5 ◦C in salmonid waters. Temperature limits were not exceeded in any sampling point for any period in any of the plants (Figure 3). Directive 2006/44/EC also states that thermal discharges must not cause the temperature downstream of the point of thermal discharge to exceed 10 ◦C during breeding periods of salmonid species, which need cold water for reproduction and only to waters that contain such species. In Galician waters, the period when salmon breed extends from November to January [45]. Temperatures registered for all sampling points between November and January was below 10 ◦C at San Xusto, Hermida, and Gomil. The temperature at the Touro plant was above 10 ◦C in all points for two of the sampling periods. However, the Kruskal-Wallis test showed no significant difference between the temperatures recorded for all the points during the study period. Conductivity (mS/cm) 0.04 <0.03 These results also comply with the requirements of Directive 2006/44/EC [44] on the quality of freshwaters needing protection or improvement in order to support fish life, by which thermal discharges must not cause the temperature downstream of the point of thermal discharge to exceed 21.5 °C in salmonid waters. Temperature limits were not exceeded in any sampling point for any period in any of the plants (Figure 3). Directive 2006/44/EC also states that thermal discharges must not cause the temperature downstream of the point of thermal discharge to exceed 10 °C during breeding periods of salmonid species, which need cold water for reproduction and only to waters that contain such species. In Galician waters, the period when salmon breed extends from November to January [45]. Temperatures registered for all sampling points between November and January was below 10 °C at San Xusto, Hermida, and Gomil. The temperature at the Touro plant was above 10 °C in all points for two of the sampling periods. However, the Kruskal-Wallis test showed no significant difference between the temperatures recorded for all the points during the study period.

**Figure 3.** *Cont.*

*Water* **2020**, *12*, x FOR PEER REVIEW 8 of 19

**Figure 3.** Temperature values measured during the survey period at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant). **Figure 3.** Temperature values measured during the survey period at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant).

#### *3.2. pH 3.2. pH*

According to the Directive 2006/44/EC, artificial pH variations with respect to the unaffected values shall not exceed ±0.5 of a pH unit within the limits falling between 6.0 and 9.0 provided that these variations do not increase the harmfulness of other substances present in the water. At San Xusto and Hermida, low pH values were recorded in all four points close to the lower limit allowed for salmonid waters (Figure 4). However, mean pH values for all points for the sampling period was 6.02 and 6.09 at San Xusto and Hermida, respectively, which is within the limits established by the directive. In addition, there were no significant differences between all observations recorded in any of the two plants. Mean pH values for all points at Touro was 6.55, falling within the limits established by the Directive 2006/44/EC. There were only two observations at Touro, in points 1 and 3 recorded during the same period, where the pH value was below 6. Kruskal-Wallis test showed a significant difference between the sampled points. At Gomil, mean pH values for all points was 6.40, and there was no significant difference between the sampled points. However, some observations showed values below 6, but these were close to the limit allowed. According to Reference Condition (Table 2), the quality of the water for the four river starches can be classified as good/moderate. The low pH values obtained for some of the points could be explained by the natural According to the Directive 2006/44/EC, artificial pH variations with respect to the unaffected values shall not exceed ±0.5 of a pH unit within the limits falling between 6.0 and 9.0 provided that these variations do not increase the harmfulness of other substances present in the water. At San Xusto and Hermida, low pH values were recorded in all four points close to the lower limit allowed for salmonid waters (Figure 4). However, mean pH values for all points for the sampling period was 6.02 and 6.09 at San Xusto and Hermida, respectively, which is within the limits established by the directive. In addition, there were no significant differences between all observations recorded in any of the two plants. Mean pH values for all points at Touro was 6.55, falling within the limits established by the Directive 2006/44/EC. There were only two observations at Touro, in points 1 and 3 recorded during the same period, where the pH value was below 6. Kruskal-Wallis test showed a significant difference between the sampled points. At Gomil, mean pH values for all points was 6.40, and there was no significant difference between the sampled points. However, some observations showed values below 6, but these were close to the limit allowed. According to Reference Condition (Table 2), the quality of the water for the four river starches can be classified as good/moderate.

characteristics of the waters in Galicia, which are generally acidic [46]. The low pH values obtained for some of the points could be explained by the natural characteristics of the waters in Galicia, which are generally acidic [46].

*Water* **2020**, *12*, x FOR PEER REVIEW 9 of 19

**Figure 4.** *Cont.*

*Water* **2020**, *12*, x FOR PEER REVIEW 10 of 19

**Figure 4.** pH values measured during the survey period at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant). **Figure 4.** pH values measured during the survey period at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant). **Figure 4.** pH values measured during the survey period at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant).

#### *3.3. Dissolved Oxygen 3.3. Dissolved Oxygen 3.3. Dissolved Oxygen*

Mean dissolved oxygen values for the study period was 11.19 mg/L, 11.17 mg/L, 9.60 mg/L and 11.03 mg/L at San Xusto, Hermida, Touro and Gomil, respectively. These values were higher than the Reference Condition (Table 2). The results are in agreement with Directive 2006/44/EC, which states that the dissolved oxygen of 50% of the samples recorded must be equal or above to 9 mg/L. San Xusto plant had 96% of values above 9 mg/L, Hermida plant had 94% of values above 9 mg/L, Touro plant had 55% of values above 9 mg/L, and Gomil plant had 97.82% of values 9 mg/L (Figure 5). Results for the Kruskal-Wallis test show no significant difference between points during the study period for any of the plants. Dissolved oxygen concentration was inversely linked to temperature, which matches results from [41]. Mean dissolved oxygen values for the study period was 11.19 mg/L, 11.17 mg/L, 9.60 mg/L and 11.03 mg/L at San Xusto, Hermida, Touro and Gomil, respectively. These values were higher than the Reference Condition (Table 2). The results are in agreement with Directive 2006/44/EC, which states that the dissolved oxygen of 50% of the samples recorded must be equal or above to 9 mg/L. San Xusto plant had 96% of values above 9 mg/L, Hermida plant had 94% of values above 9 mg/L, Touro plant had 55% of values above 9 mg/L, and Gomil plant had 97.82% of values 9 mg/L (Figure 5). Results for the Kruskal-Wallis test show no significant difference between points during the study period for any of the plants. Dissolved oxygen concentration was inversely linked to temperature, which matches results from [41]. Mean dissolved oxygen values for the study period was 11.19 mg/L, 11.17 mg/L, 9.60 mg/L and 11.03 mg/L at San Xusto, Hermida, Touro and Gomil, respectively. These values were higher than the Reference Condition (Table 2). The results are in agreement with Directive 2006/44/EC, which states that the dissolved oxygen of 50% of the samples recorded must be equal or above to 9 mg/L. San Xusto plant had 96% of values above 9 mg/L, Hermida plant had 94% of values above 9 mg/L, Touro plant had 55% of values above 9 mg/L, and Gomil plant had 97.82% of values 9 mg/L (Figure 5). Results for the Kruskal-Wallis test show no significant difference between points during the study period for any of the plants. Dissolved oxygen concentration was inversely linked to temperature, which matches results from [41].

(**a**) **Figure 5.** *Cont.*

**Figure 5.** Dissolved oxygen values measured during the survey period at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant). **Figure 5.** Dissolved oxygen values measured during the survey period at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant).

(**d**)

Point 1 Point 2 Point 3

#### *3.4. Conductivity Water* **2020**, *12*, x FOR PEER REVIEW 12 of 19

Mean conductivity values were 0.04 mS/cm at San Xusto and Hermida, 0.09 mS/cm at Touro and 0.07 mS/cm at Gomil. There were no significant differences between the points for any of the plants from the Kruskal-Wallis test. Our results are in compliance with the requirement from Reference Condition (Table 2) for San Xusto and Hermida, and are very close for Touro and Gomil. These results are also in agreement with previous studies from Costas et al. (2009) [47] and Ansemil and Membiela (1992) [46]. *3.4. Conductivity*  Mean conductivity values were 0.04 mS/cm at San Xusto and Hermida, 0.09 mS/cm at Touro and 0.07 mS/cm at Gomil. There were no significant differences between the points for any of the plants from the Kruskal-Wallis test. Our results are in compliance with the requirement from Reference Condition (Table 2) for San Xusto and Hermida, and are very close for Touro and Gomil. These results are also in agreement with previous studies from Costas et al. (2009) [47] and Ansemil and Membiela

### *3.5. Biological Status of the Water* (1992) [46].

A total of 12 samples were collected for each river. At San Xusto, IBMPW index had a lower score at point 3, located upstream the power plant, than at points 2 and 4, although the Kruskal-Wallis test showed no significant differences between points (Kruskal–Wallis test = 1.756; df = 2; *p* = 0.416), and 86.84% of the scores fell within Quality I category (unpolluted or not consider altered water). At Hermida, Kruskall-Wallis test showed no significant differences between points (Kruskal-Wallis chi-squared = 3.800, df = 3, *p*-value = 0.283), and 95.91% of the scores fell within Quality I category. At Gomil, Kruskal-Wallis test showed no significant differences between points (Kruskal-Wallis chi-squared = 0.286, df = 2, *p*-value = 0.866), and 86.95% of the scores fell in the Quality I category. These results (Figure 6) are in agreement with some authors, such as Copeman (1997) and Almodóvar and Nicola (1999), [48,49], who found there were no alterations of the benthic macroinvertebrate community, due to small hydropower plants. At Touro, on the other hand, scores were very low, only 8.69% of them falling within Quality I category. Mann-Whitney test showed a significant difference between points (*p* < 0.05). In Point 1, all scores fell within Quality V category (waters highly contaminated), whereas, in point 3 most of the scores fell within Quality III category (evidence of some contamination). It is important to point out that the low scores obtained in point 1 are probably not, due to the low biological quality of the water, but rather to the lack of representative areas when the sampling was done. *3.5. Biological Status of the Water*  A total of 12 samples were collected for each river. At San Xusto, IBMPW index had a lower score at point 3, located upstream the power plant, than at points 2 and 4, although the Kruskal-Wallis test showed no significant differences between points (Kruskal–Wallis test = 1.756; df = 2; *p* = 0.416), and 86.84% of the scores fell within Quality I category (unpolluted or not consider altered water). At Hermida, Kruskall-Wallis test showed no significant differences between points (Kruskal-Wallis chi-squared = 3.800, df = 3, *p*-value = 0.283), and 95.91% of the scores fell within Quality I category. At Gomil, Kruskal-Wallis test showed no significant differences between points (Kruskal-Wallis chi-squared = 0.286, df = 2, *p*-value = 0.866), and 86.95% of the scores fell in the Quality I category. These results (Figure 6) are in agreement with some authors, such as Copeman (1997) and Almodóvar and Nicola (1999), [48,49], who found there were no alterations of the benthic macroinvertebrate community, due to small hydropower plants. At Touro, on the other hand, scores were very low, only 8.69% of them falling within Quality I category. Mann-Whitney test showed a significant difference between points (*p* < 0.05). In Point 1, all scores fell within Quality V category (waters highly contaminated), whereas, in point 3 most of the scores fell within Quality III category (evidence of some contamination). It is important to point out that the low scores obtained in point 1 are probably not, due to the low biological quality of the water, but rather to the lack of representative areas when the sampling was done.

**Figure 6.** *Cont.*

**Figure 6.** Evolution of the water biological quality index IBMWP during the period of study: I (>120) very clean waters; I (101–120) unpolluted or not considerably altered water; II (61–100) moderate effects of pollution are evident; III (36–60) polluted water; IV (16–35) very polluted water; V (<15) heavily polluted water, at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant). **Figure 6.** Evolution of the water biological quality index IBMWP during the period of study: I (>120) very clean waters; I (101–120) unpolluted or not considerably altered water; II (61–100) moderate effects of pollution are evident; III (36–60) polluted water; IV (16–35) very polluted water; V (<15) heavily polluted water, at (**a**) Lerez (San Xusto plant), (**b**) Umia (Hermida plant), (**c**) Ulla (Touro plant) and (**d**) Mandeo (Gomil plant).

### *3.6. Ecological Status of the Riparian Zones 3.6. Ecological Status of the Riparian Zones*

A total of 7 samples were collected for each river. At San Xusto, the evaluation of the habitat according to the QBR index showed two river bank sections with extreme degradation and two with strong alteration, while in the remaining sections, the riparian habitat was classified as good quality or natural conditions (Figure 7). At Hermida, QBR index showed one section with extreme degradation, one section with strong alteration, while in the remaining sections, the riparian habitat was classified fair to good quality or even as having natural conditions (Figure 7). At Touro, the QBR index showed three sections with extreme degradation, while in the remaining sections, the riparian habitat was classified fair to good quality (Figure 7). At Gomil, the QBR index showed two sections with extreme degradation and four with strong alteration, while in the remaining sections, the riparian habitat was classified fair to good quality or even natural conditions (Figure 7). A total of 7 samples were collected for each river. At San Xusto, the evaluation of the habitat according to the QBR index showed two river bank sections with extreme degradation and two with strong alteration, while in the remaining sections, the riparian habitat was classified as good quality or natural conditions (Figure 7). At Hermida, QBR index showed one section with extreme degradation, one section with strong alteration, while in the remaining sections, the riparian habitat was classified fair to good quality or even as having natural conditions (Figure 7). At Touro, the QBR index showed three sections with extreme degradation, while in the remaining sections, the riparian habitat was classified fair to good quality (Figure 7). At Gomil, the QBR index showed two sections with extreme degradation and four with strong alteration, while in the remaining sections, the riparian habitat was classified fair to good quality or even natural conditions (Figure 7).

*Water* **2020**, *12*, x FOR PEER REVIEW 14 of 19

(**a**)

(**b**)

**Figure 7.** *Cont.*

*Water* **2020**, *12*, x FOR PEER REVIEW 15 of 19

(**c**)

### *3.7. General Discussion 3.7. General Discussion*

Renewable energies, such as the one provided by SPH are a clean source of energy. Regarding the four plants studied in the present work, it must be taken into account that their installation in the last decades was done without prior planning, and without any specific studies and analyses on the ecological quality of the river environment. Therefore, without these previous studies, it is very difficult to predict and know what impacts these SPHs have on the ecosystem. Renewable energies, such as the one provided by SPH are a clean source of energy. Regarding the four plants studied in the present work, it must be taken into account that their installation in the last decades was done without prior planning, and without any specific studies and analyses on the ecological quality of the river environment. Therefore, without these previous studies, it is very difficult to predict and know what impacts these SPHs have on the ecosystem.

For better planning, the idea would have been to carry out studies prior to the construction, and also during the different construction stages to have a better understanding of the final impacts. These For better planning, the idea would have been to carry out studies prior to the construction, and also during the different construction stages to have a better understanding of the final impacts.

studies should have been planned in four distinct phases corresponding to: Phase 0 (prior to construction), construction phase, exploitation phase, and abandonment or dismantling phase.

These studies should have been planned in four distinct phases corresponding to: Phase 0 (prior to construction), construction phase, exploitation phase, and abandonment or dismantling phase.

In this way, we would have had information on how the ecosystem would change from its original state (on a small scale and on a global scale). We would also be able to know what impacts the construction works generate. The construction can generate environmental problems, such as spills and reduction of flows, which can have specific effects on some species. Finally, we could have made an evaluation of the evolution of the ecosystem once the plant began to work, so the impact of the operation of the plant itself (water diversion, heating of the water, etc.) would have been known. In the case of the long-term impacts that SPHs can cause, analyses over several years would be recommended.

Small impacts have been detected with this study, but they are short term since the sampling data are from 2007 to 2010. This study focuses on the impacts that the SPHs operation may have on the ecological quality of water and riverbank vegetation in the short term. What we want to highlight is that there are impacts that occur in the aquatic ecosystem, therefore, it is important to know to what extent they occur. In this way, preventive measures would be taken from the knowledge that these analyses can generate. On the other hand, to know the extent of these impacts, it is very important to have a large-scale data collection. In this way, water managers, energy companies and users., can make changes in their trends, consumption and models.

### **4. Conclusions**

Results showed uniformity of the data obtained. Temperature limits were not exceeded at any point in any of the plants. Temperatures registered during the salmon breeding period was below 10 ◦C for all the plants with the exception of Touro. San Xusto and Hermida plants showed low pH values. However, mean pH at all the plants fell within the limits established. The low pH values obtained for some of the points could be explained by the natural characteristics of the waters in Galicia, which are generally acidic [50]. Mean dissolved oxygen values for the study period were in agreement with Directive 2006/44/EC, which states that 50% of the samples must be equal to or greater than 9 mg/L. Mean conductivity samples were in compliance with requirements outlined in Reference Condition for San Xusto and Hermida, and were very close for Touro and Gomil. The IBMPW results showed no alterations of the benthic macroinvertebrate community at any of the plants except for Touro, where scores were very low, although this could be explained, due to the lack of representative areas when the sampling was done. The riparian habitat (QBR index) was in general classified as good quality or close to natural conditions for all plants. The lack of data prior to the construction of the power plants is, in many cases, a limitation when determining their effects. It would be adequate to do the analyses with this information to have a complete picture.

**Author Contributions:** Conceptualization, X.Á. and E.V.; Data curation, C.A.-A.; Formal analysis, N.d.l.T.-R. and C.A.-A.; Funding acquisition, E.V.; Investigation, X.Á. and E.V.; Methodology, X.Á. and N.d.l.T.-R.; Project administration, X.Á.; Software, N.d.l.T.-R. and C.A.-A.; Supervision, X.Á.; Validation, X.Á. and E.V.; Writing—original draft, N.T.-R. and C.A.-A.; Writing—review and editing, X.Á. and E.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by TASGA. This work was possible thanks to the projects developed by the group AF4 for the company TASGA Renovables. Roberto Pérez Lodos and Juan Jesus Berzosa Aránguez from TASGA facilitated the work.

**Acknowledgments:** The authors are grateful to all the technicians participating in the samplings, as well as to Paula Rivas for her help.

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

### **References**


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

### *Article* **Evaluation of the Stability and Suitable Scale of an Oasis Irrigation District in Northwest China**

**Xifeng Zhang <sup>1</sup> , Yifan Zhang <sup>2</sup> , Jinghui Qi <sup>3</sup> and Qiang Wang 4,\***


Received: 28 August 2020; Accepted: 5 October 2020; Published: 13 October 2020

**Abstract:** Oases support human activities in arid and semiarid regions, and their stability is important for regional sustainable development and water resource management. Water consumption is the major factor affecting the stability of oases. On the basis of remote sensing images, evaporation and socioeconomic data, this study first evaluates the stability of the Dunhuang Oasis against the expansion of an oasis irrigation district and planting structure changes from 1987 to 2015. Next, it calculates a suitable area of the oasis irrigation district using water–energy balance theory. The results are as follows: (1) During the 1987–2015 period, with the expansion in the oasis irrigation area, the planting structure underwent a marked transformation from food crops to cash crops to orchards. Water consumption pattern likewise changed considerably. (2) The stability of the Dunhuang Oasis continued to weaken from 0.54 in 1987 until it reached a dangerously unstable level of 0.17 in 2010. With the implementation of water-saving measures and a water-transfer project, the stability of the Dunhuang Oasis irrigation district increased to a metastable level of 0.22 in 2015. (3) Setting the stability are 0.5 of a stable level and 0.75 of an extremely stable level, and the oasis irrigation district should be impractical and reduced by 168 and 241 km<sup>2</sup> to attain a suitable oasis ecosystem scale. Hence, at present, the water-transfer project is the most practical way to increase allocated water resource to the oasis irrigation district for improving its stability.

**Keywords:** oasis irrigation district; stability evaluation; suitable scale

### **1. Introduction**

The mountain–oasis–desert unit is the main topographic feature of inland river watersheds in China [1]. For watersheds in this arid region, water from melted snow and glacier in the mountain flows through an oasis similar to a water tower, thereby providing fresh water to people and nature before disappearing into the desert. Through this morphotectonic pattern, as an important landscape, oases support human activities and economic development in arid regions [2] but are generally water deficient [3]. Its stability is directly related to the sustainable development of regional economy and ecological security [2].

In recent decades, the expansion of the irrigation district, population growth, and economic development in the middle-stream oasis of the inland river basin has considerably modified hydrological cycles and reduced surface runoffs to downstream reaches. Consequently, this situation has led to dried-up rivers, reduced outflows to terminal lakes, and the deterioration of the ecosystem in downstream areas [4–7]. These effects have manifested in dried-up inland lakes, such as the Aral Sea in Central Asia [8,9], Lake Urmia in Iran [10,11], and the Shiyang River in Northwestern China [12]. Thus, assessing whether an oasis is stable under the present development pattern is important in the pursuit of sustainable development. Moreover, evaluating a suitable oasis scale is crucial to the government development plan.

Several studies on stability evaluation and suitable oasis scales have focused on the Endorheic Basin in Northwest of China, specifically, the Keriya River Basin [13], the Manas River Basin [2], the Tarim River Basin [14], the Weigan River [15], and the Heihe River Basin [16]. Research methods in stability evaluation and suitable scales have continuously progressed owing to the rapid development of ecological hydrology. In the past, suitable oasis cropland areas were calculated by establishing a regression equation using the total amount of runoff water resources and the cropland area of the oasis [17]. Wang et al. [18] applied contrastive analysis to different oasis landscapes and their homologous water–energy balance relationship and proposed the concept of the "green degree" to assess stability evaluation and suitable oasis scales. In addition, a suitable farmland scale model had been established from the perspective of crop water footprint and water resources available in an oasis city [19]. Recently, Hao et al. [16] developed an approach for calculating oasis and cultivated land scales by combining water–energy balance and wind–sand dynamic theories with ecological health assessments in the Heihe River Basin.

Agricultural irrigation consumes the largest proportion of the water supply of an oasis, and the suitability of a cropland scale can directly affect the stability of an oasis. Previous research on stability evaluation and suitable oasis scales has generally neglected the impact of changes in planting structure and agricultural water consumption. Additionally, no relevant studies have been conducted on the Shule River Basin, which is adjacent to the Heihe River Basin. The Dunhuang Oasis, which is the study area of this research, is located west of the Shule River Basin and an important region along the Silk Road.

Similar to the case of the middle reaches of the Heihe River Basin, the Dunhuang Oasis also has a sharp contradiction of the water resource issue between the economy and the ecosystem. In the past decade, planting structures constantly changed with the expansion in oasis irrigation areas. Moreover, water exploitation and utilization rates reached nearly 100% in the oasis, in which agricultural water consumption accounted for nearly 90% of total available resources [20]. Hence, natural oasis ecosystems are unable to receive necessary water resources. Dunhuang City has proposed a water resource plan, namely, the "Comprehensive Planning of the Rational Use of Water Resource and Protection of Ecosystem Services in the Dunhuang Basin" [21] to improve the ecological environment of its natural oasis. The main purpose of this plan is to allocate large amounts of water resources to the natural oasis in the lower reaches of the Dunhuang Basin and leave only a limited amount of water for the irrigation district of the Dunhuang Oasis.

The objective of this study is to propose a novel research idea to evaluate the stability of the oasis and calculate the suitable oasis scale by combining water–energy balance with altered planting structures. First, this study evaluates stability during the long-term expansion process of the Dunhuang Oasis from 1987 to 2015. Next, it calculates a suitable oasis scale to a certain stability degree. Finally, it discusses to what extent the current water plan can improve the stability of the oasis and provides a scientific reference for the sustainable development of the oasis.

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

### *2.1. Study Area*

The Dunhuang Oasis irrigation district is located east of Dunhuang City (92◦130–95◦300 N, 39◦400–41◦400 E) in Northwest China, with a total area less than 5% of that of the entire city (Figure 1). The Dunhuang Oasis irrigation district lies between the Mingsha Mountain (Echoing-sand Mountain) and the Sanwei Mountain in the southeast and the Mazong Mountain in the north and spreads into the

Gobi Desert in the west, thereby earning the name the "Gobi Desert Oasis". The world famous Mogao Grottoes are located southwest of the Dunhuang Oasis. The elevation of the oasis ranges from 800 to 1500 m and is high in the south and low in the north. The oasis experiences the strong sunshine with the annual solar radiation reaches to 6418.4 MJ/m<sup>2</sup> , and mean maximum and minimum temperatures are −28.6–43.6 ◦C. The oasis also experiences high evaporation, with the annual potential evaporation reaches 2486.0 mm, but the annual precipitation is only 39.2 mm. The Danghe River, which originates from the Qilian Mountain, provides the water resources for the oasis, and the Danghe River alluvial fan developed the oasis. The oasis consists mainly of gray brown desert, alkali saline, and irrigation desert soil (China soil science database, http://vdb3.soil.csdb.cn). According to the Dunhuang City National Economic Statistics Yearbook (1987–2015), spring wheat, melon, cotton, and grapes are the main crops planted in this area [22]. *Water* **2020**, *12*, x FOR PEER REVIEW 3 of 12 the annual solar radiation reaches to 6418.4 MJ/m2, and mean maximum and minimum temperatures are −28.6–43.6 °C. The oasis also experiences high evaporation, with the annual potential evaporation reaches 2486.0 mm, but the annual precipitation is only 39.2 mm. The Danghe River, which originates from the Qilian Mountain, provides the water resources for the oasis, and the Danghe River alluvial fan developed the oasis. The oasis consists mainly of gray brown desert, alkali saline, and irrigation desert soil (China soil science database, http://vdb3.soil.csdb.cn). According to the Dunhuang City National Economic Statistics Yearbook (1987–2015), spring wheat, melon, cotton, and grapes are the main crops planted in this area [22].

**Figure 1.** Study area derived by Landsat TM imagery (Bands 5-4-3). **Figure 1.** Study area derived by Landsat TM imagery (Bands 5-4-3).

### *2.2. Remote Sensing Image Data Processing 2.2. Remote Sensing Image Data Processing*

The high-quality remote images with a 30 m spatial resolution from Landsat TM (Thematic Mapper) were taken in low cloud covers during high biomass season, and all mage data were acquired from http://glovis.usgs.gov/ and http://www.radi.ac.cn/. We selected 1987, 1990, 1996, 2001, 2007, 2010, and 2015 to identify the actual irrigation district of the oasis using the software of ERDAS 9.1 and ArcGIS 10.3. Based on a 1:250,000 topographic map of the Dunhuang Oasis, all images were processed under the common universal transverse mercator coordinate system. A total of 50 uniformly distributed ground control points (e.g., roads and rivers) were used for geometric correction and georeferencing, and we used the quadratic polynomial transformation and nearestneighbor resampling methods to identify ground control points in image-to-map rectification, the root mean square error of the geometrical rectification was less than one pixel. When the images were all ready, we divided the land use types into eight classes according to visual interpretation, namely, cropland, water, high-density grassland, medium-density grassland, low-density grassland, shrub land, urban construction land, and barren land. For verifying the result of visual interpretation, the field investigation points and corresponding interpretation results were compared, and the overall interpretation accuracy of land use classification in the 1987, 1990, 1996, 2001, 2007, 2010, and 2015 images was 80.91%, 84.96%, 79.85%, 79.69%, 84.98%, 85.79%, and 89.45%, respectively, which met the minimum accuracy requirement of 70% [23]. In this study, we used only the results of the cropland The high-quality remote images with a 30 m spatial resolution from Landsat TM (Thematic Mapper) were taken in low cloud covers during high biomass season, and all mage data were acquired from http://glovis.usgs.gov/ and http://www.radi.ac.cn/. We selected 1987, 1990, 1996, 2001, 2007, 2010, and 2015 to identify the actual irrigation district of the oasis using the software of ERDAS 9.1 and ArcGIS 10.3. Based on a 1:250,000 topographic map of the Dunhuang Oasis, all images were processed under the common universal transverse mercator coordinate system. A total of 50 uniformly distributed ground control points (e.g., roads and rivers) were used for geometric correction and georeferencing, and we used the quadratic polynomial transformation and nearest-neighbor resampling methods to identify ground control points in image-to-map rectification, the root mean square error of the geometrical rectification was less than one pixel. When the images were all ready, we divided the land use types into eight classes according to visual interpretation, namely, cropland, water, high-density grassland, medium-density grassland, low-density grassland, shrub land, urban construction land, and barren land. For verifying the result of visual interpretation, the field investigation points and corresponding interpretation results were compared, and the overall interpretation accuracy of land use classification in the 1987, 1990, 1996, 2001, 2007, 2010, and 2015 images was 80.91%, 84.96%, 79.85%, 79.69%, 84.98%, 85.79%, and 89.45%, respectively, which met the minimum accuracy requirement of 70% [23]. In this study, we used only the results of the cropland in the interpretation, as the scale of the

in the interpretation, as the scale of the oasis was mainly affected by the changes in the cropland. On the basis of the interpreted cropland data, we calculated the actual crop areas of crop food, cash food, oasis was mainly affected by the changes in the cropland. On the basis of the interpreted cropland data, we calculated the actual crop areas of crop food, cash food, and orchards by multiplying the crop food, cash food, and orchard area ratios in the Dunhuang City statistical yearbook [22] to determine the water consumption of each crop and orchard.

### *2.3. Water Consumption Data Analyses*

In this study, water consumption denoted crop, domestic, and industrial water consumption.

### 2.3.1. Crop Water Consumption

The Food and Agricultural Organization (FAO)–Penman–Monteith methods were used to estimate reference crop evapotranspiration (ET0), and multiplying by the crop coefficient (Kc) can get actual evapotranspiration (ET) as crop water consumption [24]:

$$\text{ET}\_0 = \frac{0.408\Delta(\text{R}\_\text{n} - \text{G}) + \gamma \frac{900}{\text{T} + 2\text{T}3} \text{U}\_{2(\text{e}\_\text{s} - \text{e}\_\text{s})}}{\Delta + \gamma (1 + 0.34\text{U}\_2)} \tag{1}$$

R<sup>n</sup> is the net radiation at the canopy surface (MJ/m<sup>2</sup> . day);

G is the soil heat flux density (MJ/m<sup>2</sup> . day);

T is the mean daily air temperature at 2 m above the ground (◦C);

U<sup>2</sup> is the wind speed at 2 m above the ground (m/s);

e<sup>s</sup> is the saturation vapor pressure (kPa);

e<sup>a</sup> is the actual vapor pressure (kPa);

e<sup>s</sup> − e<sup>a</sup> is the saturation vapor pressure deficit (kPa);

∆ is the slope of the vapor pressure temperature relationship (kPa/ ◦C);

γ is the psychometric constant (kPa / ◦C).

Based on the above principle, ET<sup>0</sup> is computed using the CropWat 8.0 [25] tool after the monthly averages of minimum temperature, maximum temperature, humidity, wind speed, and sunshine hours are inputted.

The values of Kc vary with the crop development stages, and values of monthly Kc were adopted from FAO-56 [24] for spring wheat, cotton and grape. Adjustments to Kc mid in climates where RHmin differs from 45% or wher *U*<sup>2</sup> is larger or smaller than 2.0 m/s, were made by following the guidelines of Allen et al. [24]:

$$\mathbf{K}\_{\rm cmd} = \mathbf{K}\_{\rm cmd(Tab)} + [0.04(\mathbf{U}\_2 - 2) - 0.004(\mathbf{RH}\_{\rm min} - 45)] \left(\frac{\mathbf{h}}{3}\right)^{0.3} \tag{2}$$

$$\mathbf{K}\_{\rm cend} = \mathbf{K}\_{\rm cend(Tab)} + \left[ 0.04(\mathbf{U}\_2 - 2) - 0.004(\mathbf{RH}\_{\rm min} - 45) \right] \left( \frac{\mathbf{h}}{3} \right)^{0.3} \tag{3}$$

Kcmid(Tab) is the tabulated Kc values in the mid-seasonof Table VI-12 of Allen et al. [24]; Kcend(Tab) is the tabulated Kc values in the late-season of Table VI-12 of Allen et al. [24]; *U*<sup>2</sup> is wind speed at 2 m height over grass, the range is 1 m/ s ≤ *U*<sup>2</sup> ≤ 6 m m/s; RHmin is daily minimum relative humidity, the range is 20% ≤ RHmin ≤ 80%; H is mean plant height, the range is 0.1 m ≤ h ≤ 10 m.

This research considered the water consumption of vines. Figure 2 shows the adjusted K<sup>C</sup> of wheat, cotton and grape. Compared with wheat and cotton, the K<sup>C</sup> of grape is more complicated, and the values of the germination, growing season with shoots, flowering, and berry growing periods as well as berry and tendril maturation were 0.80, 1.09, 1.13, 1.07, 1.03, and 0.82, respectively.

**Figure 2.** Crop evapotranspiration coefficients. **Figure 2.** Crop evapotranspiration coefficients.

Actual ET rates of food crops and cash crops are estimated by Allen et al. [24] Actual ET rates of food crops and cash crops are estimated by Allen et al. [24]

$$\text{ET} = \text{ET}\_0 \times \text{K}\_\mathbb{C} \tag{4}$$

### 2.3.2. Domestic Water Consumption 2.3.2. Domestic Water Consumption

defined as follows:

Domestic water consumption represented the water consumption of urban and rural populations, tourists, and livestock, which was positively correlated with living standards. It is Domestic water consumption represented the water consumption of urban and rural populations, tourists, and livestock, which was positively correlated with living standards. It is defined as follows:

$$\mathbf{D} = (\mathbf{P}\_1 \times \mathbf{C}\_1 + \mathbf{P}\_2 \times \mathbf{C}\_2 + \mathbf{P}\_3 \times \mathbf{C}\_3 + \mathbf{P}\_4 \times \mathbf{C}\_4) \times 365 \times 10^{-8} \tag{5}$$

D = (Pଵ × Cଵ + Pଶ × Cଶ + Pଷ × Cଷ + Pସ × Cସ) × 365 × 10ି଼ (5) D is the total domestic water consumption (10<sup>8</sup> m<sup>3</sup> );

D is the total domestic water consumption (108 m3); P<sup>1</sup> is the amount urban population (10<sup>8</sup> m<sup>3</sup> );

P1 is the amount urban population (108 m3); P<sup>2</sup> is the amount rural population (10<sup>8</sup> m<sup>3</sup> );

P2 is the amount rural population (108 m3); P<sup>3</sup> is the amount tourists population (10<sup>8</sup> m<sup>3</sup> );

P3 is the amount tourists population (108 m3); P<sup>4</sup> is the amount livestock number (10<sup>8</sup> m<sup>3</sup> , sheep unit);

P4 is the amount livestock number (108 m3, sheep unit); C<sup>1</sup> is the average per capital water use coefficient of urban (L/day);

C1 is the average per capital water use coefficient of urban (L/day); C2 is the average per capital water use coefficient of rural (L/day); C<sup>2</sup> is the average per capital water use coefficient of rural (L/day);

C3 is the average per capital water use coefficient of tourist (L/day); C<sup>3</sup> is the average per capital water use coefficient of tourist (L/day);

C4 is the average per capital water use coefficient of livestock (L/day); C<sup>4</sup> is the average per capital water use coefficient of livestock (L/day);

The values for C1, C2, C3 and C4 (Table 1) are obtained from the literature for the four periods with different living standards: 1987–1990, 1991–2000, 2001–2007, and 2008–2015 [26]. The values for C1, C2, C<sup>3</sup> and C<sup>4</sup> (Table 1) are obtained from the literature for the four periods with different living standards: 1987–1990, 1991–2000, 2001–2007, and 2008–2015 [26].


C5 is the amount of water used to create 10,000 RMB worth of industrial output (m3/104 RMB).

The value for C5 (Table 1) was also determined separately for the four periods with different living

### 2.3.3. Industrial Water Consumption

Industrial water consumption was affected by industrial output, technology, processes, and the amount of water used to create RMB 10,000 (Chinese currency) worth of industrial output. It is defined as follows:

$$P = \text{Indo} \times \text{C}\_5 \times 10^{-8} \tag{6}$$

P is the industry water use (10<sup>8</sup> m<sup>3</sup> );

Indo is the industrial output (RMB);

C<sup>5</sup> is the amount of water used to create 10,000 RMB worth of industrial output (m<sup>3</sup> /10<sup>4</sup> RMB).

The value for C<sup>5</sup> (Table 1) was also determined separately for the four periods with different living standards for the period of 1987–1990, 1991–2000, 2001–2007, and 2008–2015, respectively [26].

### *2.4. Oasis Stability and Suitable Oasis Scale Model*

In this study, a stability index (H0), which is based on the theory of ecological water–energy balance, was used to estimate the stability degree of the oasis under certain water resource conditions [18]. H<sup>0</sup> is the relative equilibrium analysis between water and energy conditions of the oasis, which can not only reflect ecological evolution and degradation from the internal perspective of the oasis but also the "green degree" or "stability degree" from the overall regional view point. The greater the H0, the less affected the water stress and the higher the oasis stability, and vice versa. The formula is as follows:

$$\mathbf{H}\_0 = \frac{\mathbf{W}\_1 - \mathbf{W}\_2 + \mathbf{P} \times \sum \mathbf{A}\_i}{\mathbf{E} \mathbf{T}\_0 \times \sum \mathbf{A}\_i} \tag{7}$$

where A<sup>i</sup> is the area of the land type (I = food crops, cash crops, and grapes; km<sup>2</sup> ), W<sup>1</sup> is the total available water volume of the river basin (10<sup>8</sup> m<sup>3</sup> ), W<sup>2</sup> is the annual average agricultural, industrial, and domestic water consumption, ET<sup>0</sup> is the reference crop evapotranspiration, and P is the annual average precipitation. The H<sup>0</sup> of the oasis was classified on the basis of the characteristics of the natural environment (Table 2).


**Table 2.** Classification of oasis stability.

Based on the previous section, the calculation model of the suitable oasis scale (A) is

$$\mathbf{A} = \frac{\mathbf{W}\_1 - \mathbf{W}\_2}{(\mathbf{ET}\_0 - \mathbf{P}) \times \mathbf{k}\_\mathbf{P} \times \mathbf{H}\_0^\*} \tag{8}$$

where k<sup>p</sup> is the comprehensive coefficient of plants in the oasis, which includes planting crops, trees, and grass. H∗ 0 is used to estimate the suitable oasis scale.

### **3. Results**

### *3.1. Land Use*/*Land Cover Changes between 1987 and 2015*

Figures 3 and 4 show changes in land use/land cover between 1987 and 2015. From Figure 3, we can deduce that the cropland area exhibited an expanding pattern, and the interpretation results demonstrate the total cropland area in the oasis was 272.69, 276.41, 295.33, 327.08, 371.40, 380.47, and

389.41 km<sup>2</sup> in 1987, 1990, 1996, 2000, 2007, 2010, and 2015, respectively. In terms of the temporal characteristics of land use/land cover changes, incremental rates from 1996 to 2007 more than doubled from 1987 to 1996 and from 2007 to 2015. With regard to spatial pattern, the cropland increased by 116.72 km<sup>2</sup> . Growth is mainly attributed to its transformation from a grassland and a barren land in the fringe of the oasis [27] and to the continuously expanding urban construction land within the oasis. *Water* **2020**, *12*, x FOR PEER REVIEW 7 of 12 the fringe of the oasis [27] and to the continuously expanding urban construction land within the oasis. *Water* **2020**, *12*, x FOR PEER REVIEW 7 of 12 the fringe of the oasis [27] and to the continuously expanding urban construction land within the oasis.

**Figure 3.** Land use/land cover spatial pattern from 1987 to 2015 in the Dunhuang Oasis. **Figure 3.** Land use/land cover spatial pattern from 1987 to 2015 in the Dunhuang Oasis. **Figure 3.** Land use/land cover spatial pattern from 1987 to 2015 in the Dunhuang Oasis.

**Figure 4.** Planting structure changes in the Dunhuang Oasis. **Figure 4.** Planting structure changes in the Dunhuang Oasis.

**Figure 4.** Planting structure changes in the Dunhuang Oasis. The three different-colored bar charts in Figure 4 represent the crop areas of crop food, cash food, and orchards. The planting structure exhibited marked changes and even transformed during the 1987–2015 period. Specifically, (1) the crop food area decreased continuously from 1987 to 2007, and the crop food area in 2007 was only 1.61km2, making it too small to display. Meanwhile, cotton was the main cash crop, which experienced a substantial increase owing to its high value. A distinct transformation in the food and cash crops from 1996 to 2000 can be observed, as the crop area of the cash crops gradually became larger than that of the food crops. In addition, the trial planting of grapes The three different-colored bar charts in Figure 4 represent the crop areas of crop food, cash food, and orchards. The planting structure exhibited marked changes and even transformed during the 1987–2015 period. Specifically, (1) the crop food area decreased continuously from 1987 to 2007, and the crop food area in 2007 was only 1.61km2, making it too small to display. Meanwhile, cotton was the main cash crop, which experienced a substantial increase owing to its high value. A distinct transformation in the food and cash crops from 1996 to 2000 can be observed, as the crop area of the cash crops gradually became larger than that of the food crops. In addition, the trial planting of grapes The three different-colored bar charts in Figure 4 represent the crop areas of crop food, cash food, and orchards. The planting structure exhibited marked changes and even transformed during the 1987–2015 period. Specifically, (1) the crop food area decreased continuously from 1987 to 2007, and the crop food area in 2007 was only 1.61km<sup>2</sup> , making it too small to display. Meanwhile, cotton was the main cash crop, which experienced a substantial increase owing to its high value. A distinct transformation in the food and cash crops from 1996 to 2000 can be observed, as the crop area of the cash crops gradually became larger than that of the food crops. In addition, the trial planting of grapes in the sandy soil region began in 1996. (2) In 2010, cotton areas gradually decreased, and grape orchards increasingly involved large-scale cultivation.

### *3.2. Water Consumption*

### 3.2.1. Agricultural Water Consumption

Agricultural water consumption comprises the core section of the suitability evaluation of the oasis. In the Dunhuang Oasis, agricultural water consumption includes food crops, cash crops, and orchards. From 1987 to 2015, total water consumption increased initially from 2.293 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> in 1987 to 3.513 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> in 2007, then decreased to 2.902 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> in 2015 (Table 3). Listed individually in Table 3, food crop water consumption decreased sharply by 1.334 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> from 1987 to 2015. Meanwhile the water consumption of cash crops with high economic benefits progressively increased by 2.374 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> from 1987 to 2007 then began to decrease to 2.207 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> in 2010 and eventually decrease to a very low value of 0.774 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> in 2015. This trend in cash crop water consumption is closely related to that of grapes, which increased slowly before 2010 before rising rapidly.


**Table 3.** Agricultural water consumption in 1987–2015 (10<sup>8</sup> m<sup>3</sup> ).

### 3.2.2. Domestic and Industrial Water Consumption

Domestic and industrial water consumption should not be neglected in the suitability evaluation of the oasis. During the 1987–2015 period, the total population of the Dunhuang Oasis increased from 108,373 to 142,558, in which the rural population accounted for nearly 70% to 80%. Rural livestock increased from 349,820 sheep units in 1987 to 486,816 in 2015. In addition, the number of tourists increased from less than 100,000 in 1987 to 6,603,914 in 2015. Furthermore, industrial output increased by 84.6 <sup>×</sup> <sup>10</sup><sup>8</sup> RMB [22]. Thus, we can calculate overall domestic and industrial water consumption under increasing populations and the booming tourism industry. Table 4 shows that overall domestic and industrial water consumption rapidly increased from 0.038 <sup>×</sup> <sup>10</sup><sup>8</sup> to 0.219 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> . The percentage of domestic and industrial water consumption in overall water consumption is very small.

**Table 4.** Domestic and industrial water consumption in 1987–2015 (10<sup>8</sup> m<sup>3</sup> ).


### *3.3. Oasis Stability Evaluation*

In this study, the total available water volume originates from the Danghe River; thus, we use its perennial mean runoff, that is, 4.13 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> , as W1. From the above data, total available water quantity for the oasis in 1987, 1990, 1996, 2000, 2007, 2010, and 2015 was 1.79 <sup>×</sup> <sup>10</sup><sup>8</sup> , 1.66 <sup>×</sup> <sup>10</sup><sup>8</sup> , 1.37 <sup>×</sup> <sup>10</sup><sup>8</sup> , 1.78 <sup>×</sup> <sup>10</sup><sup>8</sup> , 0.43 <sup>×</sup> <sup>10</sup><sup>8</sup> , 0.67 <sup>×</sup> <sup>10</sup><sup>8</sup> , and 1.00 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> , respectively. From the oasis stability evaluation, as shown in Equation (7), the H<sup>0</sup> of the Dunhuang Oasis were 0.54, 0.51, 0.41, 0.39, 0.15, 0.19, and 0.22 from 1987 to 2015 (Table 5).


**Table 5.** Stability of Dunhuang Oasis.

### *3.4. Suitable Oasis Irrigation District Scale*

We take 2015 as an example and use Equation (8) to derive a suitable oasis irrigation district scale. Crops and natural vegetation are considered as having a common effect on kp, which is 0.72 [18], and the H∗ 0 is set as two values, namely, 0.5 for the stable level and 0.75 for the extremely stable level.

Table 6 shows that the suitable oasis irrigation district scale was smaller than the actual area in 2015. According to the water–energy balance model, the current oasis irrigation district scale should be reduced by 168 and 241 km<sup>2</sup> to attain stable to extremely stable levels, respectively.

**Table 6.** Suitable oasis irrigation district scales/km<sup>2</sup> .


### **4. Discussion**

The planting structure exhibited two marked changes. The first planting structure change happens between crop food and cotton. Except for cotton's high value and cultivation suitability, the rapidly increasing cotton lands, which were mainly attributed to farmers, had autonomy in terms of land use activities since the early 1980s under China's economic reform policy. Against this background, farmers envisioned market economy ideas, and production activities were closely associated with market demands. In addition, the trial planting of grapes in the sandy soil region began in 1996. The second planting structure change happens between cotton and grape; this is because of successful grape trial experiments in the sandy soil region, and from then on, a large number of farmers began to install grape trellises in fields previously planted with cotton or wheat. Hence, the planting structure of the oasis changed once again.

Although the total cropland area increased by 18.01 km<sup>2</sup> from 2007 to 2015, agricultural water consumption decreased by 0.611 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> . This result may be beneficial to the transformation of crop patterns and water-saving irrigation measures. Traditional field canal irrigation was the primary irrigation pattern used in the past, but, from the beginning of 2010, advanced water-saving irrigation models, such as micro, pipe, and greenhouse micro irrigation, were implemented in the entire oasis.

Table 5 shows that stability was in a stable level in 1987 and 1990 owing to low agricultural, domestic, and industrial water consumption. During this period, though each industry was gradually developing, stability was not extremely stable, which may have been due to limited water resources. Farmers' enthusiasm for production and cropland areas increased under reform and opening-up policies, but irrigation measures were difficult. At the same time, tourism in the Dunhuang Oasis began to flourish. When the national economy improved, the stability of the oasis fell to a metastable level from 1996 to 2000 and reached a dangerously unstable level from 2007 to 2010. Water resource exploitation and utilization rates nearly reached 100% in the oasis, in which agricultural water consumption accounted for nearly 90% of overall available resources [27]. Serious ecological problems, such as accelerated desertification and salinization, shrunken terminal lake and declining groundwater

level, accompany rapid economic development [27]. Water has become the primary restricting and bottleneck factor in the socioeconomic development of the Dunhuang Oasis.

In this case, a series of water resource plans was implemented to address this issue, in which the most important program "Comprehensive Planning of the Rational Use of Water Resource and Protection of Ecosystem Services in the Dunhuang Basin" was proposed. This program aims to reduce the croplands of state farms, implement agricultural water-saving measures, and carry out a water diversion project from the Sugan Lank Basin to the Danghe River. Approximately 0.835 <sup>×</sup> <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> of water allocated from the water diversion project is intended for the improvement of the ecological environment, which plays a crucial role in alleviating the water crisis by increasing groundwater levels and recovering the vegetation in marginal areas of the oasis [20]. The system dynamic model simulated the agricultural water consumption under different scenarios in the Dunhuang Oasis and shows that the proportion of agricultural water consumption in overall water consumption can be reduced from 92.50% in 2010 to 86.30% in 2025 [20], but, if reduced by 168 km<sup>2</sup> , to attain a stable level, the agricultural water consumption should be decreased to at least half of what it is now. In this study, the suitable oasis area is far less than the actual area, which is a common issue in the Endorheic watershed oasis in Northwest China [2,16,28]. Reducing cropland in the oasis is the most direct and effective means, which is difficult to achieve. Specifically, individual croplands in the Dunhuang Oasis should be reduced to preserve the ecological environment. However, several reasons highlight the difficulty of this solution. (1) The reduction of croplands will harm the economic interests of farmers and subsequently decrease their quality of life. Thus, the possibility of criminal problems due to poverty should be considered. (2) The reduction of croplands is not in line with the Chinese policy of farmland protection. Hence, only on the basis of maintaining water-saving irrigation, reducing domestic water consumption, improving industrial water consumption efficiency, forbidding sprawl inside and outside the oasis, and increasing the amount of water allocated to the oasis from the water-transfer project, can the stability of the oasis be improved and the sustainable development of the regional economy and ecology be maintained.

### **5. Conclusions**

This study analyzed the stability of the Dunhuang Oasis against the background of planting structure changes during the 1987–2015 period. Our main findings and recommendation are as follows:


**Author Contributions:** X.Z. and Q.W. conceived and designed the experiments. X.Z. and Y.Z. performed the experiments and analyzed the data. J.Q. provided remote sensing image of 2010 and 2015. X.Z. and Q.W. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Science Foundation of China, grant number: 41701321 and No. 41661105.

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

### **References**


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

**Yaohuan Huang 1,2, Zhonghua Li <sup>3</sup> , Mingxing Chen 2,4,\* , Xiaoyang Song <sup>5</sup> and Ping Kang <sup>6</sup>**

	- Resources Research, Chinese Academy of Sciences, Beijing 100101, China

**Abstract:** Water resource has become a key constraint for implementing the "Belt and Road" initiative which was raised by the Chinese government. Besides the study of spatial and temporal variability of precipitation, this study created a water hazard risk map along the "Belt and Road" zone through combined flood and drought data from 1985. Our results showed that South-Eastern Asia, southern China and eastern Southern Asia are areas with the most abundant precipitations, while floods in these areas are also the most serious. Northwest China, Western Asia, Northern Africa and Southern Asia are areas highly vulnerable to drought. Furthermore, the potential influence of flood and drought were also analyzed by associating with population distribution and corridor map. It reveals that China, South-Eastern Asia, Southern Asia, Western Asia and Northern Africa have the largest population number facing potential high water hazard risk. China–India–Burma Corridor and China–Indo-China Peninsula Corridor have the largest areas facing potential high water hazard risk.

**Keywords:** water security; precipitation; drought; flood; transport planning

### **1. Introduction**

In 2013, China proposed the "Silk Road Economic Belt" and the "21st-Century Maritime Silk Road" initiatives, which are collectively referred to as the "Belt and Road" [1]. The "Vision and Action for Promoting the Construction of the Silk Road Economic Belt and the 21st Century Maritime Silk Road" ("Vision and Action") was issued subsequently in 2015, which marks the formal implementation of the "Belt and Road" [2]. "Vision and Action" proposes that infrastructure interconnection is a priority mission for the "Belt and Road", due to improving the accessibility of roads not only facilitating the lives of local individuals but also promoting the sustainable development of the local economy and society [3]. Strengthening the ecological cooperation of the countries and regions along the corridor to avoid potential ecological risks and establish a green silk road is another important recommendation of the "Vision and Action" [2]. It is foreseeable that the implementation of the "Belt and Road" will have a significant impact on the transportation, urbanization and water resources in the countries and regions along the corridor. However, water scarce, flood and drought have posed a huge potential threat to the implement of the "Belt and Road" and sustainable development of society, especially in arid and semi-arid areas [4–7]. Therefore, understanding the temporal and spatial variability of water resources and water-caused natural hazards along the corridor can not only ensure the smooth construction of transportation infrastructure but also promote the sustainable development of countries and regions along the "Belt and Road" zone [8].

**Citation:** Huang, Y.; Li, Z.; Chen, M.; Song, X.; Kang, P. Spatial Variability of Water Resources State of Regions around the "Belt and Road". *Water* **2021**, *13*, 2102. https://doi.org/ 10.3390/w13152102

Academic Editors: Fernando António Leal Pacheco and Luís Filipe Sanches Fernandes

Received: 5 June 2021 Accepted: 27 July 2021 Published: 31 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/).

The "Belt and Road" zone traverses Eurasia and spans subtropical, temperate, cold temperate and frigid zones, with complex terrain and geological conditions. It is a high-risk zone for natural hazards such as flood, drought and extreme precipitation [9]. In Central and Western Asia, the construction of transportation, oil and gas pipelines is threatened by high temperatures, droughts and extreme precipitation. Frequent floods in Southern and South-Eastern Asia have also brought tremendous security risks to the operation of transportation facilities. According to the EM-DAT hazard database, there were more than 7200 natural hazards that happened in the world from 1990 to 2010 and more than 3003 times in the "Belt and Road" area, including 1131 floods and 94 droughts respectively which caused serious casualties and property damage [10]. At the same time, most countries and regions along the "Belt and Road" corridor are economically underdeveloped and agricultural dependent with relatively weak ability to withstand water stress [11]. Impacted by global warming, the risks of encounter extreme precipitation, drought and flooding are increasing [12–15], which brings tremendous potential risk to local human life, property, agriculture, economy and society along the "Belt and Road" corridor [16–19]. The fifth assessment report of Intergovernmental Panel on Climate Change (IPCC) indicated that the Lancang-Mekong River Basin has an increased precipitation during the monsoons of the past 30–50 years, while the precipitation during the dry season has dropped sharply, and with the accelerated thaw of the Himalayan glaciers and the Pamirs glaciers, the supply of glaciers to the rivers will be significantly decreased, and in a few years the billions of people living in South and Central Asia may confront the risk of losing fresh water [20]. Xia et al. found that there will probably be an increase in extreme floods and droughts in the Eastern Monsoon of China and irrigation water in the North China Plain will increase by 4% with the impact of global warming [21]. Reza et al. analyzed the data observed by more than 400 river monitoring stations distributed throughout Iran and discovered that floods caused by extreme precipitation in most parts of Iran have an obvious growth [22]. Other studies have shown that Belt and Road countries or regions, such as the Philippines, and Vietnam, Pakistan, the North-South Road Corridor and East-West Road Corridor of Myanmar, South Asia and South-Eastern Asia, also faced serious risk of floods and droughts [23–25]. However, spatial variability of water resources' condition and its related risk at a macro scale is more significant to the sustainability of Belt and Road Initiative.

Since water resources is a major constraint to the execution of "Belt and Road" initiative, it is necessary to figure out the spatio-temporal pattern of water resources on a macro scale. In this paper, we first analyzed the spatial distribution of mean precipitation in more than 60 countries of the "Belt and Road" zone. Then we used the flood data and drought data in the same period to make the water hazard risk map and graded the map with five levels based on potential water-caused risk. In addition, based on the flood and drought frequency map, we assessed the potential impact population, the corridor along the "Belt and Road" zone. We hope this paper will provide some support of water resources' state for the Belt and Road initiative to some extent.

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

### *2.1. Study Area*

The "Belt and Road" zone stretches across the continent of Asia, Europe and Africa from the Pacific in the East to the Atlantic Ocean in the west and from Indonesia in the South to the Arctic Ocean in the north (Figure 1). Over 60% of areas of the "Belt and Road" zone is arid and semi-arid grassland, desert and high-altitude ecologically fragile areas with dry climate and low precipitation caused by effect of the Himalayas and global weather patterns [4]. Central Asia, Western Asia and Northern Africa are the driest areas in the world where serious water shortage and severe land desertification pose serious threats to social and economic sustainable development. South-Eastern Asia and Southern Asia are strongly affected by monsoon, with frequent natural disasters including droughts, torrential rains and floods [26]. "Belt and Road" countries have a population amounting to more than 70% of the world's population. However, the amount

of water resources is only 36% that of the global total. It poses higher pressure in water security compared with the world average level [27]. In the initial vision for "Belt and Road Initiative" in 2015, there are 6 land economic and transportation corridors, which are a global economic connectivity program led by China. The 6 land corridors are China– Mongolia–Russia Corridor, New Eurasian Continental Bridge, China-Central Asia-West Asia Corridor, China–Pakistan Corridor, Bangladesh–China–India–Burma Corridor and China-Indochina Peninsula Corridor (Figure 1). Herewith, we zoned the study area of "Belt and Road" to 6 major zones of the Mongolia-Russia and Central Asia zone (MRCA), the South-Eastern Asia zone, the Southern Asia zone, the Western Asia and Northern Africa zone (WANA), the Central and Eastern Europe zone (CE Europe) and China [28] according to the 6 corridors.

**Figure 1.** "Belt and Road" monitoring area and corridor map.

### *2.2. Data Source*

The data used in this paper mainly include precipitation, drought, flood, population, cropland, highway and railway. The specific data sources and formats are shown in the following table (Table 1).



(1) Precipitation

Global daily precipitation data of 1985-2016 is derived from the NOAA Climate Prediction Center (CPC) Unified Precipitation Products dataset. It is created on a 0.5◦ lat/lon over the global land by interpolating gauge observations from 30,000 stations by considering orographic effects in precipitation.

(2) Drought

The Global SPEI (Standardized Precipitation Evapotranspiration Index) drought dataset of 1985–2016 is made available by Consejo Superior de Investigaciones Científicas (CSIC), with a 0.5 degrees spatial resolution and a monthly time resolution. SPEI is one of the most widely used drought indices in monitoring and quantifying droughts, which is developed by Vicente-Serrano et al. [29] based on the standardized precipitation index (SPI). The specific procedures for calculating SPEI can be found in the studies of Mahmoudi et al. [30] and Pei et al. [31] Positive values of SPEI indicate wet conditions, while negative values indicate dry conditions. SPEI classification rules are as follows: SPEI > −0.5 no drought; −1.0 < SPEI ≤ −0.5 light drought; −1.5 < SPEI ≤ −1.0 moderate drought; −2.0 < SPEI ≤ −1.5 severe drought; SPEI ≤ −2.0 especially severe drought. In this paper, the droughts severity of no drought, light drought, moderate drought, severe drought and especially severe drought is assigned to values of 0, 1, 2, 3, and 4, respectively. Then, the drought frequency was classified into 5 classes of approximately equal number of grid cells based on the accumulated 32 years of datasets of drought severity from 1985 to 2016.

(3) Flood

The flood dataset of 1985–2016 is derived from Dartmouth Flood Observatory, University of Colorado, with a 0.5 degrees spatial resolution and a yearly time resolution, which is mainly retrieved based on official reports and remote sensing sources (http: //floodobservatory.colorado.edu/index.html, accessed on 30 July 2021). The original data of floods is divided into three classes of severity based on 1–2 scale. Class 1 (large flood events): significant damage to structures or agriculture; fatalities; and/or 1–2 decadeslong reported interval since the last similar event. Class 2 (very large events): with a greater than 2 decades but less than 100 year estimated recurrence interval and/or a local recurrence interval of at 1–2 decades and affecting a large geographic region (> 5000 km<sup>2</sup> ). Class 3 (Extreme events): with an estimated recurrence interval greater than 100 years. In this paper, we assigned the Class 1, Class 2 and Class 3 to values of 1, 2 and 3, respectively, and generated the 5-classes flood frequency map similar to drought frequency map mentioned above.

(4) Population

The population data of 2015 is obtained from Socioeconomic Data and Applications Center of NASA with spatial resolution is 30<sup>00</sup> × 3000. Population input data are collected at the most detailed spatial resolution available from the results of the 2010 round of Population and Housing Censuses, which occurred between 2005 and 2015. The raster datasets are constructed from national or subnational input administrative units to which the estimates have been matched.

(5) Cropland

The Global Food Security-support Analysis Data 30 meter (GFSAD30) Cropland Extent data product provides cropland extent data across the globe, divided and distributed into 7 separate regional datasets, for nominal year 2015 at 30 meter resolution which was released by Department of the Interior, U.S. Geological Survey. Additionally, the validation dataset used to conduct an independent accuracy assessment of global cropland extent is available.

(6) Railway and Highway

The railway and highway data of 2016 was downloaded from the Environment Data Cloud Platform of Institute of Geographical Sciences and Natural Resources Research, Chinese Academy of Sciences. The data format is ARCGIS shapefile.

### **3. Analysis of Spatio-Temporal Heterogeneity of Water Security**

### *3.1. Precipitation*

Changes in precipitation is the primary driving forces of flash floods [32,33]. Similar to previous studies of precipitation in the regions of the Belt and Road [34–37], our analysis of mean precipitation along the "Belt and Road" zone from 1985 to 2016 shows that the spatial distribution of precipitation is extremely uneven. The Western Asia and Northern Africa zone has the least annual precipitation of 142 mm among all regions; the areas with the most abundant precipitation are mainly concentrated in South-Eastern Asia, eastern Southern Asia and southern China. In particular, South-Eastern Asia has the highest average precipitation of 1867 mm, much higher than the average level of "Belt and Road" zone (Figure 2, Table 2). Spatial distribution of precipitation within China is also extremely uneven. While there is abundant precipitation in the southeast coast, precipitation in Northwest China is scarce. Southern Asia and South-Eastern Asia are faced with similar situations as China. Precipitation in the Mongolia-Russia and Central Asia and the Central and Eastern Europe are both at moderate level along the "Belt and Road" zone. Areas starting from Northwest China to Western Asia and Northern Africa are faced with severe water shortages, which is consistent with the fact that this area has extensive deserts with rare precipitation and intensive evaporation. This area has the most fragile ecological systems, which means special attention should be paid to local water resources' condition and ecological environment before carrying out urban and transport planning in this area.

**Figure 2.** The spatial distribution of mean precipitation along the "Belt and Road" zone from 1985 to 2016.


**Table 2.** "Belt and Road" regional precipitation statistics table.

Furthermore, we also analyzed the inter-annual variation trend of precipitation during 24 h with at least 50 mm according to the Chinese national standard of "grade of precipitation (GB/T 28592-2012)". Figure 3 shows an increasing trend with extreme precipitation in most areas of South-Eastern Asia, South-Eastern China and eastern Southern Asia, while the change trend in other areas is not obvious.

**Figure 3.** Inter-annual variation trend of precipitation during 24 h with at least 50 mm from 1985 to 2016.

### *3.2. Droughts*

Drought is a natural disaster that has high occurrence frequency, long duration and a wide range of impacts. Drought can be regarded as regional and time-series water deficit processes, resulting in diminished water resource availability and ecosystem carrying capacity [38,39]. As shown in Figure 4a, Central Asia, Southern Asia, Western Asia and Northern Africa along the "Belt and Road" zone are all threatened by desertification and drought to varying degrees [40]. As can be seen, the drought in most areas of Western Asia and Northern Africa and Northwestern China are the most severe and the impact of drought in Mongolia-Russia and Central Asia is the most moderate. The severity of drought in eastern China, Southern Asia, South-Eastern Asia and Central and Eastern Europe is the slightest.

**Figure 4.** (**a**) drought frequency map along the "Belt and Road" zone; (**b**) drought frequency hierarchical map.

Then the drought frequency map was classified into five classes of approximately equal number of grid cells (Figure 4b). The greater the grid cell value in the final dataset, the higher the relative frequency of drought occurrence. The five drought frequency grades are low drought risk, low to medium drought risk, medium drought risk, medium to high drought risk and high drought risk respectively.

As shown in the table (Table 3), droughts in Western Asia and Northern Africa are the most serious, with high drought risk area covering 67% of the region, far exceeding the average level of "Belt and Road" zone. Droughts in China and Southern Asia are slightly better than Western Asia and Northern Africa but still not optimistic. The shares of the areas that suffer with medium drought risk and medium to high drought risk in South-Eastern Asia are both higher than the average level of "Belt and Road" zone. Central and Eastern Europe and Mongolia-Russia and Central Asia are the least affected areas by drought along the "Belt and Road" zone.


**Table 3.** "Belt and Road" drought impact area and proportion statistics table.

### *3.3. Flood*

Flood is the most common natural hazard in the world which causes tremendous losses to human life and property every year [41–43]. Based on the analysis of the distribution of floods that happened along the "Belt and Road" zone from 1985 to 2016, the flood frequency map was obtained (Figure 5a); it revealed that southern China, South-Eastern Asia and the north of Southern Asia are the regions with the most serious floods, while floods in Mongolia-Russia and Central Asia, Central and Eastern Europe and Western Asia and Northern Africa are relatively rare.

**Figure 5.** (**a**) flood frequency map along the "Belt and Road" zone; (**b**) flood frequency hierarchical map.

Then the flood frequency map was classified into five classes of approximately equal number of grid cells (Figure 5b). The greater the grid cell value in the final dataset, the higher the relative frequency of flood occurrence. The five flood frequency grades are low flood risk, low to medium flood risk, medium flood risk, medium to high flood risk and high flood risk respectively.

As shown in the table (Table 4), floods in Southern Asia are the most serious, with high flood risk area covering 67% of the region, far exceeding the average level of the "Belt and Road" zone. Floods in China and South-Eastern Asia are slightly better than Southern Asia but it is still not optimistic. The shares of the areas that suffer with medium flood risk and medium to high flood risk in Western Asia and Northern Africa and Central and Eastern Europe are both higher than the average level of "Belt and Road" zone. Mongolia-Russia and Central Asia is the least affected area by flood along the "Belt and Road" zone. There are especially few floods in Central Asia.


**Table 4.** "Belt and Road" flood impact area and proportion statistics table.

In addition, the distribution of flood in 2016 and potential affected cropland, railway and highway was analyzed. The most serious floods occurred in South-Eastern China. Southern Asia was affected by floods the most extensively. Total area affected by flood is the largest in Mongolia-Russia and Central Asia, while no floods occurred in Central Asia. Floods and disasters not only cause fatal blows to infrastructure of cities and industries, they also have serious consequences on agriculture and transportation [44]. Statistics on regions along the "Belt and Road" zone affected by floods are made based on spatial distribution of cropland, highways and railways. The results are shown below (Figure 6, Table 5).

**Figure 6.** (**a**) spatial distribution of floods in 2016 along the "Belt and Road" zone; (**b**) potential affected cropland by floods in 2016 along the "Belt and Road" zone; (**c**) potential affected railway by floods in 2016 along the "Belt and Road" zone; (**d**) potential affected highway by floods in 2016 along the "Belt and Road" zone.


**Table 5.** Potential affected cropland/railway/highway by floods in 2016 along the "Belt and Road" zone.

In 2016, Mongolia-Russia and Central Asia had the largest area affected by flood disasters, reaching as high as 2 million km<sup>2</sup> . However, since flood areas are mainly located in sparsely populated old-growth forest areas, the area of cropland affected is not large. China and Southern Asia have the largest areas of cropland affected by floods, both exceeding 1 million km<sup>2</sup> . With respect to highways and railways, China, Southern Asia and the Mongolia-Russia Central Asia are the most heavily affected areas. Central and Eastern Europe and Western Asia and Northern Africa are the least affected areas in cropland, highways and railways (Figure 7a,b).

**Figure 7.** (**a**) Statistic of floods potential affected land and cropland along the "Belt and Road" zone in 2016; (**b**) Statistic of floods potential affected railway and highway along the "Belt and Road" zone in 2016.

### **4. Hydrological Disaster Impact Analysis**

### *4.1. Water Hazard Risk Analysis*

Water security is a key element to national and social development and regional stability, and many scholars have studied the vulnerability framework which combines natural and human-related risks [45–48]. In this paper we drew the water hazard risk map by combining flood data (Figure 5b) and drought data (Figure 4b) during 1985–2016 by accumulating the assigned values of drought severity and flood severity (Figure 8). Then the severity of the water hazard risk map was classified into five classes of approximately equal number of grid cells: low risk, low to medium risk, medium risk, medium to high risk and high risk respectively.

**Figure 8.** Water hazard risk map.

According to the statistics data from water hazard risk map (Table 6), 80% of the lands in Southern Asia are threatened by high and medium to high water hazard risk, which is much higher than the average level of the "Belt and Road" zone. The areas face high and medium to high water hazard risk in China, Western Asia and Northern Africa and South-Eastern Asia have also reached 65%, 62% and 56% respectively; the security of water resource is also not optimistic. The best areas with respect to water resources condition along the "Belt and Road" zone is Mongolia-Russia and Central Asia, where the areas facing high hazard risk in water resources are basically zero. The water security condition in Central and Eastern Europe is at a moderate level.


**Table 6.** "Belt and Road" regional water hazard risk statistics table.

As shown in the Figure 9, Southern Asia faces much higher water risk than other regions, followed by China and South-Eastern Asia, Mongolia-Russia and Central Asia; Central and Eastern Europe have the best water security condition.

**Figure 9.** "Belt and Road" regional water hazard risk histogram.

### *4.2. Potential Impact Population Analysis*

The "Belt and Road" zone passes through three continents—Asia, Europe and Africa, covering a wide range of areas, with complex and diverse natural environments and highly fluctuating population density (Figure 10). Spatial distribution of population has been recognized as a fundamental indicator of various studies including ecosystem assessment [49].

**Figure 10.** "Belt and Road" zone population density map.

By integrating the population density map of areas along the "Belt and Road" zone with the map of water hazard risk (Figure 8), the results are shown as the following table (Table 7). This means that 84%, 82%, 79% and 77% of the population of China, South-Eastern Asia, Southern Asia, Western Asia and Northern Africa are facing potential high water hazard risk, which is the most severe along the "Belt and Road" zone. The population facing potential water hazard risk in Mongolia-Russia and Central Asia is the least, followed by Central and Eastern Europe.


### *4.3. Potential Impact Corridor Analysis*

The key areas of "Belt and Road" initiative mainly include six corridors, which are China–Mongolia–Russia Corridor (CMRC), New Eurasian Continental Bridge (NECB), China-Central Asia-West Asia Corridor (CCAWAC), China–Pakistan Corridor (CPC), Bangladesh–China–India–Burma Corridor (BCIBC) and China–Indo-China Peninsula Corridor (CICPC). Take NECB as an example, it runs through China and Central Asia with possible plans for expansion into South and West Asia. The Eurasian Land Bridge system is important as an overland rail link between China and Europe, with transit between the two via Central Asia and Russia. By integrating the corridors map (Figure 1) along the "Belt and Road" zone with the map of water hazard risk (Figure 8), the results are shown as the following table (Table 8). This means that 85%, 72% and 57% of the area of the China–India–Burma Corridor, China–Indo-China Peninsula Corridor and China–Pakistan Corridor are facing potential high water hazard risk, which are the most severe along the "Belt and Road" zone. China–Mongolia–Russia Corridor and New Eurasian Continental Bridge have the least area facing potential water hazard risk, followed by China-Central Asia-West Asia Corridor.

**Table 8.** Potential impact corridor area along the "Belt and Road" zone by water hazard.


### **5. Discussion**

Water hazard risk of regions around the "Belt and Road" have increased in response to continued global warming and rapid urbanization. Understanding the spatial variability of water resources state of regions is necessary to the "Belt and Road" initiative. In nations of Western Asia and Northern Africa (WANA), the major water problem is the high drought risk (67% in Table 3), which leads to a decrease in food production, forest ecosystems degradation, expansion of desert and so on. Herewith, to WANA, the steps of afforestation, forest care and management, tree species improvement, domestic water saving, water use and irrigation efficiency improvement, utilization of sewage and rainwater, industrial water recycling and sewage treatment should be included in the "Belt and Road" initiative. Whereas, in Southern Asia, South-Eastern Asia and China, the main threat of water is the high risk of flood (Table 4), which also leads to vast economic loss and damage. Herewith, various methods such as enhancing flood prevention research, strengthening hydrologic infrastructure, emphasizing the role of flood early-warning systems, raising the standard of flood control and utilizing the resources of flood water, will be helpful to these regions in the "Belt and Road" initiative.

According to the experience of China, drought can be solved by the thought of harmony between human and water, construction of water-saving society, management of water resources and strategy of connecting river and lake systems. To the flood outside the city, a large number of water projects are necessary. While loosening waterlogging in the city needs cooperation from several departments of government, including water resources, municipal administration, transportation, land and so on. Of course, many countries around the "Belt and Road" are very concerned about water problems. However, under the conditions of frequent extreme climate, lack of water resources, fragile ecological environment and complex transboundary water resources issues, water resources security and its corresponding ecological security are significant to the "Belt and Road" initiative.

In this paper, we focused on floods and droughts to represent water security, which is mainly retrieved from surface water. Whereas groundwater is important to floods and droughts for providing nearly half of the water used for irrigated agriculture and it supplies drinking water for billions of people. Groundwater levels declining will exacerbate the risk of droughts in countries around the "Belt and Road". Furthermore, the hydrological interaction between ground water and surface water will loosen the risk of floods. So, when considering the water security around the "Belt and Road", groundwater resources and their availability for exploitation should be taken into account in the future. Furthermore, apart from the quantitative aspect, water quality both of surface and ground water is also a significant issue to water security. It is the same for droughts and floods; water quality deterioration has been classified as an important water hazard. Although water quality data is scarce around the "Belt and Road", the utility value of surface and groundwater should also be analyzed in terms of its quality in future study.

### **6. Conclusions**

Analysis of water security along the "Belt and Road" zone from 1985 to 2016 shows that, (1) The areas with the most precipitation are mainly distributed in the southeast, including South-Eastern Asia, southern China and eastern Southern Asia. Precipitation is scarce from Northwestern China to Western Asia and Northern Africa, with annual precipitation being less than 100 mm in most areas. (2) To the impact of floods, Southern Asia is most serious impacted, followed by China and South-Eastern Asia; Mongolia-Russia and Central Asia are the least affected area by flood along the "Belt and Road" zone. (3) To the impact of droughts, Western Asia and Northern Africa are the most serious, followed by China and Southern Asia. Central and Eastern Europe and Mongolia-Russia and Central Asia are the least affected areas by drought along the "Belt and Road" zone. (4) To the potential water hazard risk, China, South-Eastern Asia, Southern Asia, Western Asia and Northern Africa have the largest population number facing potential high water hazard risk. China–India–Burma Corridor and China–Indo-China Peninsula Corridor have the largest areas facing potential high water hazard risk.

Water security has become a key constraint to the sustainable economic and social development of countries along the "Belt and Road" zone. Rapid urbanization has exacerbated the contradiction between water shortage and water demand and has also caused the increasingly serious floods and droughts in urban areas. Therefore, bearing capacity of water resources and the environment should be carefully considered before formulating urban and transport planning. All these considerations will contribute to the smooth implementation of the "Belt and Road" Initiative. Of course, there are still many shortcomings in this study. For example, all the countries along the "Belt and Road" zone are not included in the analysis. Only the precipitation, drought and flood in the region are analyzed. The data such as groundwater and soil moisture are not taken into consideration. The utility value of surface and groundwater should also be analyzed in terms of its quality. All of the above shortcomings will be the emphasis in our future research.

**Author Contributions:** Data curation, Z.L. and X.S.; Investigation—Methodology, M.C.; Resources, P.K.; Supervision, Y.H.; Writing—original draft, Z.L. and M.C.; Writing—review and editing, Y.H. and M.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by National Natural Science Foundation of China (Grant NO. 41822104), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA19040402), the National Key Research and Development Program of China (Grant No. 2017YFB0503005 and 2016YFC0401404) and China high-resolution earth observation system (21-Y20B01-9001-19/22).

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

**Informed Consent Statement:** Not applicable.

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

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

### **References**


### *Article* **Machine Learning to Evaluate Impacts of Flood Protection in Bangladesh, 1983–2014**

**Achut Manandhar 1,\* , Alex Fischer <sup>2</sup> , David J. Bradley 2,3,4, Mashfiqus Salehin <sup>5</sup> , M. Sirajul Islam <sup>6</sup> and Rob Hope <sup>2</sup> and David A. Clifton <sup>1</sup>**


Received: 22 December 2019; Accepted: 1 February 2020; Published: 11 February 2020

**Abstract:** Impacts of climate change adaptation strategies need to be evaluated using principled methods spanning sectors and longer time frames. We propose machine-learning approaches to study the long-term impacts of flood protection in Bangladesh. Available data include socio-economic survey and events data (death, migration, etc.) from 1983–2014. These multidecadal data, rare in their extent and quality, provide a basis for using machine-learning approaches even though the data were not collected or designed to assess the impact of the flood control investments. We test whether the embankment has affected the welfare of people over time, benefiting those living inside more than those living outside. Machine-learning approaches enable learning patterns in data to help discriminate between two groups: here households living inside vs. outside. They also help identify the most informative indicators of discrimination and provide robust metrics to evaluate the quality of the model. Overall, we find no significant difference between inside/outside populations based on welfare, migration, or mortality indicators. However, we note a significant difference in inward/outward movement with respect to the embankment. While certain data gaps and spatial heterogeneity in sampled populations suggest caution in any conclusive interpretation of the flood protection infrastructure, we do not see higher benefits accruing to those living with higher levels of protection. This has implications for Bangladesh's planning for future and more extreme climate futures, including the national Delta Plan, and global investments in climate resilient infrastructure to create positive social impacts.

**Keywords:** Bangladesh; climate resilience; flood protection; machine learning; socio-environmental impacts

### **1. Introduction**

Climate change is expected to increase the frequency and extent of extreme flood events, which will directly impact the environment and the livelihoods of people in the affected areas [1,2]. Low-lying coastal regions of the world are particularly vulnerable to these flood events and sea level rise [3,4]. The issue is compounded in countries such as Bangladesh, where about 60% of the country is lower than 6 m above the sea level [5] while more than 70% of land is used for agriculture, the country's primary

economic source [6]. Although there have been global investments on climate change adaptation, these adaptation measures can have both beneficial and unintended detrimental consequences when wider issues or longer time frames are considered [7,8]. Integration of adaptation actions and policies across sectors at different spatio-temporal and societal scales remains a key challenge to achieve effective adaptation in practice [9]. Another major challenge is the lack of consistent empirical methods linking climate change to the impacts on the environment and the livelihoods of people [10]. Although success of these interventions relies on developing principled methods to monitor and evaluate the impacts across different environmental and socio-economic factors [11], currently there is a lack of such methods in the literature. This work aims to bridge a key gap in knowledge by proposing rigorous analytical methods to evaluate the impacts of adaptation measures.

Bangladesh has had five decades of political and policy attention focused on implementing flood mitigation strategies, whose evaluation has been documented in the past literature on Bangladesh flood protection infrastructure. The total flood protection coverage area currently stands at 5.37 million ha, more than one third of the country [12,13]. Primarily aimed to protect monsoon crops and prevent damage to homesteads, these interventions have often not considered the social, economic, and environmental dimensions of water resource management [13]. While the interventions led to several positive impacts, they also resulted in considerable medium to long-term negative consequences in many places. The flood secure environment facilitated enhanced economic activities in the protected areas, e.g., increased agricultural output [14]. But when embankments failed to provide protection during moderate to extreme floods, the resulting damage was higher owing to the accelerated economic activities compared to areas outside embankment [13,15]. Floodplains were deprived of several environmental and ecological functions, including improvement of soil fertility from silt-laden inundation water, groundwater recharge, and biodiversity. The disruption in hydraulic connection between river and floodplain led to substantial damage to fisheries and local boat transports [14–18], compromising the livelihood activities of people, especially the marginalized.

Previous evaluations of flood protection investments in Bangladesh have widely suggested that it has been difficult to attain the stated objectives of the interventions based on only technical and economic viability, but without giving due consideration to the hydromorphological features of the floodplain and the socio-economic condition of its inhabitants. Despite the negative consequences of embankments, the popular demand for flood protection by people has been high. Hence, the priority of the Government started to shift from traditional flood control to flood management towards the later stages of the Flood Action Plan (FAP) in late 1980s. Here, flood control refers to the conventional method of constructing an embankment and drainage regulator whereas flood management refers to mitigating flood damage without causing degradation of the floodplain environment, which might involve implementing floodplain land use regulation that identifies floodplain zones and enforces appropriate planning and design during construction of infrastructures in these floodplain zones to account for flood and preservation of floodplain resources and environment. There was a real paradigm shift to integrated flood management, i.e., covering issues relevant to not only flood but also drainage, irrigation, navigation, environment and socio-economic development, which was subsequently reflected in the National Water Policy (NWP) in 1999 [19] and the National Water Management Plan (NWMP) [12]. A combination of structural and non-structural measures was envisioned, with full structural protection against floods in regions of economic importance (such as metropolitan areas), a reasonable degree of protection in other critical areas (such as district towns), and flood proofing measures in the rural areas. However, translating integrated management into action, particularly at the program and project levels, has remained a major challenge [20].

The Bangladesh Delta Plan 2100 also gives more emphasis to restoration, redesign and modification of existing embankments and associated structures, and a high importance to the urgency of maintenance [21]. New developments are envisioned only for protecting economic strongholds and critical infrastructure. Most are already in place but requires improvement; additional developments will be required in some areas. In order to properly evaluate the impacts of these investments, this work strongly recommends adoption of principled analytical methods early on, considering socio-economy and environment across longer time frames, so that appropriate studies are conducted at the outset, e.g., baseline and periodic survey and monitoring to assess the impact of interventions.

The changes in Bangladesh Delta are affected by many factors, upstream interventions being one of them. In this work, we evaluate the impact of flood infrastructure on a large project in Bangladesh. We took advantage of an existing data set of this project to further clarify evaluation outcomes in Bangladesh flood management using machine learning approaches. We chose the site also because it had robust historic continuous data to trial machine learning that enabled the type of analysis that other sites could not. Machine learning, generally considered a subset of artificial intelligence, is a field of study that uses algorithms and statistical models to learn patterns from data so that useful inference may be made about new data. Machine-learning approaches are useful in the context of this study because they provide robust metrics to evaluate the impacts of interventions as well as identifying the most informative indicators.

### *Context and Related Works*

The Meghna–Dhonagoda Irrigation Project (MDIP) is one of the largest flood control, drainage and irrigation projects in Bangladesh, implemented in 1988 with the objective to protect the area in Matlab North from river flooding and drainage congestion during the monsoon via embankment and regulators. The primary aim was to improve agricultural conditions in monsoon, with special reference to encouraging introduction to high yielding variety (HYV) monsoon rice (Aman), and also to provide irrigation from surface water in the Rabi (dry) season and early monsoon seasons. Located 55 km south-east of the capital city Dhaka, the study area is 184.4 km<sup>2</sup> with a population of 230,185 as of 30 June 2014 [22]. The Dhonagoda river bisects the area into Matlab North and Matlab South. The embankment protects the Matlab North area from flooding from Meghna river on the north and west and Dhonagoda River on the east and south (Figure 1).

**Figure 1.** Baris (cluster of households) are colored based on whether they are located inside (brown) or outside (blue) the embankment (red).

This study area provides a unique setting to evaluate the long-term impacts of an embankment on households living on both sides of the embankment by using an unrelated but population-wide Health and Demographic Surveillance System (HDSS) [22,23] managed by the International Centre on Diarrheal Disease Research, Bangladesh (icddr,b). The research site was created in 1963, over 25 years before the embankment's construction, to provide field-based research on cholera vaccines and treatments. The site maintains a primary research focus on studying public health interventions and demographic changes. The study site evolved around two sets of comparison zones. The primary division is a quasi-experimental maternal and child health and family planning program where half the population received reproductive health interventions and home visits from community health workers, while the other half received the government's public health services [22].

The other key division across the study area is between those within and those outside the embankment. The embankment-led agricultural and related economic variables are not rigorously tracked through the icddr,b data systems and health-focused research agenda [24]. Discussion of the embankment's impact did not appear in the icddr,b annual reports until the 2014 analysis of the socio-economic survey and is not a feature of these reports which serve as a core site output [23]. However, the embankment provided a clear division of the population between people living inside vs. outside the embankment [22,25]. The HDSS data, spanning pre–post-embankment periods provide a unique opportunity to study the long-term impact of embankment.

Multiple studies have used the HDSS data to advance analysis around the indirect impacts of embankment [25–30]. Earlier studies of child mortality in Matlab in 1981 identified the socio-economic factors as indicators of severe malnourishment, impacting child mortality, which was a basis for socio-economic survey focus on household assets, education and occupation [26]. A later study during 1983–1992 showed that child mortality was higher outside than inside the embankment, with the differences particularly significant for deaths caused by infectious diseases [27]. The study recommended comparison of long-term mortality data with migration patterns and other factors (e.g., proximity to main rivers). Another study of the embankment's impact on cholera during 1983–2003 showed that after controlling all other environmental variables, living inside the embankment area does not appear to affect whether the household experiences cholera. Counter-intuitively, among households where cholera is reported, living inside the embankment area significantly increases the number of cholera cases, possibly due to a combination of environmental factors and behavioral change [29]. In particular, a joint study by Bangladesh Rural Advancement Committee (BRAC) and International Centre for Diarrheal Diseases Research, Bangladesh (icddr,b) analyzed the impact of embankment on both environment and people [28]. Based on both quantitative and qualitative surveys conducted in 1992 and 1996 respectively, the study revealed both positive and negative impacts. The positive impacts were associated with a higher level of agricultural yields and economic prosperity, particularly for a subset of farmers, while the negative impacts were associated with lower fish catch and intake of fruit and vegetables, displacement of poorer households, and complaints about ill health due to greater demand for agricultural labor [28]. However, the core findings around welfare differences were inconclusive. The study recommended doing a follow-up study using longer-term data to analyze how these impacts would evolve over time. Another study found that the embankment resulted in a net welfare loss, which had been an outcome of higher than anticipated construction costs, lower benefits to agriculture due to loss of soil fertility over time, higher waterlogging damages and the highly negative impact on capture fisheries [30]. The embankment also caused negative distributional outcome, with big landowners benefiting from increased agricultural production, reduced property damage, and increased livestock and aquaculture production. In contrast, the traditional fishermen and the river transport workers, belonging to the poor sections of the community, suffered significantly negative impact. However, the variables related to these impacts were not monitored consistently over time, limiting ability to include in longer-term impact assessments.

Most studies occurred within a decade of the embankment completion, showing only the short-term impacts. There is a limited body of literature exploring long-term impacts outside agricultural productivity, providing a core motivation of this paper to re-purpose multi-decadal HDSS data. Using this data, machine-learning approaches can be used to empirically evaluate the impact of embankment and identify the most informative socio-economic indicators. Besides, these approaches provide robust metrics to evaluate the model's accuracy and generalizability to future examples. We explore the extent to which machine learning can provide analytical insights from detailed historical data on long-term impacts from climate resilient infrastructure to help guide future policy and investments.

### **2. Methods**

### *2.1. Machine-Learning Approaches*

Machine-learning approaches have found useful applications in a wide variety of fields, including text processing, computer vision, healthcare, finance, and robotics [31–36]. Recently, these approaches have also been applied to socio-economic [37–40] and environmental [41–45] studies.

We implement machine learning on two types of data to answer two specific questions:


To answer the first question, socio-economic variables corresponding to a household (inputs) are mapped to a binary label (output). The responses provided by households during a socio-economic survey are used to determine whether the households live inside or outside embankment. By comparing classification outputs over time, we can infer whether the embankment has caused differences in welfare inside vs. outside embankment over time. To answer the second question, a regression model is used to map time (input) to event rate (output), e.g., mortality rate per year. After learning two independent regression models for two event rates corresponding to inside vs. outside embankment, we can analyze whether the embankment has caused differences in the event rates over time. An array of approaches exists to perform classification and regression. The ones adopted in this work are motivated by their suitability to model the available data in answering the aforementioned questions.

### 2.1.1. Classification Approaches

For each household, its label *y<sup>n</sup>* (output), indicating whether it lies inside or outside embankment, is defined in terms of its responses to *D* survey questions **x<sup>n</sup>** = [*xn*1, ..., *xnD*] *T* (inputs) as follows:

$$y\_n = f(\mathbf{x\_n}) + \varepsilon\_{n\nu} \tag{1}$$

where *e<sup>n</sup>* is the residual or noise, and *f*() is the mapping function specific to the type of classifier. The mapping function and its parameters can be learned using a collection of household survey responses and their corresponding labels (termed "training" examples). Once a mapping function is learned, it can be used to predict labels for new previously unseen (termed "test") examples. By evaluating whether the labels are correctly predicted for test examples, we can determine the ability of a classifier in discriminating two classes. The simplest mapping function is a linear function of the inputs as follows

$$y\_n = w\_0 + w\_1 \mathbf{x}\_{n1} + \dots + w\_D \mathbf{x}\_{nD'} \tag{2}$$

where the weights **w** = [*w*0, ..., *wD*] *<sup>T</sup>* are known as the parameters of the mapping function. However, this simple linear discriminant function [31] cannot model interactions between variables or other

complex phenomena. Generally, classification techniques learn a non-trivial mapping function, capable of performing nonlinear classification, which are more relevant to model real-world examples.

Some standard classification approaches are Logistic Regression (LR) [31], kernel-based methods, e.g., Support Vector Machines (SVMs) [46], decision trees, e.g., Random Forest (RF) [47], and Neural Networks (NN), e.g., Stacked Auto-Encoders (SAE) [36]. Logistic Regression (LR) is one of the most popular classification approaches, which uses a logistic sigmoid function to perform probabilistic binary classification. However, without using kernels, LR is often only suitable to classify linearly separable examples. In general, we expect other above-mentioned approaches to perform as well as or better than LR.

Random Forest (RF) is a type of ensemble-based method that performs an average over an ensemble of many estimates obtained over bootstrapped subsets of data, where the ensemble of estimators can be thought of as leaves of a tree [47]. An added advantage of RF is that it also ranks the input variables by their importance to discriminate the two classes. In the context of this study, the variable importance ranking helps identify which survey questions are most informative of households living inside vs. outside embankment. When using RF, it is important to choose an appropriate number of estimators, and optimize other parameters.

Stacked Auto-Encoders (SAE) is a type of deep learning methods that first learns the lower-dimensional representation of data by constraining the hidden layers to capture the most relevant aspects of the data. Then the whole network is discriminatively fine-tuned like a feedforward neural network to perform classification [36]. Similar to Principal Component Analysis (PCA) [34], which is the most common dimension reduction technique, SAE can be used to learn useful lower-dimensional representations of data to uncover patterns in data, e.g., identify clusters of households with similar features. Unlike PCA, SAE is not limited by the Gaussian distribution assumptions of the input space, it is more flexible to model mixed data types (binary, categorical, and continuous), and being a deep learning approach, it is capable of handling large amounts of data. These are all favorable properties for analysis of socio-economic survey data. It is important to optimize SAE's parameters such as learning rate, batch size, number of epochs, number of hidden nodes, number of hidden layers, etc.

We do not go into details regarding the specifics of these approaches other than provide a general intuition for the ones that are implemented. Having experimented with a variety of these approaches, we report majority of our results based on RF, while comparing RF to LR and SAE on one specific example for reference.

### 2.1.2. Regression Approaches

Shifting from discrete to continuous output labels, for each year *x<sup>n</sup>* (input), the event rate for that year, i.e., its label *y<sup>n</sup>* (output), can be defined with (1), where *y<sup>n</sup>* is now a continuous variable. Different regression approaches exist, e.g., Relevance Vector Machines (RVMs) [48], Gaussian Processes (GPs) [49]. Gaussian Processes (GPs) are probabilistic methods (i.e., capable of giving uncertainty of model's predictions) that can be used to model time-series data as a distribution over function [49]. We use GPs because in addition to modeling event rates as time-series, they are also useful to compare differences in the resulting time-series models. In a Bayesian framework, the event rates **y** = [*y*1, ..., *yN*] for N years **x** = [*x*1, ..., *xN*] can be defined by a GP prior as follows

$$f(\mathbf{x}) \sim \text{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}^T)),\tag{3}$$

where **y** = *f*(**x**) + *e*, *e* is the noise, *m*(**x**) is the mean function and *k*(*x*, *x T* ) is the kernel or covariance function. This allows the posterior distribution of the function evaluated at a finite set of points **x**<sup>∗</sup> to be a multivariate Gaussian distribution as follows

$$p(\mathbf{f}\_\*|\mathbf{x}\_\*, \mathbf{x}\_\*\mathbf{y}) = N(\mathbf{f}\_\*|\mu\_\*, \mathbf{\Sigma}\_\*),\tag{4}$$

where *µ*<sup>∗</sup> and **Σ**<sup>∗</sup> are the posterior mean and covariance, respectively. This posterior formulation allows us to compare the differences in two event rates in a principled manner, which is described in the next section. When using GPs, it is important to choose an appropriate kernel function and optimize its parameters.

### *2.2. Evaluation Metric*

We evaluate inside/outside embankment classification using Receiver Operating Characteristics (ROC) curve based on k-fold (k=3) cross-validation [31], which prevents a classifier from overfitting to training examples, thus ensuring the results are generalizable. Intuitively, a diagonal ROC curve means the predictors (socio-economic survey variables) are not indicative of inside/outside class discrimination. On the other hand, the more the curve pushes towards the top-left corner, the more discriminatory are the predictors. A statistical significance test can be performed to compare the area under ROC curves (AUCs) of the later three years with the AUC of 1982 [50]. Whenever a comparison is statistically significantly different with *p*-value less than 0.01, the corresponding *p*-value is highlighted with a ∗ symbol.

In line with the previous studies [25,29], to compare differences in events data pre- vs. post-embankment, we divided the 32-year study duration into pre (1983–1989) vs. post (1990–2014) embankment periods. Within each of these two periods, we can compare the temporal differences in these events using a Gaussian Process-based Bayesian statistical significance test [51]. The posterior formulation in (4) allows us to model the differences in event rates, ∆**f**<sup>∗</sup> = **f1**<sup>∗</sup> − **f2**∗, as another multivariate Gaussian with mean ∆*µ*<sup>∗</sup> = *µ*1<sup>∗</sup> − *µ*2∗, and covariance ∆**Σ**<sup>∗</sup> = **Σ1**<sup>∗</sup> + **Σ2**∗. Then, we say the two event rates are equal with posterior probability 1 − *α* if the credible region for ∆**f**<sup>∗</sup> includes the zero vector or, in other words, if

$$(\Delta \mu\_\*)^T (\Delta \Sigma\_\*)^{-1} \Delta \mu\_\* \le \chi\_\nu^2 (1 - \mathfrak{a})\_\prime \tag{5}$$

where *χ* 2 *ν* (1 − *α*) is the 1 − *α*-quantile of a Chi-squared distribution with *ν* degrees of freedoms and *ν* is the number of positive eigenvalues of ∆**Σ**<sup>∗</sup> [51]. This test can be used to compare two event rates within each (pre or post) embankment period.

### **3. Data**

Part of the HDSS dataset is available for this study, covering roughly one third of the study area geographically.

### *3.1. Socio-Economic Survey Data*

Socio-economic survey data are available for four years, roughly a decade apart. The number of questions in surveys increased over the years from 21, 40, 62, to 146 in 1982, 1996, 2005, and 2014, covering 14791, 19448, 22799, and 25840 households, respectively. Considering 1982 as a baseline pre-embankment period resulted in only four variables equivalent across all four years. Those variables were agricultural land ownership, primary drinking water source, number of cow/buffaloes/goats owned, and boat ownership. Considering only the later three years (1996, 2005, and 2014) resulted in slightly more (12) variables equivalent across all three years. Those variables were agricultural land ownership, homestead land ownership, primary drinking water source, number of cow/buffaloes/goats owned, boat ownership, household assets—sofa, chair/table, showcase, radio, TV, bike/bicycle, primary rood structure, and sanitation facility type. For a consistent comparison, we only consider households that are common across all years, i.e., 10563 and 14276 households common across all four and the later three years, respectively. We note that the results are similar when all households per year are used. For a separate analysis specific to socio-economic survey data from 2014, whose aim was not to compare with prior years, we use all available variables and households.

### *3.2. Events Data*

Events data, collected monthly to bi-monthly, correspond to birth, marriage, divorce, death, inside/outside migration, and inward/outward movement [23]. Here, *migration* refers to the migration from anywhere outside the study area into either the embanked part or the outside part of the study area; whereas *movement* refers to the movement of the study area inhabitants either into or out of the embanked area. We hypothesize that if the flood-protected area were to provide increased socio-economic benefits and stability, people would more likely move to inside the flood-protected area. Although experiments were performed with all data, we only provide results from analyzing death, internal movement, and migration data, which are presumed to be the most informative of the effect of embankment. The event counts in each group (inside/outside) are normalized by the inside/outside annual mid-year population. Since the mid-year population was unavailable for one particular year, we only use data from the remaining 31 years.

### **4. Results**

### *4.1. Socio-Economic Survey Data Analysis*

ROC comparison across four years in Figure 2a seems to suggest there are differences in inside/outside discrimination over time based on socio-economic variables. However, one of these variables, the ownership of a boat, is mostly a consequence of differences in habitat due to embankment rather than an indication of welfare. Indeed, when we remove the ownership of boat variable, Figure 2b shows that the discrimination does not change much over time except for 1996. To understand the increased discrimination in 1996, we rank the three variables (agricultural land ownership, primary drinking water source, and number of cow/buffaloes/goats owned) by their importance in Figure 2c. We observe that agricultural land ownership is the most discriminative variable across all four years. When we classify households based on only agricultural land ownership (results not shown), we observe that although the discrimination increases temporarily during 1996, it falls back to pre-embankment period over time. The statistical significance test in Table 1 supports that only the 1996 AUC metric is significantly different compared to 1982. The results suggest that although there was a short-term increase in agricultural land ownership for inside residents, the difference seems to have evened out over time. Despite the statistical significance test, it should be noted that all AUCs, including 1996, have low values, implying that the equivalent variables are not very discriminative. The results also highlight the need for more data on agricultural productivity and markets to further investigate the differential impact over time.

Setting aside the temporal comparison for a moment, if we were to make use of more variables for a particular year, we could learn a better model and use the classifier's outputs to further investigate the most informative variables in details. Figure 2d shows that when 122 relevant out of 146 variables from year 2014 are used, the classifier performs strong inside/outside discrimination. Figure 2e shows the top 10 informative variables, ordered by their importance—homestead land size, number of households sharing drinking water source, number of chicken/ducks, agricultural land size, fuel source is wood/gas, primary income source, primary drinking water source is Arsenic-safe tubewell/pipe, and if zakat received. Each of these variables can be visualized by aggregating the corresponding household responses at bari-level and comparing each plot to inside/outside locations of baris in Figure 1 as a reference. Figure 3i,ii show that indeed households inside the embankment appear to own more homestead land and owning homestead land appears to be correlated with owning agricultural land. Consequently, land ownership is identified as the most important discriminative indicator. However, for several other variables, the difference appears to be a consequence of the households' proximity to significant infrastructures (e.g., hospital, gasline, pipeline, deep tubewells) instead of a direct consequence of the households living inside/outside embankment. Availability of gasline (Figure 3iv) in certain parts of the study area is complementary to those households not using the more ubiquitous wood as fuel source (Figure 3iii). Likewise, availability of pipeline (Figure 3vi) is

complementary to those households not using the more ubiquitous tubewell as drinking water source (Figure 3v). Similarly, larger number of households sharing a drinking water source (Figure 3viii) appears to be complementary to those households using deep tubewells (Figure 3vii).

**Figure 2.** Inside vs. outside embankment classification outputs based on socio-economic survey data using (**a**) four equivalent variables across four years, (**b**) same as (**a**) excluding boat ownership variable, (**d**) 122 variables from 2014 survey, (**c**) top three, and (**e**) top ten discriminating variables identified by Random Forest classifier in respectively (**b**) and (**d**). (**f**) lower-dimensional representation of SES 2014 data, and (**g**) the location of households, corresponding to the cluster inside the blue circle, in the study area.

**Figure 3.** Selected variables from socio-economic survey data. Responses aggregated over bari. (**i,ii**) Homestead vs. agricultural land owned, (**iii,iv**) primary fuel source is wood/wood dust/paddy husk vs. gasline, (**v–vii**) primary drinking water source is green shallow tubewell, pipeline, vs. deep tubewell, and (**viii**) number of households sharing drinking water source.

**Table 1.** Comparison of AUCs across four years based on Wilcoxon statistic. The ∗ symbol denotes the AUCs are significantly different with *p*-value less than 0.01.


Referring back to Figure 2d, we note that we expect Stacked Auto-Encoders to perform as well or better than Random Forest when Stacked Auto-Encoders' parameters are exhaustively optimized. Rather than obtaining the best discrimination, the primary goal of using Stacked Auto-Encoders was to show its value in learning useful lower-dimensional representation of data. Figure 2f shows a 3-dimensional representation of 122 survey questions, which identifies a cluster of outside residents who are very different from inside residents based on their responses to the survey questions. A geospatial plot (Figure 2g) reveals that this cluster of outside households, in fact, are all located close to the icddr,b hospital, which incidentally, is associated with factors such as availability of gasline. These results suggest that although there was a significant relationship between introduction of embankment and change in land ownership patterns, recent development progress is linked to other variables (e.g., electricity/energy and proximity to services). Spatial analysis of these variables might provide useful insight to track differential progress in the study area.

### *4.2. Events Data Analysis*

Figure 4 summarizes the results of events data analysis for three cases—inward/outward movement, inside/outside in-migration, and inside/outside mortality over time. The first column shows the individual event rate data modeled by independent Gaussian Processes. For each of these plots, the dots represent the normalized event rates per year, the solid lines represent the estimated means, and the shaded regions represent the 99% confidence interval. For each of these three cases, the second column shows the estimated differences in two event rates. The more the two event rates differ, the more the curve deviates from the zero-line, represented by the dotted black line. Finally, the dashed red line marks the pre- vs. post-embankment periods. For each of the second-row plots, within each embankment period, the confidence interval region helps us visually determine if the difference in event rates is statistically significant. If the confidence interval excludes the zero-line, the difference in event rates in this period is statistically significant. These test results are also summarized in Table 2.

Figure 4a,b show that the inward movement rate increased leading up to the embankment's construction, peaking at 1988, and continued to remain significantly higher than the outward movement rate. In recent years, the difference seems to have evened out. Figure 4c,d show that after embankment's construction, the in-migration rate to the area outside of embankment has been increasing compared to inside. However, the difference in in-migration rates is not statistically significant. Based on available data, it is unclear if the decreasing inward movement and the increasing in-migration inside the embankment area are due to growing competition for land inside the embankment or due to falling motivation to move inside. In contrast, Figure 4f shows that the embankment has not caused differences in mortality inside vs. outside over time. The mortality data was disaggregated into vulnerable (under 5-year-olds or over 70-year-olds) vs. non-vulnerable population, and male vs. female. The results (not shown) are similar to Figure 4e,f. Similar analysis done using out-migration data (not shown) shows no difference inside vs. outside. These results suggest that apart from internal movement within the study area, significant differences are not observed in migration and mortality patterns inside vs outside over time.

### *4.3. Hydro-Climatic Data Analysis*

To fully understand the impact of embankment, it is important to link the socio-economic impacts with the hydro-climatic events. With about 47% of the land being low lying, the Meghna–Dhonagoda project area within Matlab North used to be regularly flooded in the monsoon up to a depth of 2–3 m in the pre-project condition [52] (Saleh et al. 2000) (Figure 5b). In contrast, there is a lower proportion of low-lying land area in Matlab South (Figure 5b). While the low-lying areas get inundated during an average annual flood, the relatively higher lands are inundated only during moderate to extreme floods. Since Matlab South is only exposed to the Dhonagoda river, the area is relatively less vulnerable compared to other areas which are directly affected by floodwater from the large Meghna River.

Figures 5c,d show the inundated areas in Matlab in two major flood years, delineated from LANDSAT 4-5 TM images using Normalized Difference Water Index (NDWI) [53]. Although designed based on 1 in 100-year flood level, the embankment, suffering major breaches, could not protect the project area during the 1988 flood (corresponding to a 30-year flood). After repairs, the embankment successfully withstood floods during subsequent years, but frequently suffered many problems. Although the river water level was higher than the design 100-year flood level during the 1998 flood, the severest in recent history both in terms of magnitude and duration, the embankment was able to protect the project area from inundation (Figure 5d) (Saleh et al. 2002; [52]). The flood did; however, cause substantial damages to the embankment. No inundation due to riverine flood inside the embankment has been reported in later years. However, some areas inside the embankment are subject to waterlogging induced by rainfall [30,54]. Two pumping stations inside the project area are operated during the monsoon to drain out the accumulated rainwater.

Figure 5b shows that the inward movement rate behaves like an impulse function during the embankment's construction, propelling a substantial inward movement. However, subsequent spikes in the inward movement rates and the overall trend do not appear to be correlated with the monthly tidal data observed at a monitoring site. Further analysis of relevant data (e.g., financial loss due to floods) may allow us to study the relationship between hydro-climatic events and socio-economic behavior.

**Table 2.** Comparison of event rates within each (pre/post) embankment period based on GP-based. Here, *χ* 2 (*ν* = 7, p = 0.01) = 1.24, and *χ* 2 (*ν* = 24, p = 0.01) = 10.86. The ∗ symbol denotes the event rates are significantly different with *p*-value less than 0.01.


**Figure 4.** Analysis of temporal difference for three specific events— (**a**,**b**) inward/outward movement, (**c**,**d**) inside/outside in-migration, and (**e**,**f**) inside/outside mortality during 1983–2014. For each case, the left column shows the individual time-series modeled by GP, and the right column shows the estimated differences in two time-series. The red dashed line marks the pre- and post-embankment periods.

**Figure 5.** (**a**) Land elevation map of Matlab North, Matlab South and Daudkandi Upazilas based on 20 m × 200 m National DEM data (Source of data: WARPO), (**b,c**) inundated area in Matlab North, Matlab South and Daudkandi Upazilas in two major flood years [delineated from LANDSAT 4-5 TM images using Normalized Difference Water Index (NDWI)], and (**d**) (top) monthly inward/outward movement of households within Matlab, and (bottom) monthly high-tide observed at station located in (**a**).

### **5. Discussion**

The study has limitations that are largely due to either unavailability of relevant data or some inconsistency in the data collected over multiple decades. Data from only one third (geographically) of the study area were available for analysis. The study may not generalize due to small study size and differences in hydroclimatic conditions. The results based on the classification approaches should be accepted with caution due to lack of enough common indicators across different years, especially the pre-embankment baseline year. We acknowledge other infrastructures (health, roads, electricity, education, etc.) have impacts in relation to the embankment. The spatial plots of a selection of socio-economic variables suggest recent social developments may be tied to these and other interventions. Similarly, data on agricultural yields, fishing, and other water-dependent economic activities would provide vital information. The limitations in data highlight the challenges of evaluating

long-term impacts of interventions such as the embankment, where relevant indicators may not be identified during the baseline study (if performed), and may not be consistently monitored over time.

Although other health event outcome related data were collected as part of HDSS, such data were unavailable for this study. The available data only has ICD-9 (until 2001) and ICD-10 (post 2001) codes, but otherwise no specific details on identifying water-borne disease morbidity data. A medical expert's help was sought to use the ICD codes to narrow the causes potentially related to water-borne diseases. Based on this indirect method of obtaining water-borne disease morbidity data, similar analysis was performed as was done using all-cause mortality data in Figure 4e,f. The results were similar and showed no significant difference in inside/outside population. However, we note that this indirect method of obtaining water-borne disease morbidity data has its own limitation in addition to inconsistent coding schemes before and after 2001. If further health-outcome specific data were available in future, including records related to water-borne diseases or nutrition in children, we could perform a focused comparison of morbidity rates corresponding to water-borne diseases or malnutrition.

Machine-learning approaches provide further opportunities to integrate socio-economic data with other types of ancillary data, e.g., investment in public/private water infrastructure, data on fishing and farming productivity, data on loan amounts and types, hydro-climatic data, etc. Inclusion of these datasets combined with spatio-temporal analysis would constitute an interesting work to further study the embankment's impact on the environment and people. By reformulating the problem and restructuring the data, the machine-learning approaches and the framework implemented in this study can be used to answer questions relevant to related studies, e.g., identify the most informative variables indicative of owning a deep tubewell, using solar panels for electricity, contracting certain water-borne diseases, etc., which can provide valuable information to future developmental interventions.

### **6. Conclusions**

Overall, the available socio-economic indicators and mortality, migration data do not provide a strong predictive value of the location of inside vs. outside residents. The study reinforces findings in the past literature around the immediate impacts of the embankment but does not find those are continued three decades later, with certain well-being indicators evening out across the two defined areas. The proposed approaches are particularly suitable in the face of large numbers of variables and samples. These methods are applicable to evaluate the environmental and socio-economic impacts of human interventions in general beyond the context of Bangladesh.

Moreover, the proposed approaches provide a new framework to identify indicators that are relevant to evaluate the impacts of interventions. Results show households inside the embankment owned a larger proportion of agricultural land within a decade of the embankment's construction. Similarly, the embankment is associated with the differential movement of people, with more people moving inward vs. outward within the study area. However, the difference appears to be evening out in recent years.

This work performs a quantitative analysis of the impacts of embankment using machine-learning approaches. By providing rigorous analytical tools, these approaches may be relevant to tackle the global challenges of flood risks. Such principled analysis is key towards evaluating both short-term and long-term as well as intended and unintended consequences of interventions, providing evidence to support future actions and policies related to combating climate change.

**Author Contributions:** Conceptualization, A.M., A.F., D.J.B., R.H.; methodology, A.M., A.F., D.J.B., R.H., D.A.C.; software, A.M.; validation, A.M., A.F., D.J.B., D.A.C.; formal analysis, A.M., A.F., D.J.B., R.H.; investigation, A.M., A.F., D.J.B., M.S., R.H.; resources, M.S.I., R.H., D.A.C.; data curation, A.M., A.F.; writing–original draft preparation, A.M., A.F., M.S.; writing–review and editing, A.M., A.F., D.J.B., M.S., M.S.I., R.H.,D.A.C.; visualization, A.M., D.J.B., M.S.; supervision, D.J.B., R.H., D.A.C.; funding acquisition, R.H., D.A.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This document is an output from the REACH program funded by UK Aid from the UK Department for International Development (DFID) for the benefit of developing countries (Aries Code 201880). However, the views expressed and information contained in it are not necessarily those of or endorsed by DFID, which can accept no responsibility for such views or information or for any reliance placed on them.

**Acknowledgments:** The HDSS data was kindly provided by icddr,b. icddr,b acknowledges its donors who provide unrestricted support to icddr,b for its operation and research. Current donors include the Government of the People's Republic of Bangladesh, Global Affairs Canada(GAC), The Swedish International Development Cooperative Agency (Sida), and the Department of International Development (UK Aid).

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

### **Abbreviations**

The following abbreviations are used in this manuscript:


AUC Area Under ROC Curve

### **References**


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